Skip to content

Uddeepta Deka / Updated Charged Lens

Go to a project
Commit 14c793f3 by Uddeepta Deka
Options

final edits

parent 16de5cd2
Inline Side-by-side
Showing with 604 additions and 11 deletions
notebooks/Fw_computation.ipynb
This source diff could not be displayed because it is too large. You can view the blob instead.
notebooks/deflection_angle.ipynb
...@@ -28,7 +28,7 @@ ...@@ -28,7 +28,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 2, "execution_count": 7,
"id": "847a11cd-1be9-4800-b631-59046e99fb79", "id": "847a11cd-1be9-4800-b631-59046e99fb79",
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
...@@ -39,7 +39,7 @@ ...@@ -39,7 +39,7 @@
"import warnings\n", "import warnings\n",
"from plotsettings import *\n", "from plotsettings import *\n",
"from cycler import cycler\n", "from cycler import cycler\n",
"style = 'Notebook'\n", "style = 'APS'\n",
"plotsettings(style)\n", "plotsettings(style)\n",
"plt.rcParams['axes.prop_cycle'] = cycler(color=dark2)\n", "plt.rcParams['axes.prop_cycle'] = cycler(color=dark2)\n",
"warnings.filterwarnings('ignore')" "warnings.filterwarnings('ignore')"
...@@ -186,8 +186,69 @@ ...@@ -186,8 +186,69 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 8,
"id": "a2c4720f-aa20-41d1-b0f2-dbeb6ef4feb9",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAUsAAADiCAYAAADH0+UMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABVx0lEQVR4nO3deVSc6X7Y+W/tG9TKvlMFYhNaQOp9u33hdt/FfW03tGLH8ThxGiXxJD6TScTck/HxTWbpgyaZjJNxjlE7ju1k7Kjh2r6L7wa67tt7twTaFwQUW7FDLRRF7fXOH4gySCCVEJvE8zmH000t7/u8b6l+POvvkUmSJCEIgiDcl3y3CyAIgvA4EMFSEAQhBSJYCoIgpEAES0EQhBSIYCkIgpACESwFQRBSIIKlIAhCCkSwFIQnmNfrxev17nYxnggiWD5muru7sVgsOBwOWltbaW1tpbm5md7e3pSPcfr0aU6ePIlMJqO5ufmRy+RwOPbUF7K3t5eTJ08+1D2B5fvS3NxMa2vrlpWlu7ubkydPbtnxHkZvby+lpaVcuHDhnsc3c3/2PUl47DQ0NEgtLS3J3wcHByWz2Sy1t7c/8L1tbW3J17W1tUlNTU2bKkNPT8+6/79X2O12qaurK+XXt7W1SR0dHZIkLd/PrTI4OCjt5tesrq5u3fvwsPdHkCRRs3xMmc3m5P/b7XbeffddTp48+cAa3vnz57Hb7cn3bYbT6VxTI62rq9vUcbbT6vuTivPnzyffs9n7sh6r1bplx9pKD3t/BNEMf2I0NTUBy80+gDNnznD69GkaGxuTza3u7m56e3vp6OjgzJkznD17ds0x1nsPQGdnJ2fOnKG1tZUzZ87Q29uL0+nkzJkzvP/++9TX1ydff+bMmeRxVpqf3d3dOBwOOjs7qa+vp7GxccPraG1tpbu7m+bmZpxO5wPf393dzenTp6mvr8disXDmzJl7jrnRda0+xsp96ezs3PA6Vs7f2dmJxWK55w/T3ffp7ucaGxvXlH29a13vHBtd40bXtbrsK8e9n43+rWx0z0+fPp28nn1lt6u2wsNraGiQTp06dc/jZrNZamtrk7q6uqS2tjZJku5tajc0NCSbXx0dHcnnNnpPT09P8nGPxyM1NTVJHo9HMpvNyWPa7Xapp6dH6urqkhoaGpKPNzU1Jd8LJJu3ZrN53aZuT09P8rxtbW3J997v/Svl6OnpWdPcXWl+3u9e3H1PV5rhG12Hx+ORAKmnp+ee8q93n1b+f3XZ7Xa7NDg4uOG1rneO9a7xfp/X6rJv1NxO5f6sd8+7urqS//b2WzNe1CyfIF6vF7vdTkdHB/Pz8+vWsjay0Xva29uTzWyz2UxHR8c9711p0nV1da1pkjc2Nq6pva5u/q9X46mrq6Ojo4Pe3l4GBwcZHBxc8/zd71890rtRV8Bm7sVG17FynXV1dfc01R90n+4u+0bXevc5NrrGja7r7Nmza2p8D2puP+j+3F3uY8eO0dnZSXNzM8eOHbvvsZ80yt0ugLA1VppPDQ0NnD17lvr6elpaWlJ+v9vtXvc9brc7pabcitVNU6vV+lB9dl6vl7fffpt3332X+vp6enp67vt6s9mM3W6nt7eXurq6ZFfE3eV/2HuxUpYVqVzHZu5TKte60TVudF1OpxObzZZyOR72/pjNZnp6enj77beT5d4v/Z+iZvmYWv1lXvnitbe3YzabaWxspK2tLfmalX44WP5yrHecjd7T2NhIa2trMhCs1EDuPj/AiRMneO+995KPnz9//qGmJq0c22w231Or3Eh7eztnz57lzJkz69Z673cvNrLRddxv8Gyj+7Qer9e74bWud471rvF+n1d7e/ua4z2o3A9zf1ae7+jooK6u7p5pSU+03e4HEB5OV1eXZDabJbvdLp06dUo6deqU1NTUdM/0naamJgmQ6urqks+tTGNpaWlJ9qvZ7XbJ4/Fs+B5JkqSWlpbk4yt9WA0NDVJTU5P03e9+VwKS/Vjt7e1SS0uL1NHRkewL6+jokACpo6MjOc1p9dSnFT09PZLdbpdaWlqk9vb2ZF/o/d7f0NAgmc1mCUj2z/X09EhmszlZpo2ua/V5zWaz1NTUlLy+9a6jvb1dAjacorXefVp5T1dX15qyb3St651jvWu833U1NTVJdXV1UktLy7r926ncn43ueUdHR7J/d3Wf8n4gkySRKV14PDmdTrq7u5NNyJUa26lTp3a5ZFtnP1zj40I0w4XHVltb25rf3W73npzz+Sj2wzU+LkTNUnhsOZ1OTp48mRxdPnHixLqDPI+z/XCNjwsRLAVBEFIgmuGCIAgpEMFSEAQhBSJYCoIgpEAES0EQhBTsm+WONTU1OByO3S6GIAh70ODgINevX7/va/ZNsHQ4HHzve9/b7WIIgrAHvfHGGw98jWiGC4IgpEAEy3VIkoSYfioIwmoiWK7jpmeSn7n6CMaiu10UQRD2CBEs7zIb9HNx1sXkko8fjlxjamlht4skCMIesG8GeFIRikX5aGIQCYmr8+NoFSpCsSiHbPkctOUhk8l2u4iCIOySfRMsh4aGqK6u3vD53/qt3+JXf/M3SLDcV3lpdowb7kleLagkLiWYDi7wfK4DnVK9U0UWBGEP2TfBsrS0NKWpQ18vPsjHk4P83Yqn+RtXH91jNxlamOPvHDiONxzk+VwHuQbTDpRYEIS9RPRZ3kWrVPFqQQX1mUU0FFbxds2LzAYX+c5AL6F4lHOuW1yec4nRckHYZ/ZNzfJhyGQyamx5ZOrS0SlV/PbhV4km4gB4w0tcmhtjNugXzXJBSIHX68XtdmO1Wh/rzc1EsFxHPOBFplCSpU/na8UH+XTKyXjAS0KS+KMbn6BXqfmV8uVm+XO5DvJEs1wQ1tXZ2cn58+eB5R1IVwfN48ePP1bbY4hgeRcpkWDpxjkSIT86xzNo8yp5Jf8ANz2TXJx18UuOI/z57fP8h8s/463yY4TiUQ5a8ziUkY9cJno1BGG1pqamJyazuwiWdwmPXiIemEeKRQj0/ZzIrBN9xUtUW5eb5R9NqPntw+m819/Df7n5Ca8XLY+wz9xplhtUml2+AkHYW7xeL++88w42m435+XkcDsdD7+O+F4hguUp80U14pBeAhY//K1I4QPpzv0Z8YRqd41ky8yr5anENn045+R+qnuWD8X5yDUZgOVj+9fBVnsmxU5Ru3c3LEIQ9o7u7m9bWVs6dO5fsr6yvr8dut9PQ0LC7hXtIIliuItemocoqIzJ9G63jafyf/H+4v/u/kXb8TRKxSLKW+aWCCm64J1HK5STurCN/b6CHKksukUScMlMmx7KKUcoVu31JgrBrvF4vzc3NawIlgN1up6OjQwTLx5lMqUZf9QqqLDtytQFVloPF8534P/n/CI9cwtTwT/Cf70DreIbqvCqydOl8NDmAJ7REPJHgz25/wW1vMW+UHmI2uMgLuQ4sWsNuX5Yg7Iq3336bY8eO3bN1r9Pp3KUSPRoRLNehshWhfKqZ4MCnyNU6NMV1xH2TyGRyEtEwS30fEJ1xYq18ma8VH+SL6WF+5cBxymey+N7QFYYX5vmVA8fxR0MczSiiypqz25ckCDuuu7v7nn3PYXlU/MSJE7tQokcjguUGZEo1+sqXUWUu1zITkQAAS9d+QmTyNsbnf424fwat/WleyK8mz2BCJVdQYrTx329foN87TX6amZ7ZESaXfDybY0enVO3yVQnCznA6nXi9Xo4dO7bm8c7OTgAxwLNTWltbOX78OE6nc9vnaalshcu1zMHPiEzeQmkrJtj3Ie6/+t9IO/bLJGIRorNOSipeJqOklo8nB/jHtS8jv5N04+LsKP5ICHcoIOZkCvuO3W5f8/s777xDS0vLYzk5/bGbGNjZ2YnD4aCpqYnBwUF6e3u3/ZwypRp9xUsYDn0NXelxrN/8HTT24/g/+3O8P/09IrNO/Oc70Ez385XCqjtzLmWE41F+NHKd/+fSOS7NjfEz1y16ZkaJJxLbXmZB2E12ux273c6FCxeSj50+fRqr1Up7e/sulmzz9kzN0uv10t3dvWYCa3d3N11dXcmNxlpaWjh79iwnT54ElvfV6e7uvqcDebuorAXJWqZcrUNbUkdo+CIypQYpEWNp4GOUMwMcrniZ3IJKPply8s8OvUrHQA9/fPNTns91EEvEmVry8VyuA4tGvyPlFoTd0NXVRWtrK729vczPz2Oz2ejq6trtYm3angiWTqeT1tZWvF7vmmDZ2tpKT08PAI2NjRw7dgyv14vVujyP0W63b8vNn5r0k5amJi393gnmK7VMVaYduSYNdV4VAOGxqwQu/xDj879GYnEeS+FhvlZYy2czI/xG1bN8PDnIj0auo5IreL24hh+PXOdIRiGVlmyRJ1N4Iq3ULhsaGu6p0Hi93seuKb4ngqXdbufEiRNrqudnzpxZ09/R2NjIO++8g9lsxu12A8tBdqu3tw2FoowMu5EksFh0FBSaMaTdmyxDZS1AebyJkPMLwhM3kGsMSNEQ7u+/g772NRKxMMq5IV488BKDBjNKuQK7KROzRgfAZMBHNBFnPODh2Ry7WPkjPHE6Ozs5ffo0TqeTEydOJOdVnjlz5rFaE75iTwTL9fT09KwJlna7nbNnz/Ktb30rOU9rcHAw2SR/kEQiwcLC2i0iNBoNGs3aIDXu8iFJ4PEE8S+E8HiCWK168gtNGAxrg6ZMqUZ34AVUWWUEb3+A0lZI4MqPWbr6E8LDvZi+/I+JL3kpyqsmM/8gn8yO4Q0vEY7H+INrPydHb+Kt8nrcoQDHs0soNWZs5lYJwp7j9Xrp6uqira2N8+fP09raitvt5q233hJ9llvN7XavqaabzWacTidNTU20trYmB3pS7a/s7+/HZFo7Ev27v/u7fPvb307+HgxGmZtdniJ068YMkxMLFBaZqarOwu1ewmrTU1BgQn9X0FSac0g79ibhkYvIlRq0JXUErvwYhX65/CHXNZRzIzSUPcd1vYmbnkn+bsXTvNffw7+/dI5ftB8mHI/hWvTwVHYpGsWe/VgEISVms/mxDYob2bPfyrunHHi93uRj6010fZDy8vJkqqgVd9cq1WoFRcUWJsYXeOqZQoaHPPTdnME15qWwyEzt4Vzc80vYMgwUFJrQ6f523qRMrkBbegxVpp1g3wcoLfkAxHzT+Lp/n7SnmklEAlRm2snLP4hBqSHPYOJ7Q1c429/DzJKf14prmAku8myOXUwxEoQ9Zs8GS4fDsWbwxu123zPB9WHI5XKMRuN9X6NQyMnNM5KVncb0lB+NRklJqYUhp5uRIQ8KxfJAzNTkAvNzATIyDeQXrA2aijQrhrpvEnFdIzR0HplShcKYie/cf0JTUk/602+hd7toLD3GZV0uOqWaKksOWbrlsvnCS/zMdYsKczZHMwvF+nJB2CP2bLBsaWlZU4Nc6f/YrFQ2LPut3/otYDlo5uWbyMpOTwZNR5kNmUxGMBil68e3KSwyU1GVuSZoarXLQVMmk6EprEWVUcLS7Q+RN1gIO8/j/+I95v/q32B6+TdRxyMcNuWQn1eDRq4knIgRS8T5T1d/Tqkxg68W1zC1tMBzuXZs2rRNX7cgCFtjTwRLp9NJe3s7Fy5coLe3N9kP2d7eTmtrKw6Hg8bGxnua5g8j1Q3LVlMq5eQXmMjOSWdq0s/U5AIqlYLK6ixu35plZNhDSamFiqos5mYDZGSmkV9gTAZNuS6dtMNfIzLVj1ylRZ1fxWLPX6EwZgMQdbuwLMzQmFtFj8bEeGiRp7JL+OHINQZ8M7xVdoyFSIgaay61tnwU8sduDYEgPDFk0j7ZeeuNN9546GB5t1gszuSEn6lJP+FwlMGBefr75sjINPDMc8UAyOQyMu/UNDWav/1blIgECQ18SmRmYPn3aAjP9/5PNPbjGGpfR5FmZTKrgouxOBMBH2f7LzAZ8PFacTUv5x/ArNHzXI4dq8hiJAhbLpX4sG+CZW1tLfF4fMPnVzfDHyQajSdrmqFQjGg0gV6vYmLcx9xsgAOVmej0arKy0sjLN64JmtH5MYL9HxFfdBO4ujzNSGGwkv7sr6DOqyJiLqDXkMl0LMb747fJ0qdTa8tHkiTkMjkHbXnU2vLEFhaCsIVSCZZ7ohm+EzbTDN+ISqWgsMhMTm46kxMLTE/5icclwuE4w0Mehpxu7GU2DlRkMjOzSHZWGnkFJtRqxXJiDnMzoZFe5CoN2tKVNeb/AX3ta6TV/yLP+KdxWYpRFVaSACRJ4r/3XyBTl0YsEce16OG5HLvIlSnsC+sthV5tvWXR22HfBMvtoFItTzXKzTUyMbGAXC4jv8BE/+1ZBvvncQ64efHlUqSExMzMIlnZaeTlm1CrlejsT6HOLid4+0MUpmxCg5+jSFtexpkIeMiPhcnQmrlozGNWoSZDm8bPXH1cd0/yVlk9vkiQWms+NbZcUcsUnlgbLYVebb1l0duRL2LfBMuHGQ1/WCq1guISC7l5RibGfWi1SsrKMxgadGMyawFwjfkIhWLMziySlZ1OXp4RlcGC4cgvEJ26jVylRYqFkSQJ3/vvIsWjpD/7qzwTXsCVlo0y/wDV1lzeG+jh/73yNzQWVpGQJMYWPTybaxdJOYQn0npLoVfbaFl0R0fHlpdl3wTLrWyGb0StVlBSaiUvz8jE+AJanQopIRGNxuk5P4YkQfmBDMoOxJiZ9pOdk05unhF1bgVKWzGhwc+ITN/GcPQb+D/9czzffwdd9avkH/k6tsAcV9NzyT70CudcfcSk5TRv7nCAH41c47CtgCprjqhlCikLxiIMeGd35Fxl5kx0yntzLDyqjZZFb4d9Eyx3klqjpMRuJTd/OWjOzCzS+PoBbt+ape/WLAP9c5RXZBKLJZieWiQnN53cvHT0Va+gzq0gePtDVBmlLF3vJnD5h0TGr2N943/lKe8YY/4ZtJmFxHTpALzvug1IRBNxRvxuns0pFX2ZQkoGvLN89fv/cUfO9aNf+KfUZuRv+XE3Wha9HUSw3EYajZJSuzXZPNfpVJRXZHL71iyBxQgymYxoNM7IsJvpqZWaZvbyOvPRy8vrzEvriXkmkMnlJCJB8vwz2KIBbmjMzNhKCCdivO/q49r8BG+W1eGNLFFtEfMyhQcrM2fyo1/4pzt2ru1wv2XRW00Eyx2g1SqxO2zk5RsZd/nQ61UkEssztsZGvVy9PMmBykxCoRjTU35yco3kFBxBlV1G8PZHKNKX/6EtXe9m6epP0Nc0UH/wdVyBOTTmQqpqX6ZzoJf/eOVveCX/ALFEgrFFD8/klJJ5pwYqCHfTKdXbUtvbrDNnziRX6Z08eTKlNG5bvSz6fsQ8yzseZYDnYQWDUcZdPubnAgQCEW7dnGVkyI1araSiMpNShxWNRklunpGc3HQS7lGCA58QD7gJXP0pS1d/srw66KlmojlVDCDHlVVO99wkM0E/v3rgeDKhcKU5h8OZBajEGnPhMdXZ2Ul7e/uaoNjZ2ZkcHXc4HAwODgLQ3NxMW1vbQ9cuxaT0VbZiBc9WWwmac7N3guaNGUZHPNQdK6C4xAIsT0/KzTOSlakl5rpMeOwyUd80i190EHFdxfrGv0Jhzmc8uMCgPgN/VjkoVVx3T+D0zfFaUTUWjYGnc0pFJiPhseN0Ojl58iQXLlzg3LlzySlBFoslObizFfMsRbBcJdVgGb05TXzKj0whB4UM5HJkSjko5cg0SmQ61fJ/DWpkaerl1z2ipaUI42M+5ueXWFwMo9erkctlXL08iSFNTXGJBa1WRV6+EVt6gojzE2IeF9H5MVS2QiQpQfDWB0gl9ThjMaYyy/gwEuP7Q1dIU2n4ZcdRys1Z2I0Z1GcVi3yZgnAXsYJnExILYRKzAeJTfhILoTuBUoFMo0Bu0iLTrt37W6ZXIUvXIE/TIDfrkFt0yEzahwqier2a8opM8gMRXC4f7vklJEkiHI7Rf3uOvluzVFZlEQ5HmdCoyC94AXPOHHLnZyTCi8Tc4yye/w7yaz+l/HgTGTIZBpmKisqnOOvq5z/f+Jj6rCLeKD3ERMDHU9klFKVbt/rWCcITbd/ULFPtswx/Mkzc5SP4gxtEL7jWvEb9TBHa1yuJT/sJ/eAmcpv+zo8BuU2PIvvOYIpchtyoQW6981ymAXla6nvsBBYjuFxePO4gCwshbt2YwTXmIy1NzZe/Uo5CIUetVpCba8AU7Cc6fuVO0/w9Iq5rqPOqUR97kxGZgiljLj+Xa/l8bpyWgy8m+y4L0yw8lV2yLXPfBOFxI2qWq6Q8KT2+/LdD23gAzcsOiMWRogmkUAyZ/k6tUiZDbtYRn1kkemsGQjFkFh3pv/0iAKGf9CG36lHkGZFnpyOTy5BplcgzDMgz01DkpCNfZ+fIFYY0NRWVWSwuhhkf82E0aqmoCuGeX0KhkBOPJxgc8BIKxdDq8sgpLiDdexGlMZPw2BX8n7+Hevo2VZWvYAv6SAtO8HxeBWGZnOmlBX7m6uMbJbVMLy1wJLOQclOW2GFSEB5g3wTLVCkdNuQ56ZBILAfOeAIpnkCKxCEUQwrHkGmVKN6sBZaTXEhLUaSl6PLv0Tix27Mk5gIgARolikITul+qRQrFiI15kclkyAxqFNnpyHPTUWSmIVPfO1qdlqahoiqLRX+YsTEvJtPy0smZmUXOfz7GzevTVFZnES4yo9UeITvbgUGdhjq3EuRKkIH+1s+oSstkXCYxNz/MWHouA94Z/t3Fbr5aXEMoHsPpm+Pp7BIxmV0Q7kMEy7so8oykMslGisaRFsMk/GGkxQgJX4iENwj+MGm/9TxSJE58coH4iIf4uA/Zna0nlv5bLzKZDEWZDaUjA3mGHplchtxmQJFvRJFvuqfJnpauoao6G78/jGvUC8CXGsq4eX2aC1+46Ls5S+2RXMI56WjUL5ChcZHuv4EkxYj75whf+TGZGSXYjryBLuTHbsugMwp/6bxEz+wob5XVMx8KUGXN4ZAtX2xlIQjrEMFyk2QqBTKLHrllbQILKRIn4QuS8ARJzAVIVGQihWLJ55VlGcRuzxLuuk34x33ITFr0v3IUJEjMBYhenkRu0qLIN6HIN645fnq6hqqabBZ8IVxjPiwWHW73EjevzxCLLvfH+hdjhFR5aLRZZMQHMb70m0QrXmTx8/dIdP8H7I5nMdW+xm9GIryUlcefBwKo5AokJG64Jxn1u3kqq4S8NPOO3EdBeFyIYLnFZGoFisw0FJlpcGB55U3CHyYxFyA+u4hWq0R6thgpEiM27CE+MI/cogMg+JfXkMIxlJVZqA5kIruhQqZXLQfOAhPyDAMymQyjSUu1SYvPF8I15sVq1bMyTtd7wUVgMUJlTRbhvAo0yjxsmRbM33AQ7v+Y+OI8DnM2GaEAOu8UDrWOxaCHBYWSP775Ka8WVrAYDVOSbqM+q0gMAAnCHWI0/I6dXMGT8AaJT/mJT/tJzAbgztLHyOejRK9OEnf5QCZDUWJB+3pFcpRdplehKDSjLLYgN+uSx/N6g7jGfCz6w8zPBbhxfZrZmQAms5bqmmxyctJQhWfIiPZh0oSQyWDp1s8JXj9H6NDXGDfm4NGk88fBKFdDSxzJKODrJbXYtAYxACQ8drxe75rkGqkQo+Gr7ESKtlTJzTrkZh2qyiykWILEzCLxCR8yjQL100Uk/GFit2aI3ppBpl+u2UU+HwW5DGV1NrG+WWTpGpTFFhRFZsxmHWazDq8niGvMiy3DwOzsIjevzfD5p6O8/vUKZLocJtU2Zv2jZEpDGLLKCY9cRPXxn+DIq2au8hV+W53O5XQDf+Ke4N95pnmz7CiRRJwh3xxP5ZSKnJnCjktldY7X68VisSR/HxwcfOhgmYp9Eyz3KplSvjyolGdEqi9Ybq5PLKDISUd9vDD5uvjEAtErk/DDWygcVlQHc0nMBZBdm0Ju1aMssWAqMmM+lIvHE8SQpiYjw0AgEEGrVRGLJfjis3EcZVlEzdmooplkvFCAaa4X//lOLH/zB1ga/kdUMjmVGvhBHNKUy4NS00E/Pxq5RqUlh1pbvlhnLuyYVLKgv/fee3g8nm0JkKuJYLmHyGSyv+3vPJy33Fwf9xEb9aL7pYNovnKA2I1polcnCf3VNRQ5z6LISSc2OEd8dhHZ5QkUuUaMpVYO1ubg8QQZH/MRCEQIh2JEwjE+/nAYm01PVY2DsCIXtd5I1muH0Ex+gjavAktCYnzoAk2GDOITV1mQqvlPrgGK062EYlGGF+apzyqiON2227dLeESJ8BKRyVs7ci51biXyh2yZpJoFva2tjZ6eHpqbm2loaNiS8q5HBMs9LNlcr8kh4Q0SG/GgyDSgPl64vBQzXYMkSSz9+SWkQATVoVxUR/KIu3zItErSSqxUl5jxxSRcY15eftXB9NQiN65P89EHQ9gdVo4cPY7LPYY6PY2cUAh9dAr9R3+C3ZzLXMWXkAUXaAhH6HCPc3HOxS+U1hKIhskzzHI8uwSjWrvbt0nYpMjkLUa/fXxHzlX07fNoSx5uX5xUs6C3trbS1dVFY2MjbW1tKaV224x9M8CzF7MObYYkSctN9TEv8TEvUjhOfMpP9NI40SuTSEtR5Lnp6E8cSQ4Cya16FKUWfAYN45N+lpYiTE36UakUZGQa8HmDREMhjNIU2sgMGfFBuPxnxGadxHOrmCx/AY/awE9DEbpQU5Nj50T5MeQyGTXWPGqsuWJu5mNor9csm5ubsdvtyRyX3d3dNDc34/F41n19d3c3jY2Nm2qSiwGeJ9Dqprp0JJ/ElJ/YsBtFnhFNwwFiA3PEbk4ju7OcMvyhE3lOOsq5AHqNkspCEwvZ6ej0akLB5VVHTqeboUE3WdlpVJZVEcSIrs6O1X8B+eX/TulMP/oDL/NLyiDHgwGm40FIxAlKCXpnRxlamONYVjEFaZb7FV3YY+Qa/UPX9nbSw2ZBb2hooKGhgQsXLmxLc1wEy8eYTC7728GhcIzYqBdFlgFVZRYAUixB9Po0iXMDyExa1HUFqI7mYTBqqbDo8Ft1TCTgyNE8srLSuHF9mg8+niY7x0xteTZLcR36546QYwpRmaZi7moX1uACHo2W+M1z/GVMzhfxBG+UHmIxGqYgzcKxrGLSVKknDRGEFXdnSt9MFnS73b5t20rsm2b4Xppnud1W+jfjIx4SwSiJ8QUiPS6i1yZBJiP9X7ySXIsuKeUsmrRMqhQEFXJcYz76b8/x/IslaNRylqZHkPknSNdEMQ7/JfT9gITBxuyBl5k25fNBwMc5RRoFOWX8QmktNm0atbY8qiy5Yg8g4ZFtlAV9JVO61+sFSDa7m5ubN7UNrkj+u8qT0mf5MKSEtLw+3ekmPrmAFIoSH19A6bAhReIE/ugLVNXZqI7mI0tTs6hUMK1TsWTUglxGJBLnJz/sIyNDR3luCEN8Dn1smrSBs6imeohklDJ+5A1coSA/WQpyNS2L337qG6SpNBhVWo5ll4js7MIj2Wie5Uqm9N7eXlpbW2lqasLhcPDWW29tagqRCJarPC7BMp6IEYouEkuEicYjxBNR4tLy2nIZsuU+S5kKpUKNSqFBKdegVuoeuF94IhAhPuQm5pxHCsVILIYJd/cTvTYFcQllRSbq44Uo7Fb8kTizOhWLFj2j04vcujFDIBAhL9dAWaYfPT70AScWz8fI619jbNFDKLjARDyOzX6cxcwyRoMLFKfbKEqzUpdVJJrmwp4mBnj2uHgihi84zUJoDn94nkDYSyQeWvMaSZJISPFkoJStExRlyNAodWhV6WhVaWiVBgwaMwa1Ba1qOe2a3KBGfjAHZXX2cm1zcB55mgbtaxVEr0wSuTBG+P1BDA4b6WoFacEIi94Q6So5Jc8WMeQNcuvGDDdi6Tx/LJ/QrI6JNAemxTilwSkiH76LrfQZfBods8OX+EEggLXwEF8rqWU84KXGmku1GDUXHmMiWO6wRCLObGCUWf8I3uAUCSlBMOpnbnEMz9IE/tA81bkvYdHncm3ib7gy0Y0kJZLvL7Ee4XnHCZYiPt7v/1M0Sj0GtYU0jZU0jZUi68E1tUyVXI1BY8GgMWPUZmLSZqLON6HMNy3XNp3LiTxUTxUi3Rkdj494WPqzi6iP5FN0rIBAOEZaTMJxJAe/UYciTcN8QM3w7UnsFjdeeSWa4l8lffAvMLouoyt/kb+fls3n/R/wJ5N91JY9RSgeZdA3R31WkdjSQngsiWC5QyKxEOPem0z5nURioWRiiu5bf8i0f7kDW6M0YNRmEo2HAcgxOlArdSjlakBCkiTSNMuBRoaMDEMhodginqUJxjzXSEhxiq3LSYnfv/0nSEhY9HlY9LlY9XmkaazIZDK0Sj1GbSZGXSaWihy0NTnLK4UG5kjMBpDbDKifLiLa4yLyxShKh438Z4sJKmV4XT6CmWkkkDHnS+CaSqMg24Q968sEM4+hGTmH+epPOVz1HHkFh6h0T/HFxe8zIH2VqsIaPpjoJ0dvpD6rWKw1Fx4rIlhus3gixrj3JmPemyxF/DjnLjAwe54vV/wmerWJsszjlGc9TYahEL3alAyiKrmaQks1JbbDKOQqkCQkJEAiloiiVmh5uuSX7jy2LBoPJ99vNeQzHxjDOddDMLoAQEPF22Qb7Yx7+3BxC4s+D4VciVZpwKzPwXI8F2M4E+XQAnKzDs3LDqLXp4h8MYo07sNUlkGaUo7fHUAXilNUbGYwLnFlzMf4jJH6A0YyK36RqdznsGQYyNb4sERuU2UwoZzrJ5wI0Y0Ge1Yx00t+KszZ1Gbki90mhceC+Fe6jXzBGfqmPyEQ8dE3/SnXJ/+GWCJKkbWWhLQ8janUdgSjNgOjNoM0rQ292oRWaUAhT+2jicZDhGNLBKOLhO78LEUXqCv8KrHEcrM6FF1kPuDCZigA4PrkzxnzXkcpV5OVXkqusZxCSw0GjRmA9Gwr1sIcjHMGtOYS1IfzkOLLXQHRz0eRfzpC5uFcItU56BQyyvLScMYlsirykMf9zNweJroYQBfWoxhdoNh7jnBBGXOOZ1HMTvIjrZHssmcJFh9k2D/HkYxCHKZMkQZO2NNEsNwGkiQx6rnGqPsaEhJfjHyXobmLlGU9RW3eq+hU6Vh02WQbHVj0uagUmx8pVim0qBTaZPN8tVA0wGLYTSDixR+aYyE0RywR5XnH38GzNMm0f5BJ3wAXXT9CozJQqjnC7OIIw/NXyDMdQK3RoqnSYwlZMU1pMMyoUT9dhEypIHJ+FHnPOBllGUSeKkSrVhC7OE4410jfmIZAWElhRhR71d8nOP8s2tEu8ib/nL9TcpDP00xcu/pD3hu9zMHqVwjFY9z2znA8u5hMXfqj3HpB2DYiWG6xRCLOrelPmFkcJhD2kK7N4HD+V6jJeRmjLpPs9FKKLDXo1MZtL4tWZUCrMpDBcqo3SZJYivjwhWZYCM6SZyqnJjdELB6BO7W6ce8trk++j0wmJzvdToG5ikJLDfpCE4psGWZ3GmaTEePzLxC/Nk3k01GMejWmPCPeYTcLQ/N806hhOJbg0hy45i0U2uqpqLczMX4NbWCQA45CitKt3JifIHGzG436q7gt+fxk9AalxgzqMgtFhnZhzxHBcgvFEzFuTH3A7OIYHw38Gd7gFL9Q+z+jVxsxaTMpyzyebOruBplMtjylSGMmz3QAgEDYi2dpEk9wkoXgLEcKXqM882lc3pu4vDfoHfshkiRRmfM8PsmN2zSNKSMDmXsEi96A5VAZelk6cpkMdVc/toUQsYM5qIstlFi1jIbjTMa1aDLNyFVaQtF6lOE4BvdFvub8GH/FMwSHvuBqX4RLadk8f+AZXIseDlrzqLTkiFVA+0wqyX5Tec122DeT0rd7uWMiEefa5PvMLo7y8/4/ZS4wxktlv0a+qYJiay2Flpo93yeXSMTxhWZxB1zMB1yEYkt3Ru6Xm/u9oz/k5vSHmHU5FFsPUWQ5SFrIgHw6jCVgweo2oflsgdilCUhAoiITX20ucbWCRJoGlyzBB04vxRlhitOmiY9+jsnTg8EY5XbhQT6JRJjSGMmtfoW6olrMGh11mWKq0X5SX1+/JtlvW1vbPcl+U3nNwxKT0lfZ7m0lbs98xlzAxc/6/ghfaJpXD/wDck1l1OS8jFmfvW3n3UpyuQKLPgeLPgdH5jECYS/zgTHmA+P4w24OFzSSmV7MiPsq1yff5/L4T6kr/DpVlS8wvxhidnoOtRUyXs3HdBGUV33kl1rxBSIEht3kp6l5Bbg+IuNjfSEFmVnIcg7hHfucgus3eKOmml6VAmfPX/KjwfN86ZlmPoiGydalU59VjFXsa/5ESyXZb6oJgbfDvgmW22nce4uZxRF8S1Msht28euA3yTU6qMl7BaM2Y7eLt2krTfYiay3h2BJzi6NY9DkUWmqIxSNM+G5j1ucAcMv/Ba7ITUpzDxOIxVHoY6Q/YyRLNoklakH11zeRzDoMVVnk5JuY8gS55ZQxX15K4ZEsAjMHWdTaKFX5qFVfoD+6QPHAhwSzy+lLz2U66MduzOCo6M/cMpFInJnpxR05V1Z2Gmr1/VdvpZLsN9WEwNtBBMtH5A/NMzR/EQBbWiHfPPQv0aj01Oa/uu4I9eNKo9STb64k31xJKBpIBk5/2A2ASZfF1MIA58e/R49MQYG5mgrNMRZ9iyAlML2tw/x5BMPHQxgMGtKrsskpz0AelxFZNHA5rCMuycg16YgNRjmacCIPePD45+iZdaEtrmOx5kuM+t3UWPOosuaIpZOPaGZ6kf/wf3+4I+f6Z//8RQoK759Uxe12r0mCYTabcTqdD/2a7SKC5SNISAluz3zG7OIYfdOf8FTxL6JUqKnKfv6JCpR306oMFFiqKLBU3VmqOYpebSTfXMlSxMfQ/CWG5nqJZiZQFeQwN+3Eq/PizTaiegXMF0PYfIuobEX4vUF0MzIKPUoGFsJcMlnJLvy7xBY/QzVyk/yxz/l6rpn3h3v4wnUV24EXCJU/zYBvlqOZhZQYxV5Am5WVncY/++cv7ti5HiSVZL8PmxB4K4lg+QjGvTfxh9x8NvQdZDI5cpmCIksNVkP+bhdtx+hU6RRaaii01OAPzTPjH8Kky6I656XlF8hk3Fz8nNGlq1g0udizayh8qYRZ2TwW2S1Mk0p0fzmNwW4jtzyT8aUIozfgtuklDtdU4BnvReaP8my5HVdokOGrP2XKPYKutpGPYmH6PFPUZxWToXvwl1FYS61WPLC2t5NSSfa7mYTAW2VT8zK+853v8K1vfSv5u8/n27ICPS4isRCjnuvcmv4YX3CGZ0p+GaPORtGdtdn7UbrWhiPzGE+X/BIH814mM60IuUzG8/YTvFT299BrTfQuneP7sj9hzjSPRzbH8AEXfSe8zEn9pJ+7QvXgPM/Z9DyjkpM9byEj8wUWMl5DEc2h0JfFm7ML1C25MQ98zFzv97gxdp0fj17n44kBAtHwbt8C4RG0tLTQ29ub/L2rq4vW1lYAOjs7H/ia7bapmuWbb76J2WzmxIkT/KN/9I/40pe+xB/+4R/yD//hP9zq8u1Zo56rLEV8XJv8G8oyj2M15FOe+fQD80ruB3KZHJuhAJuhgGg8xOziKEZdBoWWapYiPgbnesjMrkSVUHFl4KcoDkgUO8qZH19Ed9mNWSkjI7MUWdhAdHYJRhaZMlrQ51SgZgLt5Ul0pnGWjCpGnb248ipYqv0KY4seqqy5YgO1x1h7ezutra04HA4aGxuTTey3336buro67Hb7hq/ZbpuaZ/mVr3wFmUyGJEnJzlWZTEZ/f/+WF3CrbGXy33BsifMj38U5d5HzI9/lG7X/nFLbYQ5kPbMlx39SLYbdTC0MMu0fIp5YTmj8ifM9hucvI5fJKZA5KE1UkUEuyoSa9K4gtgEtYXs+ThXM+qMsmNQYsqeweH+OKrLAdHEel+Re4iTIKDlK+cEG9HojRzMKKDVm7Pm5rcLesG3zLE+ePMmbb76Z/N3r9XLu3LnNHOqxNO69RUKSKLEdId9chUapo9h6aLeLteelaayUZVoptR1ldnGEqYVBnrO/xdHCr+Kc62Vw9jwj4b/i69rfQB+WMf90hBnrFKZbI5T5s3DkFTMYi6GZM6LP/SpT8VFkBhuH5TGkuc+Ydl0nQwFLWQ4+CQW4aZiiLrOIXLG1hbAFNhUsnU4nfr+f9PTlpAdmsxmLZX9sgxqLR5hcGMC7NEW6NgOVQkOO0YFGKXIzpkohV5JjdJBjdBAIe5laGMCosVGd8xLupQnMhnzi3gBdgx1kHM/BfrAC3+1Z1IMujKVPka7KR+b3oxpXEXFpUBQGsc1GyZkNIw/1seiZ4erFH1NU9TKe0BHy06wczSzELPJnCo9gU8GypaWFN998E4vFwvHjx4HljtZXX311Swt3P06nc8f6KlabXRwhGg/zs9t/RIn1CPVFXyffXLnj5XhSGDRmHJnHKLUdZS4wxoTvNguhOaR0FfnZNQzMnGdAd5WswwWUHawhiynmFyfI/txPsbwIlVbOzHUV8+lvITcOY3V9imVykVKTxJXA97ENfs5SzZeZWPTiMGdyOKNATGoXNmVTwdJkMvHTn/6U73znO8lh/D/4gz/Y0oLdT29vL1/+8pfxeDw7ds4V034nk75+glE/JbbD2Az56FQirdijkssVZKWXkJVewmLYzYSvH41SR23elxn1XOX21KdcCX/O65IDVbqM6aNB5EM9GN06apV2ZpaszPiLMOYbWUpcJSMhUaEtZNTTz8WP/ox6x1FGyp9jeGFeTGoXNuWR5lm++eaba/oud0pdXR1W685P+l6K+FgIzTM4dwGLLheLPo9s487Xbp90aRorB7KeptR2hOkFJ2kaM6W2o4RjS6gTarzjw/xY/gMKyuyUjhezcOsG6riO3LIjWBeNmEP13JrRYlREqFf40PoW0Bg8qPmYMYWGS47j3PblcdiWL5IOCynb8UnpXq+X7u5umpqako/tVsqlhzW7OEo4FmTce4ujBV9Fo9Rh0eftdrGeWCqFhgJLFfnmSjxLk0z4buNemsBQmM8RdSN9058yWthPRkEujmANptg4XkWUPCdkJfKYWkhD0tYS0alR37hE0GVkRB8m4OylqOQwocoXuWXOpS6zSOxvLjzQjgZLp9NJa2srXq93TbBsbW1dk3Lp2LFjj5xyaTvMB8YIRhew6PMosh4kM71YzKvcATKZDKshD6shj1B0kUlfPwa1iYq8Fxl3X+fm+IfMKaYpkg6wFA4xkTmJ1uuiKJZBKJzPTKiKkLYEU/RDji/46cst4fbtfibHblBadowFx9NkW/KoyyoSm6gJG9rRYGm32zlx4gTt7e3JxzZKufTuu+/yzjvvrHm/w+HYtVpnKLrIYtiLWZfN69X/BICMO3vaCDtHq0qjNOMoRdZaZvxDpGnMFNpqScSj4A4zrrzGtcxPKQzYKb4RI316Fos5G7kmG2vsWXwLXrKXrOgyF4jNXccTu0FpeBF/Zik/9B3AnlHAYVsBepUYBHpceL3eNck1tsuurw3fKOWS2Wymra1ty86TSCRYWFhY85hGo0GjSW3/m/nAOPFEDO/SJFZDPmqljnRt5paVT3g4CrmSXFM5uaZyPEuTjHtv4VZMYjc+R3Q0wQA9DB3vIydcQMmSAV1cQ3QsgG7ehHxIg3UiC61iGoW7n8CiGnf2ArM3P4SK5xktOEhlZhHV1lxUYhBoR6XSJef1etdMVRwcHNwfwXIzKZd6e3txOp0PNX2ov78fk2ltv9Tv/u7v8u1vfzul93uWJpnxD/Gz23/E12t+m4rsZ0UTfI+w6HOx6HNZiviY8N0mTWvlYPwrDI+d59bcJyg0CYxRHQuTU0iKWaz6UeLxPLyhKlSRKtKWesiaGSVg1PGZ92eUOi9AxfP051VSm1lCuTlTfNY7JJUuuffeew+Px7MjAXK1XQ+Wm0m5VFdXx8Ou0iwvL+f8+fNrHku1VilJEguhWab9Q2iUBky6bCz63Ic6v7D9VvZhL7YeYmphAIPaiKPoWSR/iMT0Ip+rbuGpmaV4vJjSG0EsChfutAI0shI03jRKZTICRgMjrhlccz/maftNrhYd4VZOBYcziyhOt4qR822Uahb0trY2enp6aG5upqGhYcfKt+vBcqdSLsnlcozGze2oGIh4iSWizPiHyEovRSaTYdJlbXEJha2iUmgotNSQb65ifnGMcd0tFozz1Hve4ObYzxkovM5ggZzC+WLyJStzYT/mYRk+TwmOW+mUp8tQME1cWkIVvkFg7ApfFNdxI6eCuuxicvTbvzPnTpCCUWJD7h05l7LUikynuu9rUs2C3traSldXV3L/nVOnTm15edez68GypaVlTd9kV1fXlvZVrhgaGqK6unrD5++3YZkvOEM8EWU+MMbRgq+iU6WJ5Y2PAblMTmZ6MZnpxSwEZ3EZbmG25HN4wU3f6IcMya5RFTUQcydYCjmxKDwElUYCS8XoQ1/GsHCTmOs2TqOCpWEXB4oL+cDxFLbcCuoyC7E85nsCxYbczDf96Y6cy9b566iq778XVapdci0tLbS0tNDd3U1jYyMtLS1PXp+l0+mkvb2dCxcu0Nvbm+yL2ImUS4+yYdlCaJaliA+D2kxWeqmoVT6GjLpMqnWZBCMLjPv6SDNlcCj4OkwHCKk9/CDrPKZFI45b+WTOegmY05mSFVA+L/GC2835hJnzSx4yxr/LgdJSflpSR0FeJYczCkhTpdads9coS63YOn99x871IA/bJdfQ0EBDQwMXLlzYkeb4jk8dWt3kXrFy0XuVPzRHujaDNw79CwCMYhT8saVTG+/0a9Yy4etnwnAbsg08NfoaN2Wf8cXxy6QFjRTP2DFLMsZCMXTTeeRPlZFlChByX2dycYIDgRju6X5+mFeNI6+CGmseWuX9m5l7jUynemBtbydtpkvObreLbSW22mab4dF4mFBsiVg8glKxPPcu/THesVFYplJoKbbWUmCuYtrvJC3dRkngKSbHbnDT/QlThdPk+o4i88wiC02CzkskmI8pUEXxzDCh+XEm8vwEBm4hd1QyVHSEqrwDVFiyxZrzTbpfl1xnZydNTU14vV6AZLPb7XaLYLnVNtsMX4p4Afjxjd8n31xJfdHXROKMJ4hCriTPdIBcYznzgTFMxmzyFmuITHlBFWRMF+f8oc8oHs3DPrBILD2NUXLImo1SMzlBf4aCDz23cQz3Iy+voa/wCIdyy7CbMsR0o014UKb03t5eWltbaWpqwuFw8O677+5Y2fZNsNysxbCXeCLKQmiOCo0NvdokvgRPIJlMRkZaERlpRfiCM4ybbjGXM0L2lJyyqWkGS64yVOyicDwH20KYOZWGJY8CizefI650hrxzTMxc5kipi6tFDm4VHuJwjoPC9Cd3l8/tsFGX3EqGMbvdvmap9E7aN8Fys83wQMSLLziDRAKLPgeDen8kOd7PTLosTLoslmw+xq19GPPyqJl6mf7Jz+gruIRqaZFCbxbasXkiikli8gjlHgtZY0ss+QIE/ZMERga4UFrBjcIjHM11kKUXrZHH3b4Jlo/SDPcEpwAZJl02BrV5y8sm7E16tYnyrKcottYybusjLTeXiukXiE/5iSniXDl6hYA0R3lfEeqolRGbGf28njzXELczZPS5/JQX3uIjew3WkjqOZJeKbO2PsX0TLDdDkiQCYS++pWnSNVZUCg0GjahZ7jdqpY5S2xEKzdVMZvTjyr1JeMZNnsvB9cgcnz53CbPPiGOgiDlTAfG4n8IZGfqBDPrdCRSTF6kYuY27pIoc+zEOZxVjeEynG+1nIljeRzi2RFyKc7TwdapyXwTAoBZ5D/crpUJNoaWGPFMF0xlO9DlZFE/X4XJd4UbaF/TUX6N+OguZ10TQGyWoCVC4ZMV2xcPChA/8I8xNDPGTwnKKSo9Tk1mA7jGbbrSf7ZtguZk+y3AsAIBMJkenSkcuk6NW6ra1nMLetzKCnmMsYy5jlPS8fAqnDjPvcqI2GViK+Lh05GdkT5nJmixg1mpAu6TF+sEoMzl6XCNfkBjtw1Vag91+nEprLmrFvvkqPrb2zSe0mT7LUCyAJEmc6/tDqnJepDzrqW0qnfA4ksvkyX2D3LZxRvOu450cRelSol40cd02SN+BEexD+ahmSum3QL5nmtqLJnonFejGPycx0sdQaQ0HHE9xwJyFQi5mWuxV+yZYbkY4GiAcCzDtd3Ig61mxHlzYkNWQj9WQj882w1j+dXQTmfhGR7ke/4L+8ltoSmY5OPocC+iJG8AaTaC4LsPlnMQ4G+X2WB8DJbVUO+opNdrE9LQ9SATL+wjFAgTCy/O70jQWtMq0XS6RsNetTDsK2I4yVnAds6sIv2uSOc8kWpOR2MIYE7m3KBwpJqA3Y4jKyfxsknBeNpNjnxAbuU6f4wgHS49SJOZo7ikiWN5HOBYgcGcFj0FjRqt6vLPMCDvHoDFTmfM8JbbDuIpvMjF6nci4mxGjm3HVHCNFkxSMZ1M0XMyALYMs7zVyXOkMDwbQjHoID13hVlkdtSVHyN1Hm6mtt6Hharu5ueG+CZabGeAJRZeDpVKuRq3Qo1GKYCk8HK0qjbLM4xRZahkvuYV2NIP8UTt93l4Gcy8xVjBF3mQlereVqNqNCh3R6zA4PItp9ByBkiuYDxznUPFBbNonu2Wz0YaGq+3m5ob7Jlg+7ACPJEmEYwFyTQfQqdKRyWRoRbAUNkmt1FJqO0KBuYqJ4tukj+ZTOXKM256L6FTpaDQyYvIAOmmKmCoTTdSPrQfiExpmJn7GB4W9ZFQ8w6GCKkyarZmREYkFmVoY2JJjPUiOseyBM0nW29BwtVQzqW+XfRMsH1Y0HiIhJTDrsjHrltNYaUQzXHhEKoXmTrajSiaK+0kfzScwOk1g2s/5vEvMGSbJmLVQ6ixhTp+FfLEf6+dKJm5m4B+aY87xBTlVz3Mo78AjT2yfWhjg//jJN7boyu7vX732A4qstY90jFQzqW8XESw3ELozx3Jg9gssulwy04tQK8QcS2FrKOQqCi3V5JkOMFU8yOjoRV4eNjA4d40BYw/nn76IyWsid7wahdJHWmQO79ViFobiLDn/iumyHAqrX6Qmx7HpPJo5xjL+1Ws/2OIr2/hcj2ozmxtuJREsNxCNhwC47OriQPaz5JsrxGZVwpZTyJXkmyvINZUzU+LENFKIY6iW4bkbjKhvYNKmoQxqiMv96KQgi/EgiZ44wTE5zuG/wuXIobT6BSqzSh56YrtaqXvk2t5O2szmhltJBMsNxOKR5X7L+BIapS6Z+FcQtoNcJifHWEbWQTtzJaOYR4ooHa5hMe7Dp5jg0+L3SVvUY3eWENbksRToI+OigfmBEKHBCUYOFFFW/RIHMgqe2IntO7W54Ub2TbB82NHwWCJCLBFBkhKoFXqUcpH4QNh+K6uCMmuKmS8ZY2SsF/egiaenFPTLL3Dl8A30ASeZ42XIfQuoQ6P4L+ayMBgleNvFYFUpB6pfpMya88RMbF/Jkr5TmxtuZN8Ey4cdDY/GI0RiQWC5uaJSiGAp7JxkMuKqItxFE2SM9VAwUMHE5AD9XCCQ5SYjXo56JkxUGccbHyV6WYXWGSbYP8JApYPK6hcoNWU9Nt1HG21ouJIl3W6378jmhhvZN8HyYcUSYQDyzVUY1CYUcpEdRtgdVkMe1so8vEXTjLh6yesvY9E1T0DvYyrvFkOFPZQOFaMKFhEJ30a6aMQ/lODCrUFuVx+guuoFikx7f9+ojTY0XMmSDru7uaEIlhuIxSMYNGZeKV/eKlTULIWNOJ3OHanhmPXZmA98lYWi44yM9TJ5+wYylxxPeIm+ikFU0RFKhgsJLKqxLn5I4kI2gYEQ/lv93KqpoqbyefLTRT7WzXoyOjW2QTQRJhaPEIouZx5SyndvgKezsxOZTEZ9fT2nT5+mtbWVkydPUl9fn3yNw+FI7ny31/T29nLy5El6e3t3uyj36OzspLGxEYfDwZkzZ+55/u4yd3Z20trayunTpzl58iQymYy2tjZ6e3t37P4btRnUln+F5xr+B442vMELhU284P9FrMFc+suHCGVNEVco0CX8zEWHmL44zORf9fLJX/xnzvX8iOnAwo6U80kjaparSJJEOLzc/A4EFumfvMAXI9/nRP3vEo9KhEKhbTmvRqO5b79SU1MTdXV1nDhxglOnTiUfb21tTf5/R0fHmjlo2211n9KD1NXV0d3dTXNz866W4+73vfPOO8mmX29vL/X19WuWz50+fXrN/T558iRut3vNipHGxkbcbjd1dXWcOXNmR9cqp2ks1DhepaSwnpGxXrJulTI7OkQ0HkGmDOIsvYEmrEU1W0pg6Saqi2bmnVE+uH6DzEOHqa14BpvuyV5CuZVEsFwlHA7zp3/6pwAMz19mxH2VCV+Yn459wjXjDGna7ckC8+u//utotdqHft/qkcCdWh8Ly83O5uZmBgcHU37PdgTyzZQDSAbuc+fOJe9bXV0dTU1NnD17lrq6Ok6fPr1mfXJnZyfvvfceQ0NDa461ei/rt956i9bW1h0doYXl7P3Vji9RUljHqOsSozcusTAyRRwnzvxR5CVjFI8UIJtTo1i8iOqCkdDtJWYrrpBzpI7aA0+JvYFSIJrhG4gnYsQSEZR3Bnbk8r31d2V1k3GlVtTb20tnZyf19fV0dnZisVjwer2cOXOGM2fOJJuOK1pbW5OBY6OVECvva2xsTDZJe3t7cTqdnDlz5p6mZyrH3Mjp06exWJb71Lq7u5HJZDidTjo7O5PN5Pr6eurr6/F6vfeUY6PrXM3r9dLc3My777677h+Yletpb29f0w/59ttv861vfWvdoL/ymNlsTpZpN+jVJirtL/Pia/+Ag6+/zrPFv8wLvm+S7S9kqHSMa/Ufk5BHUMjG8QUvMXOxn9GOD+nuPMOnNz/CH9meltOTYt8Ey5V5lhv9/P7v/37ytZIkEZfiy8HyzmR0uUyxW0VPOnv2bLK/cnUwqKurS37JGxoa6O3txW6309PTw4ULF+jo6KClpYVTp07hdrs5ffp08kvd0NDA8ePH6ezsvOd83d3deL1eTp06lUxasHIOs9lMS0vLmuCRyjHvZ3UTtqGhIRmsGhoaksft6enBarXS2tq6phwbXefd3nnnHaxW67pZbVb/0VkdKJ1OJ16vd80orNfrTQbl1ddpt9t3vW9Wp0qn0vEyL732mxz76ps8U/zLvOj9RRzztcjTLKjDFpbMs8hjAaaX+vCev8bY2Z/zk852vrj9OUvRyK6Wf6/aW9WlbZTKPMuVPsmEFAMgFo8mB3YUe6BmubrPcvXgDqyt3cDfNsvb29vX1KAaGxtpb2/n1KlTdHR00Nvbu2EzdqUfdL2Bj/XU1dU98JibsXJNKwGsubn5nsw0XV1dG17nat3d3etOPenu7sbpdPLWW2/R3d295rmVmuLqALpSpu7u7jVl2en1yvejVaVRYX+R4sI6xlyXcF47z8LwJNPqIYYLXEilYxSO5mGYLsXjv4ah18LA7QVGK85TevwZakoOb3rd+ZNo9yPAHqLRaPj1X/91glE/F0YNfEV6jngihlKh5rnSN7dtrqVG8/DTklZqYV6v94H9gaubylarFavVitfr5e233+bdd9+lvr4+mSNwNbfbTX19fcqDFqkccyvY7Xas1nv7j9e7zvWsd79aW1s5deoUZrP5nuc3mhZ0/vz5e2qoXq83mZh2r9CqDJSXPk9R4VFGxy7ivHYBy3A+Q9JlxgtGGCsaJ3O8AvuUh5B/FFlPAdE+D8NVX+B4+nmqC2vEhmrso2Z4KmQyGVqtFqVKjkajRqvVYtCnodNqMOjT0Wq12/KT6gqL+fn5Nb93dnYmazErgeLuPsQTJ07w3nvvJX8/f/48zc3Nydqi2WzesBbY2NhIW1tb8pirm5urz7PyeCrHvPs4q5nNZrxeb/Jnpfl7t66uruTI+srzG13n3U6cOHHP+Zubm7Hb7cmBmWPHjuF2u5PP2+12Ghoakt0Qq6+jsbFxzWNut3tHV5U8DI1ST3np83zpqy08/Y23OF7yTZ73vkG+u4w0mZao1kRaUE0s5mcs2M/8F5e59cff46//4g+5NnqNWCK+25ewq8Sfi3XEEst9Nh/0/zcKLTVUZD+7q+Xp7Oy8Zx6f0+mku7sbj8eT7Cs8e/ZsslazMo2lrq6OtrY2Tp48SWNjIzabjZaWFnp7e2lvb0/O11yZK7i6KdvS0kJXVxcWi4W6ujreffddYDmoNTQ00NzcTFtbW3I5WkNDw7rHXClvV1cXDQ0Na5av3e3UqVPU19fT1NREQ0MD3d3da6bymM3m5DUAa8qx3nWud3z42wDpdDppbGxc81qz2XxPrbSrqyvZV7wS1Nva2u5p0t/dt7kXqZU6ykqfpajwKC7XZWxXSvA6XSxpvfTlXGY2Y4b88Rzik6V4Fm+QdsHC9ZszDBzMpfLplyjPLX9ik3Xcj0ySJGm3C7ET3njjjZTXhs/6R7g5/TFne36XQ/mN1Bd9jWNFO5MkVVifTCZjJ/+prvwxeph5kysDYhttibBXReMhxsau0H/5U+aGhhjhOpPmYaKqCNlTuVhmC9DGoyjIJs2Yj/5gITXPvYI9q+SJSdaRSnwQNct1RBNhEonl0XC1Urerq3eEe7sWdoLdbk+ObKcyh3WlO+RxC5QAKoUWe8lTFBYexuW6QsbFT5h1DjHGDSatY2R68jEEEkQkHy5FCOPnkwSvjXHzUBEHn3uVkoyixyZZx6MQwXId8USM6J1EGiq5JjnXUtgdK32RO71CpqGhIeVAbbVa92xfZapUCg2lxccpLFgOmn0XP2bOOUxQ4yeuDnLx0GfY5q1ox+3M+9yEP5vmi6tj3DhSyqHnXqXQmr/bl7CtRLB8gOW/mE/+X829rKWlZUeD5GqprjzayaWm202pUFNSfIzCwsOMjV2hr/cjZpxOCubKmbQMMZX7GTmTmUSmweCfxf/RBIuXhzAeLePQs6+Sb8nZ7UvYFiJYbkApV/F0yS9j0eftdlEEYVco5CpKiuvvNM+vcuvCB+Q5B5ngNpPWIRbS+jh0Mw9lNMSE5MP/4RQLFwcwH6vg8DOvkmPK3O1L2FL7Jlg+bKZ0hVxFWebxnSiaIOxpCrmS4qKjFBTU4nJdo+/Cz5kc7Gcx4iaarkJa8ODJ6Ucx4WDJ6yHw/iS+nj4sx6s58syrZKZvT06FnbZvguXDZkqPJ6IMzV8ix1iGVdQuBeFO0DxCQcFBxl3XuHn+b5hxDjGj9+Iz+ZnO+YzsqUwU4yVEFnwsvT+N5/wNbM/UcuTpL2EzmHf7Eh7JvgmWDyuWiPD58F/wUtmvUWSp2e3iCMKeoZArKSo6Qv6doHnji5+R4SxgUjbApG2Y6ZzzFAxWke/2MBudxN89zfzn18h87giHj7+MVW/c7UvYFBEsBUHYlLuD5vXPz5HtHGRS5cQYs6KIh5GpAgSisywsLLD440lmP71E9nNHOfLUlzBpDbt9CQ/lyZhRuo90d3dvmH5st7Kl369MwpNvJWh+5c3f5su/0sKRggbStRai6ekMOsa4Xvc5nqwbBKPTeLy3GPrpj/nJ7/17Pvzkx/jDwd0ufspEsHxMrCwbtNvtG2YB2uls6SuOHTuWcmYi4cm1EjRfb/qfeOVX3ib3YA217hcpdVXhMwW4Un+BofzbhEIuZj2X6P/r7/HD3/t3fPz5OQLhvZ9LUzTDN+Ce8ROfzmBc8qJy3yYyZcRisVBaWkooFOLGjRv3vGdlpUdfXx+BQGDNcyUlJVitVmZnZxkbGwMgNzeX3NzcB5ZldUbwjTLprD7/TnuS5hgKj04uV1BcfITCwlpcrqtc/fgnZA8PMa0cxqfzoZQlMHvizKVP4/MuEvj+JOMffErhy89x+OgL6FV7c8WcqFlu4MedH/Of/6deWv/uf+LE68uJIX7nd34HAJfLlczYvfpnxW/8xm/c89wPf/hDYHk1yspjd+dk3Mh6mclXMt6sZL1ZnbgWlpNOrJcVZ8WZM2eSmchXspKv5HFsbW1N1hTXyz6+Xjb21WV1OBz35IQU9h+5XEFR0RG+euJf8OVf/UcczP8SFf4jJNJsuG1RrtZfwJ11nUB4irn56/R97y/4we/9X3zR+wHBPZiAWNQsN/ALJ77EM186jEwmx6zLoizzWHLLg4KCgvvmavzjP/7jdWuWsLxPy7PPLmcxSqVWCWszk68Eprq6Orq6unA4HDidzjXZ0ru7u5mfn0/mZ1xPS0tLcsOzhoaGNe9fydjT3d1NR0dHci/n5uZmTp8+TUtLSzJrT09PT/IcK5mPVj8mCHK5gqLiIxQU1uIavcqlj3/E3IgK+1iMSdswM3UXyJrOwDRXTJpnjsBfTDD8/kfYX32JQ7XP7JlcmnujFHtQuk3DT0b/My+V/RrV+Qc4mPe3TVytVnvfJm9FRcWGz2VmZpKZ+egrG1bWIa+kGbPb7ckAdezYMU6ePInT6UymVVtPS0sL7e3tyTyOZ8+eXZMI4kHZx+++B83NzbvWbyrsfXK5gqKSIxQU1TI2cpXLH/01GaNFzKlGmbANI4+6yRjXQDTCpORn8TvjDP3sIxyNL1Nb8xQq+e5u7SKC5RPIbDbT09PD22+/ncxYvl4AW8k76XA46OjooL6+njNnzqwJmKlmH185XnNz87ZlSBeeDHK5guLSIxQW1zI2fJVLH/6ADFchkViEhCrKTOYIi4YRpMlSlmbdLL43xqDt5xxo/BI1VcdQ7lLQFH2Wj4n7TQm6O0v6Sibwjo4O6urquHDhwrrvW0lD1tPTg91u56233qKjoyNZa90o+/hGZWlpacFut4tpREJK5HIFxfYj/MLf+xZf/jv/hLySMjCnIVeYWExf4uaR80zl9OGPTuCeu0bvf/+v/OV//L+4dvMi8URix8srapaPgdWZyQ8fPgws90va7XYuXLiQDHorfYbHjx+nubmZkydPcvz48ftm7v7Wt76VDI4rTfcVG2VZXz3409LSkgzOnZ2dtLW14XA4MJvNG24dKwiryeUKShxHKSo9xOjQFS598H2s4/nMq1xM2Ya5friXg1cOIg/7mI7Ps/hnY9zKKqL69a9SWXZwxxIQP5aZ0pubm+nu7uZb3/rWPbv3beRhMqWPeW5wa/pjOi/+77xU9msczm/kYN4rj1BiQRBSlUjEGR28xMUPf8D8xBge1RxGnwbzbIJ+xzgGbx4KpZ50eTrG7FJqXnudCkfNIyUgfiIzpa+M0Hq9XiwWS8rB8mGpFDp+6fD/glqp35bjC4KwPrlcQUl5PUWOI4z0X+LiR9/HHR3HnetnPsPDWPE4GTM2mCgiPOnn0z8d4nq2ncNf/xqO4spty9q+48HS6/XS3d29ZhChu7s7OQ0GuG+i15Um5UrTdLvIZXL0atO2HV8QhPuTyxWUVtRTXH6EkYFL9P78uygnDXhUE0xmDHP9yEUs0wU4xjOYnvTy4R85uZpbypFvfJPSwvItL8+OBkun00lra+s9mzq1trYmR1AbGxs5duzYA1ejeL3eDSdcb4VoPMxnQ51U5bwoUrQJwi6SyxWUHqinuOwII7cv0vPBdzFP5uFWT5CQEihiEXSBJWZMIQITi3jfHeJSfjl13/gmxfmlW1aOHQ2WdrudEydOrFm5cubMmTV7lzQ2NvLOO+/w7rvv3rNPs8PhSNY633vvvW1rggMkpBijnmuU2I5s2zkEQUidXK6gtPIYxRV1jNzq4cIH38U7PUU8LcpYppPh0lEyZmxkTRYRHl/kZ2cGKDr4NF9u/ntbc/4tOcojWJm2smJlVNdsNif3gl75WQmUq7coTTXLTiKR4N/+23/LwsJC8iccXt6U7Pd///e39qK2wXaV8VGO+zDvTfW1D3rd/Z7f6Ln9+vnu1Geb6uu36rOVy+SUVh3nzZZ/Q1CXhyk/H1usHPtINUuGIDdrLxGXR4nIA+iM2tQv4gF2PVi63e41v5vN5jXTV+7W2dlJc3MzDocjudQvFf39/fzLf/kvMZlMyZ+Vmut+/TI96nFFsNw6Ilg+/Gcrl8n5g//SSdPJ/50v/3ILJcajVE48S8VQHUqZBo02jeOvfvPBhU/Rro+G3719qNfrve+Wok1NTZvam7m8fLnD9/z588nHNBrNQx9HEIS9RS6T4zj4DKXVxxm89hm9H/4A3/wMVXWvoNrCDEa7XrO8u3bodrs5duzYlp9nZGSE4eFhnnnmmeTP0aNHqa6uZnBwcM1frmgkwk9/1EVtbgMmbdaWl0XYXeFwmG9/+9vJbhjhySCXKyg/9Dxv/ObvMBtJ5+Bzr23t8bf0aJvQ0tKSTCsGy8kbVrLhbKXS0lIcDgc3bty458fhcKzZ2TESjfKDH/yQAxnPYdQ9Wdt5CsvB8l//638tguUTKhqJ8b/+H79HPL61SyJ3NFg6nU7a29u5cOHCmgDZ3t6ezKHY2Nh432b4TpHJYWKhj2BkYbeLIgjCHrDjU4dWciOu1tDQsK0TzGF53/DBwcF19w5faYav2TdcJePTkbO8pPo18s2V21o2QRD2vl0f4NkppaWlxOPxdbeDqK6uXhMoBUEQ7vZYJtLYjJqaGuRyOaWl987oHxoaWvN4IpGgv7+f8vJy5PJd79YF7i3jXjjuw7w31dc+6HX3e36j5/br57tTn22qr9+JzxY29/kODg5y/fr1+75m3wRLQRCER7E3/qwKgiDscSJYCoIgpEAES0EQhBSIYCkIgpACESwF4Y5Uk7II+5MIlilqbW2ls7OT06dP73ZRhG3Q29tLfX39bhdD2AbNzc1YLJZH/u6KYJmCzs5OHA4HTU1NDA4OrlmqKTwZ6urq7rsnuvB4Wtmza2ho6JFzTuybFTzrSXU/oLNnzyb3wnY4HHR3dz9w2wth9z3qfk/C3pXqZ7uVe3bt22D5MPsBeb3eZK1jo/Xtwt6ylfs9CXvLZj7brdiza982w1f2A1pto/2AzGZzMqO70+lM/uUS9q6H+XyFx8tmPtut2LNr3wbL9Wy0H9CJEyeSI6WDg4PbniFJ2B4bfb7C4+9+n+1m9uxajwiWq2y0H9DKwM7KQI9otj2e7rffU29vL06nUwTPx9RGn+1m9+xaz77ts1zP/fYDamtr240iCVvofp9vXV0dIqfM42ujz3aze3atR9QsV9mp/YCE3SE+3yfXTny2IliuslP7AQm7Q3y+T66d+Gz3bT5Lp9PJyZMnuXDhAufOnUv2Q4p5eE8G8fk+uXbrs923wVIQBOFhiGa4IAhCCkSwFARBSIEIloIgCCkQwVIQBCEFIlgKgiCkQARLQRCEFIhgKQiCkAIRLAVBEFIggqUgCEIKRLAU9p3Ozs7k5mQnT57kzJkzu1wi4XEgljsK+5LFYqGtrU2sDRdSJmqWwr7U0NAgdnMUHooIlsK+5PV6OX/+/G4XQ3iMiGAp7DudnZ20tbUlt5IQhFSIYCnsG729vdTX11NXV0ddXR1Op3NNwlhBuB8xwCMIgpACUbMUBEFIgQiWgiAIKRDBUhAEIQUiWAqCIKRABEtBEIQUiGApCIKQAhEsBUEQUiCCpSAIQgpEsBQEQUiBCJaCIAgpEMFSEAQhBf8/OdGwt+8YLSQAAAAASUVORK5CYII=",
"text/plain": [
"<Figure size 340x200 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"q_list = [-1.0, -0.5, 0.0, 0.5, 1.0]\n",
"x_vals = np.logspace(0, 2, 200)\n",
"plt.figure(figsize=(figWidthsOneColDict[style], figHeightsDict[style]))\n",
"col_list = dark2\n",
"\n",
"def deflection_angle_x_inverse_correction(q, x):\n",
" return (1 / x) - q / x**2\n",
"\n",
"# First legend items\n",
"legend1_handles = []\n",
"\n",
"for ii, q in enumerate(q_list):\n",
" alpha = deflection_angle_x_inverse_correction(q, x_vals)\n",
" # Create the first legend handle for each q value\n",
" l1, = plt.plot([], [], c=col_list[ii], label='%.1f' % (q))\n",
" legend1_handles.append(l1)\n",
" plt.plot(x_vals, alpha, lw=3, c=col_list[ii], alpha=0.4)\n",
" plt.plot(x_vals, alpha, lw=1, ls='--', c=col_list[ii])\n",
"\n",
"# Second legend items\n",
"l2, = plt.plot([], [], c='k', lw=3, alpha=0.4, label=r'Eiroa et.al. upto $\\mathcal{O}(G)$')\n",
"l3, = plt.plot([], [], c='k', lw=1, ls='--', label='this work')\n",
"\n",
"# Get the current Axes object\n",
"ax = plt.gca()\n",
"\n",
"# Add the first legend\n",
"legend1 = ax.legend(handles=legend1_handles, loc='upper right', title=r\"$\\widetilde{Q}$\")\n",
"ax.add_artist(legend1) # Ensure the first legend stays on the plot\n",
"\n",
"# Add the second legend\n",
"ax.legend(handles=[l2, l3], loc='lower left')\n",
"\n",
"# Plot settings\n",
"plt.xscale('log')\n",
"plt.yscale('log')\n",
"plt.xlabel(r'$x$')\n",
"plt.ylabel(r'$\\alpha$')\n",
"plt.title(f\"Deflection angle for charged lens\")\n",
"plt.savefig('../plots/deflection_angle/ref_response.pdf',format='pdf', bbox_inches='tight')\n",
"plt.show()\n"
]
},
{
"cell_type": "code",
"execution_count": null, "execution_count": null,
"id": "319bbe99-9a2a-4c6e-8328-aa8c7c7fa6cc", "id": "93db4315-6fc1-4ec8-aca5-4cd611d091d9",
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [] "source": []
...@@ -209,7 +270,7 @@ ...@@ -209,7 +270,7 @@
"name": "python", "name": "python",
"nbconvert_exporter": "python", "nbconvert_exporter": "python",
"pygments_lexer": "ipython3", "pygments_lexer": "ipython3",
"version": "3.10.15" "version": "3.9.19"
} }
}, },
"nbformat": 4, "nbformat": 4,
......
notebooks/lens_class.ipynb
...@@ -142,7 +142,7 @@ ...@@ -142,7 +142,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 6, "execution_count": 5,
"id": "c30757f1-b4bc-4bcf-833b-7d7ed650bd0e", "id": "c30757f1-b4bc-4bcf-833b-7d7ed650bd0e",
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
...@@ -288,6 +288,150 @@ ...@@ -288,6 +288,150 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 6,
"id": "2935824b-e767-4763-b0e8-45e4ab20b559",
"metadata": {},
"outputs": [],
"source": [
"y_arr = np.linspace(0., 3., 100)\n",
"qeff_arr = [0.4]\n",
"\n",
"num_images = np.zeros((len(y_arr), len(qeff_arr)))\n",
"\n",
"for ii, y in enumerate(y_arr):\n",
" for jj, qeff in enumerate(qeff_arr):\n",
" lens = ChargedLens(y1=y, qeff=qeff)\n",
" ## using triangle mapping\n",
" # possible_images = np.array(find_images(lens, method='T')).T\n",
" # gtd1_im, gtd2_im = lens.gradTimeDelay(*possible_images.T)\n",
" # gtd_im = np.sqrt(gtd1_im**2 + gtd2_im**2)\n",
" # idxs = np.where(gtd_im < 1e-1)[0]\n",
" # x_image = possible_images[idxs]\n",
" # num_images[ii][jj] = len(x_image)\n",
" # num_images[ii][jj] = len(lens.imageLocations(mag_cut=True, gtd_cut=False, mag_tol=1e-6)[0])\n",
" num_images[ii][jj] = len(lens.img)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "0b896591-b871-4672-b4f3-06cbb4429c05",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[0.],\n",
" [0.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.],\n",
" [2.]])"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": []
},
{
"cell_type": "code",
"execution_count": 8, "execution_count": 8,
"id": "4ec7d737-bfb6-4722-befc-3c026be9bd83", "id": "4ec7d737-bfb6-4722-befc-3c026be9bd83",
"metadata": {}, "metadata": {},
......
notebooks/mismatch.ipynb
This source diff could not be displayed because it is too large. You can view the blob instead.
notebooks/plotting.ipynb
This source diff could not be displayed because it is too large. You can view the blob instead.
scripts/lensed_waveform_exec_hybrid.py 0 → 100755
#!/home1/uddeepta.deka/softwares/miniconda3/envs/igwn-py39/bin/python
__author__ = "U.Deka"
"""
This HTCondor compatible code is used to generate microlensed waveforms
and amplification factor in the time and frequency domain for the
charged lens model using contour integration method.
"""
from sys import argv
import sys
sys.path.append("../utils/")
import numpy as np
from lenses import ChargedLens
from generate_ampfac_charged_lens import Ftilde_hist, Ftilde_contour, ReconstructPeak, ampFac, Fw_geom
import time
import os
from scipy.interpolate import interp1d
from astropy.cosmology import FlatLambdaCDM, Planck18
import lal
import pycbc.psd as psd
from pycbc import detector
from pycbc.waveform import get_fd_waveform
import matplotlib.pyplot as plt
from plotsettings import *
from cycler import cycler
style = 'Notebook'
plotsettings(style)
plt.rcParams['axes.prop_cycle'] = cycler(color=dark2)
def Ftilde_hybrid(lens, dt_fac=5.0, wmax=100, t_max=100):
""" Computing F(t) using a hybrid contour+hist method. """
# dt = min(np.pi / wmax / dt_fac, lens.img_td / dt_fac)
dt = np.pi/wmax/dt_fac
# tcut = np.max(np.abs(np.diff(lens.img_td)))/10
tcut = min(0.5, np.max(np.abs(np.diff(lens.img_td)))/10)
t_hist, Ft_hist = Ftilde_hist(lens, 20*tcut, dt, dx_fac=1, verbose=False)
print("Hist done")
hybridize_at = tcut * 1.8
dt *= 10
t_contour, Ft_contour = Ftilde_contour(lens, t_max, dt, dx1=0.05, dx2=0.05, verbose=False, t_min = tcut)
print("Contour done")
Ft_hist_selected = Ft_hist[t_hist < hybridize_at]
Ft_contour_selected = np.array(Ft_contour)[t_contour >= hybridize_at]
Ft_hybrid = np.concatenate((Ft_hist_selected, Ft_contour_selected))
t_hybrid = np.concatenate((t_hist[t_hist < hybridize_at], t_contour[t_contour >= hybridize_at]))
Ft_new = Ft_hybrid
t_new = t_hybrid
WMAX = 200
DT_FAC = 50
DT = np.pi / WMAX / DT_FAC
num = int(2**np.ceil(np.log2((max(t_new) - min(t_new))/DT)))
t_hybrid = np.linspace(min(t_new), max(t_new), num)
Ft_hybrid = np.interp(t_hybrid, t_new, Ft_new)
return t_hybrid, Ft_hybrid
#===========================================================================
#============================ lensed waveform ==============================
#===========================================================================
class MicrolensedWaveform:
"""
Class that generates lensed (and unlensed) waveforms.
Parmeters:
----------
y : float, optional
Dimensionless impact parameter [default: 0.1]
q : float, optional
Dimensionless effective charge [default: 0.0]
Ml : float, optional
Lens mass in units of solar mass [default: 100.]
zl : float, optional
Lens redshift [default: 0.5]
dect : str, optional
Detector config. `Aplus` or `ET` [default: `ET`]
f_hi : float, optional
Maximum frequency (in Hz) in the detector band [default: 1024]
zs : float, optional
Lens redshift [default: 2.0]
ra : float, optional
Right ascension of the source [default: 1.7]
dec : float, optional
Declination of the source [default: 0.6]
pol : float, optional
Polarization angle of the source [default: 4.8]
m1 : float, optional
Mass of component 1 binary in units of solar mass [default: 20]
m2 : float, optional
Mass of component 2 binary in units of solar mass [default: 20]
chi1 : float, optional
Spin of component 1 BH along z-direction [default:0.]
chi2 : float, optional
Spin of component 2 BH along z-direction [default:0.]
iota : float, optional
Inclination angle in radians [default: PI/8]
phi0 : float, optional
Coalescence phase in radians [default: PI/4]
approx : str, optional
Waveform approximant [default: `IMRPhenomXP`]
delta_f : float, optional
Frequency step (in Hz) used to generate the waveform [default: 1/32]
t_max : float, optional
Maximum dimensionless time of integration for F(t) [default: 30]
dt_fac : float, optional
Factor to control the time step for F(t) [default: 7.5]
"""
def __init__(self, y=0.1, q=0., Ml=100., zl=0.5, dect='ET', f_hi=1024,
zs=2.0, m1=20, m2=20, chi1=0., chi2=0., iota=np.pi/8., phi0=np.pi/4.,
ra=1.7, dec=0.6, pol=4.8, approx='IMRPhenomXP', delta_f=1./32,
t_max = 30, dt_fac = 7.5):
"""
Class Initializer
"""
# Cosmology
cosmo = FlatLambdaCDM(H0=Planck18.H0.value, Om0=Planck18.Om0)
# Lens parameters
self.lens = ChargedLens(y1=y, qeff=q)
self.Ml = Ml
self.zl = zl
self.Mlz = self.Ml * (1 + self.zl)
# Detector
self.dect = dect
if self.dect=='Aplus':
psd_fname = '../data/PSD/AplusDesign_O5.txt'
self.f_lo = 18.
elif self.dect=='ET':
psd_fname = '../data/PSD/et_d.txt'
self.f_lo = 10.
else:
print("Detector unknown... Reverting to ET...")
psd_fname = '../data/PSD/et_d.txt'
self.dect = 'ET'
self.f_lo = 10.
self.f_Sh, self.dect_asd = np.loadtxt(psd_fname, unpack=True)
self.f_hi = f_hi
self.ra = ra
self.dec = dec
self.pol = pol
# Waveform
self.m1 = m1
self.m2 = m2
self.zs = zs
self.ldist_s = cosmo.luminosity_distance(self.zs).value
self.chi1 = chi1
self.chi2 = chi2
self.iota = iota
self.phi0 = phi0
self.approx = approx
self.delta_f = delta_f
self.f_final = 0.2 / (self.m1 + self.m2) / lal.MTSUN_SI
# Unlensed waveform
self.f_UL, self.hp_UL, self.hc_UL = self.get_unlensed_fd_waveform()
self.f_plus, self.f_cross, self.hh_UL, self.dect_psd = self.detector_setup()
# ampfac params
self.t_max = t_max
self.dt_fac = dt_fac
self.ts = 8 * np.pi * lal.MTSUN_SI * self.Mlz
self.t, self.Ft, self.Ft_recons, self.w, self.Ff, self.Ff_recons = self.get_lensed_ampfac()
self.Ff_geom = Fw_geom(self.lens, self.w)
# Lensed waveform
self.hh_L, self.hh_L_recons = self.get_lensed_fd_waveform()
def get_unlensed_fd_waveform(self):
"""
Returns unlensed plus- and cross-polarized waveforms
"""
hp, hc = get_fd_waveform(
approximant=self.approx, mass1=self.m1, mass2=self.m2,
distance=self.ldist_s, spin1z=self.chi1, spin2z=self.chi2,
inclination=self.iota, coa_phase=self.phi0, delta_f=self.delta_f,
f_lower=self.f_lo, f_final=self.f_final
)
hp = hp[: int(self.f_final / self.delta_f) + 1]
hc = hc[: int(self.f_final / self.delta_f) + 1]
f = hp.get_sample_frequencies().data
return f, hp, hc
def detector_setup(self):
"""
Returns frequency domain unlensed waveform and detector PSD
"""
if self.dect == 'Aplus':
inst = detector.Detector('H1') # change this to relevant detector
else:
inst = detector.Detector('V1')
f_plus, f_cross = inst.antenna_pattern(self.ra, self.dec, self.pol, 10.)
hh_u = f_plus * self.hp_UL + f_cross * self.hc_UL
psd_arr = self.dect_asd**2
dect_psd = psd.read.from_numpy_arrays(freq_data=self.f_Sh, noise_data=psd_arr,
length=len(self.hp_UL.data.data), delta_f=self.delta_f,
low_freq_cutoff=self.f_lo)
return f_plus, f_cross, hh_u, dect_psd
def get_lensed_ampfac(self):
"""
Computes and returns the time and frequency domain amplification factor
"""
wmin = self.ts * min(self.f_UL)
wmax = self.ts * max(self.f_UL)
dt = np.pi / wmax / self.dt_fac
t_vals, Ft_vals = Ftilde_hybrid(self.lens, dt_fac=self.dt_fac, wmax=wmax, t_max=self.t_max)
# t_vals, Ft_vals = Ftilde_contour(self.lens, self.t_max, dt,
# dx1 = 0.05, dx2 = 0.05, verbose=False)
# # dx1=0.025, dx2=0.025, verbose=False)
# interpolating
dt_ = np.pi / max(wmax, 200) / 50
num = int(2**np.ceil(np.log2((max(t_vals) - min(t_vals))/dt_)))
t_dense = np.linspace(min(t_vals), max(t_vals), num)
Ft_interped = np.interp(t_dense, t_vals, Ft_vals)
# applying peak correction
rp = ReconstructPeak(self.lens, t_dense, Ft_interped, peak_window_frac=0.05)
Ft_recons = rp.get_peak_reconstructed_Ft()
# doing the FFT
w_vals, Ff_vals = ampFac(t_dense, Ft_interped)
idxs = (w_vals > wmin) & (w_vals <= wmax)
w_vals = w_vals[idxs]
Ff_vals = Ff_vals[idxs]
_, Ff_recons = ampFac(t_dense, Ft_recons)
Ff_recons = Ff_recons[idxs]
f_vals = w_vals / self.ts
Ff_abs_interp = interp1d(f_vals, np.abs(Ff_vals), fill_value='extrapolate', kind='linear')
Ff_phase_interp = interp1d(f_vals, np.angle(Ff_vals), fill_value='extrapolate', kind='linear')
Ff_recons_abs_interp = interp1d(f_vals, np.abs(Ff_recons), fill_value='extrapolate', kind='linear')
Ff_recons_phase_interp = interp1d(f_vals, np.angle(Ff_recons), fill_value='extrapolate', kind='linear')
Ff_interped = Ff_abs_interp(self.f_UL) * np.exp(1j * Ff_phase_interp(self.f_UL))
Ff_recons_interped = Ff_recons_abs_interp(self.f_UL) * np.exp(1j * Ff_recons_phase_interp(self.f_UL))
return t_dense, Ft_interped, Ft_recons, self.ts*self.f_UL, Ff_interped, Ff_recons_interped
def get_lensed_fd_waveform(self):
""" Returns the lensed waveform """
lensed_hp, lensed_hc = self.Ff * self.hp_UL, self.Ff * self.hc_UL
lensed_hp = lensed_hp.cyclic_time_shift(3)
lensed_hc = lensed_hc.cyclic_time_shift(3)
lensed_hp_recons, lensed_hc_recons = self.Ff_recons * self.hp_UL, self.Ff_recons * self.hc_UL
lensed_hp_recons = lensed_hp_recons.cyclic_time_shift(3)
lensed_hc_recons = lensed_hc_recons.cyclic_time_shift(3)
hh_l = self.f_plus * lensed_hp + self.f_cross * lensed_hc
hh_l_recons = self.f_plus * lensed_hp_recons + self.f_cross * lensed_hc_recons
return hh_l, hh_l_recons
#===========================================================================
#================================ Plotting =================================
#===========================================================================
def make_plots(wf_obj, outdir=None):
max_td = np.max(wf_obj.lens.img_td)
x_max = 1.2 * np.sqrt(2 * max_td)
ngrid = 1000
x1_ = np.linspace(-x_max + wf_obj.lens.y1, x_max + wf_obj.lens.y1, ngrid, endpoint=True)
x2_ = np.linspace(-x_max + wf_obj.lens.y2, x_max + wf_obj.lens.y2, ngrid, endpoint=True)
X1, X2 = np.meshgrid(x1_, x2_)
timedelay_map = wf_obj.lens.timeDelay(X1, X2)
gradtd_map_X1, gradtd_map_X2 = wf_obj.lens.gradTimeDelay(X1, X2)
gradtd_map_norm = np.sqrt(gradtd_map_X1 ** 2 + gradtd_map_X2 ** 2)
fig, ax = plt.subplots(figsize=(figWidthsOneColDict[style],figHeightsDict[style]))
cs1 = ax.contourf(X1, X2, timedelay_map, 20, cmap='RdPu')
cs2 = ax.contour(X1, X2, gradtd_map_norm, levels=[0.0, 0.01, 0.1], colors=['r', 'r', 'r'])
ax.scatter(wf_obj.lens.y1, wf_obj.lens.y2, marker='^', label='source')
cbar = fig.colorbar(cs1, ax=ax)
ax.set_xlabel(r'$x_1$')
ax.set_ylabel(r'$x_2$')
for img in wf_obj.lens.img:
ax.scatter(img, 0, c='k', marker='+', s=30)
ax.set_xlim(np.min(wf_obj.lens.img)-0.5, np.max(wf_obj.lens.img)+0.5)
ax.set_ylim(-0.5, +0.5)
ax.plot([], [],ls='', marker='+', c='k', label='image')
ax.legend(loc='best', frameon=True)
plt.savefig(outdir+"verif_plot.pdf", format='pdf', bbox_inches='tight')
plt.close()
plt.figure(figsize=(figWidthsOneColDict[style], figHeightsDict[style]))
plt.plot(wf_obj.t, wf_obj.Ft, label='vanilla')
plt.plot(wf_obj.t, wf_obj.Ft_recons, label='reconstructed')
plt.axhline(1, ls='--', c='magenta')
image_time_delays = wf_obj.lens.img_td - np.min(wf_obj.lens.img_td)
for td in image_time_delays:
plt.axvline(td, ls='-', c='grey')
plt.plot([], [], ls='-', c='grey', label='images')
plt.xlabel(r'$t$')
plt.ylabel(r'$\widetilde{F}(t)$')
plt.xscale('log')
plt.legend()
plt.xlim(left=1e-2, right=wf_obj.t_max)
plt.savefig(outdir+"Ft_plot.pdf", format='pdf', bbox_inches='tight')
plt.close()
plt.figure(figsize=(figWidthsTwoColDict[style]*1.5, figHeightsDict[style]))
plt.subplot(131)
plt.plot(wf_obj.f_UL, np.abs(wf_obj.Ff), label='vanilla')
plt.plot(wf_obj.f_UL, np.abs(wf_obj.Ff_recons), label='reconstructed')
plt.xlabel(r'$f$[Hz]')
plt.ylabel(r'$|F(f)|$')
plt.xscale('log')
plt.legend()
plt.subplot(132)
plt.plot(wf_obj.f_UL, np.angle(wf_obj.Ff), label='vanilla')
plt.plot(wf_obj.f_UL, np.angle(wf_obj.Ff_recons), label='reconstructed')
plt.xlabel(r'$f$[Hz]')
plt.ylabel(r'arg($F(f)$)')
plt.xscale('log')
plt.legend()
plt.subplot(133)
plt.plot(wf_obj.f_UL, np.abs(wf_obj.hh_UL), lw=2, alpha=0.6, c='grey', label='unlensed')
plt.plot(wf_obj.f_UL, np.abs(wf_obj.hh_L), label='lensed')
plt.plot(wf_obj.f_UL, np.abs(wf_obj.hh_L_recons), label='lensed reconstructed')
plt.xlabel(r'$f$[Hz]')
plt.ylabel(r'$|h(f)|$')
plt.xlim(left=wf_obj.f_lo, right=min(max(wf_obj.f_UL), max(wf_obj.f_Sh)))
plt.xscale('log')
plt.yscale('log')
plt.legend()
plt.tight_layout()
plt.savefig(outdir+"strain_plot.pdf", format='pdf', bbox_inches='tight')
plt.close()
return
#===========================================================================
#=============================== Driver code ===============================
#===========================================================================
def usage():
print(
"Computes the lensed waveform due to charged lens."
)
print("\nIt uses the class 'MicrolensedWaveform' to generate results.")
print("The following options can be used to overwrite parameters.")
print("For example, to change the impact parameter to 0.5, run:")
print("python lensed_waveform_condor_exec.py y 0.5")
print(MicrolensedWaveform.__doc__)
print("\t folder_path: absolute path to the output directory")
def main():
print(" ---------------------------------------------------------------- ")
print(" :: lensed_waveform_exec.py :: use flag -h for list of options :: ")
print(" ---------------------------------------------------------------- ")
kwargs = {"y":0.1, "q":0.0, "Ml":100.0, "zl":0.5, "dect":"ET",
"f_hi":1024.0, "zs":2.0, "m1":20.0, "m2":20.0, "chi1":0.0, "chi2":0.0,
"iota":np.pi/8., "phi0":np.pi/4., "ra":1.7, "dec":0.6, "pol":4.8,
"approx":"IMRPhenomXP", "delta_f":1./32, "t_max":30., "dt_fac":7.5,}
folder_path = "/home1/uddeepta.deka/lensing/charged_lens_updated/data/lensed_waveform_data/results_hybrid/"
for i in range(len(sys.argv)):
if sys.argv[i] == "-h":
usage()
return
if sys.argv[i] == "folder_path":
folder_path = str(sys.argv[i+1])
if sys.argv[i] in kwargs:
if (sys.argv[i]=="dect") or (sys.argv[i]=="approx"):
kwargs[sys.argv[i]] = str(sys.argv[i+1])
else:
kwargs[sys.argv[i]] = float(sys.argv[i+1])
tStart = time.perf_counter()
mw = MicrolensedWaveform(**kwargs)
runtag = 'y_%.3f_q_%.1e_Ml_%.1f_zl_%.1f_dect_%s_fhi_%.1f_zs_%.1f_m1_%.1f_m2_%.1f_approx_%s_deltaf_%.5f_tmax_%.1f_dtfac_%.1f'%(
mw.lens.y1, mw.lens.qeff, mw.Ml, mw.zl, mw.dect, mw.f_hi, mw.zs, mw.m1, mw.m2, mw.approx, mw.delta_f, mw.t_max, mw.dt_fac)
outdir = folder_path + runtag + "/"
isExist = os.path.exists(outdir)
if not isExist:
os.makedirs(outdir)
print("Created new folder to store data.")
## Saving data to files ...
Ft_datafile = outdir + 'Ft_data.txt'
waveform_datafile = outdir + 'h_lensed_data.npz'
header = f"Created using lensed_waveform_exec_test.py.\nUsage: \nimport numpy as np\ndata = np.loadtxt('{Ft_datafile}')\nt = data[:,0]\nFt = data[:,1]\nFt_recons = data[:,2]\n"
with open(Ft_datafile, 'w') as f:
np.savetxt(f, np.array([mw.t, mw.Ft, mw.Ft_recons]).T, header=header)
np.savez(waveform_datafile, w=mw.w, f=mw.f_UL, Ff=mw.Ff, Ff_recons=mw.Ff_recons,
h_ul=mw.hh_UL, h_l=mw.hh_L, h_l_recons=mw.hh_L_recons)
make_plots(mw, outdir=outdir)
tEnd = time.perf_counter()
print("Data stored in : ", outdir)
print("Time taken to generate data : %.2f s."%(tEnd-tStart))
return
if __name__ == "__main__":
main()
\ No newline at end of file
scripts/submit_lensed_waveform_charged_lens.sub
Universe = vanilla Universe = vanilla
Executable = /home1/uddeepta.deka/lensing/charged_lens_updated/scripts/lensed_waveform_exec.py Executable = /home1/uddeepta.deka/lensing/charged_lens_updated/scripts/lensed_waveform_exec.py
Arguments = Ml $(Ml_value) y $(y_value) q $(q_value) t_max 200 Arguments = Ml $(Ml_value) y $(y_value) q $(q_value) t_max 200 dect Aplus
Getenv = True Getenv = True
Output = /home1/uddeepta.deka/lensing/charged_lens_updated/data/lensed_waveform_data/outfiles/output_Ml$(Ml_value)_y$(y_value)_q$(q_value).out Output = /home1/uddeepta.deka/lensing/charged_lens_updated/data/lensed_waveform_data/outfiles/output_Ml$(Ml_value)_y$(y_value)_q$(q_value).out
Error = /home1/uddeepta.deka/lensing/charged_lens_updated/data/lensed_waveform_data/errors/error_Ml$(Ml_value)_y$(y_value)_q$(q_value).err Error = /home1/uddeepta.deka/lensing/charged_lens_updated/data/lensed_waveform_data/errors/error_Ml$(Ml_value)_y$(y_value)_q$(q_value).err
Log = /home1/uddeepta.deka/lensing/charged_lens_updated/data/lensed_waveform_data/logs/log_Ml$(Ml_value)_y$(y_value)_q$(q_value).log Log = /home1/uddeepta.deka/lensing/charged_lens_updated/data/lensed_waveform_data/logs/log_Ml$(Ml_value)_y$(y_value)_q$(q_value).log
notify_user = uddeepta.deka@icts.res.in notify_user = uddeepta.deka@icts.res.in
Notification = Never Notification = Never
requirements = (Machine != "node42.alice.icts.res.in") && (Machine != "node51.alice.icts.res.in") && (Machine != "node76.alice.icts.res.in") requirements = (Machine != "node18.alice.icts.res.in") && (Machine != "node43.alice.icts.res.in") && (Machine != "node50.alice.icts.res.in") && (Machine != "node63.alice.icts.res.in") && (Machine != "node72.alice.icts.res.in") && (Machine != "node76.alice.icts.res.in") && (Machine != "node82.alice.icts.res.in") && (Machine != "node94.alice.icts.res.in")
request_cpus = 4 max_idle = 100
request_memory = 1024 request_cpus = 1
request_memory = 2048
queue Ml_value y_value q_value from params.txt queue Ml_value y_value q_value from params_pos.txt
\ No newline at end of file \ No newline at end of file
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or sign in to comment