{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/me-manu/gammaALPs/blob/master/docs/tutorials/mixing_single_cell.ipynb)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Mixing in a homogeneous magnetic field\n", "\n", "This tutorial demonstrates the most simple use case of the `gammaALPs` code, namely the mixing between photons and ALPs in a homogeneous magnetic field. The magnetic-field parameters are probably not very realisitic, but the results nicely illustrate some of the main features of photon-ALP mixing. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you haven't installed `gammaALPs` already, you can do so using `pip`. Just uncomment the line below:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "#!pip install gammaALPs" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We first import the relevant modules:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "from gammaALPs.core import Source, ALP, ModuleList\n", "from gammaALPs.base.transfer import EminGeV, EmaxGeV\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from astropy import constants as c" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And set the ALP parameters:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "m, g = 10., 3.\n", "alp = ALP(m, g)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's define an energy range for this example, in GeV." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "EGeV = np.logspace(0.,8.,1000)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For the initial polarization, we use a fully polarized photon beam," ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "pin = np.diag((1., 0., 0.))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And we define some dummy source:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "src = Source(z=0., l=0., b=0.)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With the defined `alp`, `src`, initial polarization `pin` and energy range `EGeV`, we can now initialize our module list" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "ml = ModuleList(alp, src, pin=pin, EGeV=EGeV, seed=0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We add a single propagation enivironment, mixing in a cell-like field in a galaxy cluster. \n", "By setting `eta` and `beta` to zero, we force the electron density to be constant, so \n", "that neither the magnetic field nor the electron density change with propagation distance. \n", "We choose arbitrary values for B-field strength and electron density that give nice results for illustration. " ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "\u001b[0;36menvirons.py:\u001b[0;35m 257\u001b[0;0m --- \u001b[1;36mINFO\u001b[1;0m: Using inputted chi\n", "\u001b[0;36menvirons.py:\u001b[0;35m 270\u001b[0;0m --- \u001b[1;31mWARNING\u001b[1;0m: r_abell <= L0: assuming one domain from 0. to L0\n" ] } ], "source": [ "ml.add_propagation(environ='ICMCell',\n", " order=0, # order of the module\n", " B0=1., # B field strength\n", " L0=10., # cell size\n", " nsim=1, # one single realization\n", " n0=1e-3, # electron density\n", " r_abell=10.001, # full path, chosen that we only have a single cell\n", " beta=0., \n", " eta=0.\n", " )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we explicitly set the angle $\\psi$ between the transversal magnetic field and the $x$ direction to $\\pi / 2$ so that the $x$ polarization will fully mix. " ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "ml.modules[0].psin = np.ones_like(ml.modules[0].psin) * np.pi / 2." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We run the photon-ALP calculation to get the final oscillation probabilities into the two photon polarization states and into the ALP state" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "\u001b[0;36m core.py:\u001b[0;35m 652\u001b[0;0m --- \u001b[1;36mINFO\u001b[1;0m: Running Module 0: \n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "(1, 1000)\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/Users/manuelmeyer/Python/gammaALPs/gammaALPs/base/transfer.py:799: UserWarning: Not all values of linear polarization are real values!\n", " warnings.warn(\"Not all values of linear polarization are real values!\")\n", "/Users/manuelmeyer/Python/gammaALPs/gammaALPs/base/transfer.py:802: UserWarning: Not all values of circular polarization are real values!\n", " warnings.warn(\"Not all values of circular polarization are real values!\")\n" ] } ], "source": [ "px, py, pa = ml.run()\n", "print (pa.shape)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And plot the results:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0, 0.5, '$P_{a\\\\gamma}$')" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEKCAYAAADXdbjqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAA1oklEQVR4nO3deXxc5Xno8d8zo937IttavBKDcYJtLLFmaUhKCtQLWWgguYSkfLDJDUnpbW7jprfFbu/tpVma2/RSwE5ISG+BkoXEdiEkISEhC2DJOMbGGIxXLZbkTbs0mpnn/jFn5NF4JM12zkia5/v56DNzznnPmfeMzswz73ZeUVWMMcaYTPhynQFjjDHjnwUTY4wxGbNgYowxJmMWTIwxxmTMgokxxpiMWTAxxhiTsYJcZyAXZs+erYsWLcp1NowxWdDU1ERlZWWus5EX6uvrT6lqeaJteRlMFi1aRF1dXa6zMaHV19dTU1OT62yYPCAiNDU15TobeUFEjg23zdNqLhG5QUQOisghEdmUYPvHRWSv8/dbEVk52r4iMlNEfioibzqPM7w6HzO82traXGfBGOMhz4KJiPiBB4AbgeXAbSKyPC7ZEeAPVHUF8PfA1iT23QQ8p6pLgeecZWOMMR7ysmRyJXBIVQ+ragB4Algfm0BVf6uqZ53FF4HqJPZdDzzqPH8UuNm9UzDGGJOIl8GkCjgRs9zgrBvOncAzSew7V1WbAZzHOYkOJiIbRKROROra2trSyL5JxX333ZfrLJg8Ydfa2OBlMJEE6xLeZVJEriMSTL6Q6r7DUdWtqlqrqrXl5Qk7I5gs2rx5c66zYPKEXWtjg5fBpAGYH7NcDVzQBUNEVgDfANar6ukk9m0RkQpn3wqgNcv5NmmwrprGK3atjQ1edg3eBSwVkcVAI3Ar8LHYBCKyAPgBcLuqvpHkvtuBO4D7nccfuXkSJjnNzc2evI6q0hMI0d47QO9AiD7nrzcQjjwOhBgIhQmGlVCiP408BkOR56gOFnmjszNE15xfvnBb3APRqR00bn1UtKgtEn2U8+uj65wng2mGPD+/bbDYHnuMwW0XHmNwmwwt8MceO3rcAr+PQr9Q6Pc5f0JR9HlBZLm4wMfk4kKmlBQwpaSASUUF+HyJKhPcYd38xwbPgomqBkXkHuBZwA88oqr7ReRuZ/tDwN8Cs4B/dS70oFM1lXBf59D3A0+KyJ3AceAWr87JuK+jb4DXmzs5drqbhrO9NJztpelcL2d7ApG/7gECoXCus2liiMDkogKmlRUyb2oJc52/yuklXFQ+mbfNmUzV9NKsBZz6+nornYwBko+TY9XW1qr9mnFXTU0N9fX1Ke2jqhw+1c1vDp3it4dOs6+pnYazvaPuV1LoY1ppIWVFBZQU+ikp9FFa6Kek0E9poZ9Cv+D3+SjwCX6/4BfB74v8FfgEX/RRIn+Q+Jd85PH8F2CiUkLCfeK2p1Liif14JirtjHSs6EKiklZ8iWnIcZznYVWCYWUgFGYgGHkMhMKR5ZCzPhSmbyBMV3+Qzr4BOvuC9ARCjKa00M87qqZyxaKZXLF4JtcsmUVJoX/U/RIREfLxeywXRKReVRMOIrNgYnLurbYufvRKIz/c08TxMz1DthUV+Lhk7hQuKp9E9Ywy5s8spXJ6KTMnFTFzUhEzyorS/hIy7giGwnT3hzjTE+Bkex8tHZG/xnO9HGrt4s3WLto6+4fsU1bk57pL5rB2ZSV/eOkcCvzJN+daMPHOSMEkL2+nYty3YcMGtm7dOux2VeXlI2f41+ff4pdvnO+qPWtSEde+bTbvvGgWNQtnsHj2pJS+WEzuFfh9TCvzMa2skMWzJyVMc7Y7wO7jZ3n56Bl+c+gU+xo7+M9Xm/nPV5upnFbCJ65dxB3XLKK0yH4ojBdWMjGuGOnX4pFT3WzZsZ/nD0aCSGmhnzUrKvjg5VVctWQWfg8bb83Y0HC2h2dePcljLx/nyKluAOZNLeG/feBibqmpvqCzQCwrmXjHqrniWDBxX6IPuKryyG+O8o/PvE4gFGZKcQGfetdiPnntImZOKspRTs1YEg4rv3yzja/+5CD7GjsAePfS2Xz5IyuZN60k4T4WTLxjwSSOBRP3xX/Au/qD3PvEHn52oAWAj9RU84UbllE+pThXWTRjWDis/HBPI3+38zXO9Qwwo6yQrZ+o5YpFMy9Ia8HEOyMFE6uMNq5obGwcfH66q5+PbXuRnx1oYWpJAQ/fXsNXbllpgcQMy+cTPrS6mp/c+x7evXQ2Z3sG+Pi2l9j++wtvNR97rZncsWBiXBHtFnymO8BHt77I3oZ2FswsY8dn38UfvX1ejnNnxos5U0v41iev4I5rFhIIhbn3iVfYERdQUu2CbtxhwcS4Yt26dfQEgvzpt3dxqLWLi+dO5nt3X8PCWYl79xgznAK/jy3r38Gf/+HFhBX+/D/28KuYHoAj9Ro03rFgYlzzl9/by54T56iaXsq/3XkVc6YmbkA1Jhmfe//b2PieJQTDymcff4VjpyO9vnbs2JHjnBmwYGJctHNvM5OK/Dz6p1cy1wKJyZCI8IUblvGHl86hvXeAjf9WT38wxNq1a3OdNYMFE+OCt9q6mHPjZwH4hw9dxtvmTM5xjsxE4fMJ//TRVSycVcbrJzv5+nNvsnPnzlxny2DBxGSZqvLXT71K6Yo/4sOrq1m/aqT5z4xJ3dSSQr56y0pE4MHn38p1dozDgonJqu/VN/Di4TMc+8c1/I8/vjTX2TETVO2imdz5zsWEozeotHEmOWfBxGRN30CIr/zk4ODyDBvVblx07/UXM3tyZKxSovEnxlsWTEzWPPrbo7R09PP2yqm5zorJA5OLC/j8By4G4Es/PkggaPPa5JIFE5MVnX0D/KtTf/2XNyxjzZo1Oc6RyQe31M5n9vJraDzXyw9fsZHwueRpMBGRG0TkoIgcEpFNCbYvE5HfiUi/iHw+Zv0lIrIn5q9DRO51tm0WkcaYbTd5eErG8R+7TtDeO8CVi2bynqWzre+/8YTfJ2z7f98F4MFfvkUobG0nueJZMBERP/AAcCOwHLhNRJbHJTsDfA74SuxKVT2oqqtUdRVQA/QAT8Uk+Vp0u6o+7dY5mMSCoTDf+s1RADa8ZwkiYn3/jWe+8TcbmT+zlCOnunlmX3Ous5O3vCyZXAkcUtXDqhoAngDWxyZQ1VZV3QUMjHCc9wNvqeox97JqUvHj/SdpPNfLktmTeN+yOQDW9994ZuPGjWx49xIAvvM7+1rIFS+DSRVwIma5wVmXqluBx+PW3SMie0XkERGZkWgnEdkgInUiUtfW1pYoiUlT9AP8qXctxmcTWxmP1dTUcPPlVZQV+Xn5yBkOtXbmOkt5yctgkuhbJqUKThEpAtYB341Z/SBwEbAKaAa+mmhfVd2qqrWqWlteXp7Ky5oRHDvdzctHzlBa6OeDl9sAReO9qqoqppQUsn5VJQCPvXRilD2MG7wMJg3A/JjlaiDVzuE3ArtVtSW6QlVbVDWkqmFgG5HqNOOR79U3AHDjZfOYXFwwuN4GkRmvfezKhQB8f3cD/cFQjnOTf7wMJruApSKy2Clh3ApsT/EYtxFXxSUiFTGLHwT2ZZRLk7RQWPm+E0xuqZk/ZJvdFtx47bLqaVxaMZX23gF+9capXGcn73gWTFQ1CNwDPAscAJ5U1f0icreI3A0gIvNEpAH4b8D/EJEGEZnqbCsDrgd+EHfoL4nIqyKyF7gO+HOPTinv7Tp6hqb2PqpnlHLV4qHTqW7cuDFHuTL5bM2KyG/LnXttRLzXCkZPkj1Ot92n49Y9FPP8JJHqr0T79gCzEqy/PcvZNEn68b6TAKxZUWkN72ZMWLuiki8/e5CfvdZCbyBEaZE/11nKGzYC3qRFVXl2fySY3PAOm4bXjA0LZpWxsnoa3YEQvzjYmuvs5BULJiYtexvaaW7vo2JaCSuqpl2wffv2VJvDjElP/LX2x05VV/THjvGGBROTlmecKq4/evu8hFVcNTU1XmfJ5Kn4a+39l84F4JdvtNntVTxkwcSk5RevR6oQPrB8bsLtVVU25sR4o7a2dsjyktmTWDCzjHM9A+w5cTZHuco/FkxMylo6+jjY0klZkZ+aRQlvOGCMZ5qahvbcEpHB2/r84nW724VXLJiYlL3wZqQP/9VLZlFcYL1lTG5t3rz5gnXvvSRylwtrhPeOBROTshfejPzae/fS2cOmueuuu7zKjslzW7ZsuWDd1UtmUVLoY39TB62dfTnIVf6xYGJSEg4rv3ZKJu9eOvw9zmwEvMmlkkI/VyyKDKR9+ciZHOcmP1gwMSk5cLKD090BKqeVcFH5pGHTWW8uk2vRuzK8ePh0jnOSHyyYmJREf+VdvWQWIsOPet+9e7dXWTImoauXRG6Y8dJhK5l4wYKJSUnd0UhXyyvi7sVlzFizono6JYU+3mzt4lRXf66zM+FZMDFJU1V2HY38yrtilC7BFRUVI243JluGu9aKCnzULIxcp9Zu4j4LJiZpJ8700trZz4yyQi4qnzxi2vi+/8a4ZaRr7erF0aouazdxmwUTk7S6Y5FfdzULZ47YXgKJ+/4b44aRrrVoyWTPiXPeZCaPWTAxSdsVbS9JYtR7or7/xnjtsuppiMBrzR30Ddjsi26yYGKS9srxSDCJ/tozZiwYqWQypaSQpXMmMxBSXmvu8C5TecjTYCIiN4jIQRE5JCKbEmxfJiK/E5F+Efl83LajzoyKe0SkLmb9TBH5qYi86TzaN50LegMh3mjpxCfw9soLbzlvTK5UVlaOuH1l9XQAfm9VXa7yLJiIiB94ALgRWA7cJiLL45KdAT4HfGWYw1ynqqtUNfY2oZuA51R1KfCcs2yy7LXmDsIKF8+dktTsdXV1daOmMSYbmpubR9y+asF0wNpN3OZlyeRK4JCqHlbVAPAEsD42gaq2quouYCCF464HHnWePwrcnIW8mjivNpwD4LIEE2EZM5atmj8dsGDiNi+DSRVwIma5wVmXLAV+IiL1IrIhZv1cVW0GcB7nZJxTc4G9je0ArKhOLpjEzzFhTK5cMncKJYU+jp3u4Wx3INfZmbC8DCaJ+pKmMg3aO1V1NZFqss+IyHtSenGRDSJSJyJ1bW02x0Gq9jnB5B1WMjHjTIHfx/KKqQDWCO8iL4NJAzA/ZrkaSHpkm6o2OY+twFNEqs0AWkSkAsB5TDiBgapuVdVaVa0tLx/+brfmQj2BIIdauyjwCZc6H0pjxpPodXvAgolrvAwmu4ClIrJYRIqAW4HtyewoIpNEZEr0OfABYJ+zeTtwh/P8DuBHWc214bWm843vJYXJTYZ13333uZwrYyKSudYutZKJ6wq8eiFVDYrIPcCzgB94RFX3i8jdzvaHRGQeUAdMBcIici+Rnl+zgaecUdcFwGOq+mPn0PcDT4rIncBx4BavzilfRH/NLa9MvlRiI+CNV5K51gaDSZMFE7d4FkwAVPVp4Om4dQ/FPD9JpPorXgewcphjngben8VsmjgHWzoBWDZvStL7VFZW2v25jCeSudaWzZuCCLzV1kUgGKaowMZrZ5u9o2ZUB09GgsklKQST0fr+G5MtyYxpmlRcwKJZkxgIKYdauzzIVf6xYGJGpKq8nkYwMcYr9fX1SaW7tCJy/VojvDssmJgRNbf30dkXZOakIsonFye93+rVq13MlTHnrVu3Lql0l86zRng3WTAxI4pWcV08d/Kot52PleyvRWO8ssxphH/DaQM02WXBxIwoWsW1bF5q40s2bNgweiJjPPS2OZEJ3azNxB0WTMyIDp6MVAmk2l6ybds2N7JjTNrmzyilyO9zqm5Tuf2fSYYFEzOiN1oiv+IunmuN72Z8K/D7WFI+CYC32rpznJuJx4KJGVY4rBw+FQkm0SoCY8azi6yqyzUWTMywTnb00TcQZvbkIqaVFqa0b2Njo0u5MmaoVK61t5VbMHGLBRMzrMNOVcCS2amXSqw3l/FKKtfa0rnRYGI9urLNgokZVrSKK1rPnIpk+/4bk6mtW7cmndZ6dLnHgokZVrRksnh26sHEGK/s2LEj6bSLZ0/CJ3D8TA99AyEXc5V/LJiYYb3VFi2ZWOO7GbvWrl2bdNriAj8LZpYRVjh62np0ZZMFEzOswTaTNKq5Hn744Wxnx5iEdu7cmVL6aEn76KkeN7KTtyyYmIT6BkI0tfdS4BMWzCxLeX8bAW/GqoWznGBiJZOssmBiEjp6uhtVWDCzjEJ/6pdJKvfxMsZL0ZLJMQsmWeVpMBGRG0TkoIgcEpFNCbYvE5HfiUi/iHw+Zv18EfmFiBwQkf0i8mcx2zaLSKOI7HH+bvLqfCayTKq4jBnLFs6KlLStmiu7PJtpUUT8wAPA9UADsEtEtqvqazHJzgCfA26O2z0I/IWq7nbmgq8XkZ/G7Ps1Vf2Ku2eQXw47je/Wk8tMNINtJlYyySovSyZXAodU9bCqBoAngPWxCVS1VVV3AQNx65tVdbfzvBM4AFR5k+38dPhUtGSSXk+uNWvWZDM7xgwr1WutanopBT6hub3PugdnkZfBpAo4EbPcQBoBQUQWAZcDL8WsvkdE9orIIyIyI6NcGgBOnIlUASxMo/EdUuv7b0wmUr3WCvw+5jvX9bHTVtWVLV4Gk0QtsprSAUQmA98H7lXV6HRpDwIXAauAZuCrw+y7QUTqRKSura0tlZfNS8edYDI/zWCSSt9/YzKRzrU22G5iVV1Z42UwaQDmxyxXA03J7iwihUQCyb+r6g+i61W1RVVDqhoGthGpTruAqm5V1VpVrS0vL0/rBPJF30CIlo5+CnxCxbSStI6Rat9/Y9KVTjf0RbOsR1e2eRlMdgFLRWSxiBQBtwLbk9lRIv1MvwkcUNV/ittWEbP4QWBflvKbtxrO9gJQOb2UgjS6BRvjpZqampT3WeSUTI5Yj66s8aw3l6oGReQe4FnADzyiqvtF5G5n+0MiMg+oA6YCYRG5F1gOrABuB14VkT3OIb+oqk8DXxKRVUSqzI4CG706p4kq2l6SzmBFY7xWVVWFako15ix0enQdP2Mlk2zxLJgAOF/+T8eteyjm+Uki1V/xfk3iNhdU9fZs5tFk3l4CpPzhNsZL82dEru1oKdxkzuowzAXOB5PStI+Rym3BjfFa9YzItd10rpdQ2H74ZIMFE3OBbFRzbdxotY1m7Cop9DN7cjEDIaW1sy/X2ZkQLJiYCxy3NhOTB6KlkxNnrKorGyyYmCFU1RrgTV6Itgk2nLUeXdlgwcQMcaY7QHcgxJTiAqaVFqZ9nO3bk+r1bUzG0r3WoiUTa4TPDgsmZogTzgdr/syyjG4jn07ff2PSke61dj6YWMkkGyyYmCGy1V5SVWX34TTeqK2tTWu/aqd7sLWZZIcFEzPEiSx0CzbGS01NSd+VaYjBksk5K5lkgwUTM0TTufO3UjFmPNi8eXNa+1U513jzuT6CoXAWc5SfLJiYIbIVTO66665sZMeYUW3ZsiWt/UoK/cyZUkwwrJzssLEmmbJgYoZobo98qKoyDCY2At6MB9ajK3ssmJghGp2SSbq3no+y3lxmPKi2e3RljQUTM6izb4DOviDFBT5mTirK6Fi7d+/OUq6McU+0o0m044lJnwUTMyi2iiuTMSbGjBcV05xG+HYrmWTKgokZ1JjFnlwVFRWjJzImCzK51iqnR6pzoz+kTPosmJhBTVlqL4H0+/4bk6pMrrXzJRMLJplKK5iIyCXZzojJveZzkQ9UNkom6fb9NyZVmVxr0R9Ozed6bUK3DKVbMrlNRDaJSEozNYrIDSJyUEQOicimBNuXicjvRKRfRD6fzL4iMlNEfioibzqPM9I8p7wXLZlk2i0Y0u/7b4yXppUWUlropzsQorM/mOvsjGtpBRNV3QzsAB4UkeskidZaEfEDDwA3EpnX/TYRWR6X7AzwOeArKey7CXhOVZcCzznLJg1NTiNkxfTMq7mM8UomJRMRGbzeoyVzk550q7m2AH9OZF729wEPJ7HblcAhVT2sqgHgCWB9bAJVbVXVXcBACvuuBx51nj8K3Jz6GRmApixWcxnjlcrKyoz2j1Z1NVmProykVE0V4xeq+nyK+1QBJ2KWG4CrsrDvXFVtBlDVZhGZk+gAIrIB2ACwYMGCFLKdH8JhHeweWTkt82BSV1eX8TGMSUZzc3NG+0cb4U9aI3xG0m0zea+IPCMi3xCRzyS5T6KqsGRbvDLZN5JYdauq1qpqbXl5eSq75oVT3f0MhJQZZYWUFvlznR1jPFMZ0whv0pduMJkOvAj8LyDZnl0NwPyY5Wog2T59I+3bIiIVAM5ja5LHNDGyXcWV7hwTxnhtnlMyabKSSUbSDSZnAD+RL+4zSe6zC1gqIotFpAi4FUh2vs2R9t0O3OE8vwP4UZLHNDHs1vMmX0Ub4K2aKzNptZmo6t+JSCXwdWBfkvsEReQe4FkigegRVd0vInc72x8SkXlAHTAVCIvIvcByVe1ItK9z6PuBJ0XkTuA4cEs655Tvstkt2JjxpHKwZGLVXJkYNZiIyB3AV4mUYnYCn1HVTlVtAu5M5cVU9Wng6bh1D8U8P0mkCiupfZ31p4H3p5IPc6FoNVc2Rr8D3HfffVk5jjGjyfRai+0arKp2X7o0JVPN9TfA9cAy4BjwD67myOREizM50LwsBRMbAW+8kum1NqW4gElFfnoHQnT02sDFdCUTTDpU9RVnDMjfEBnzYSaY6Exz86ZmJ5hk2vffmGRleq1FBi5aVVemkgkmFSKyQUTeLSLlQKHbmTLei5ZM5mYpmGTa99+YZGVjTFO0etca4dOXTAP8fcAK4OPAZcBkEXka+D2wV1UfdzF/xgOqSmtHP5C9YGKMV+rr620U/BgwajBR1SGTeYtINZHgchlwE2DBZJw72zNAIBRmaklB1gYsrl69OivHMWY069aty/iOv4O3orf7c6Ut5a7BqtpAZBDhBT2rzPiU7SouiPxaNGa8iE6SZSWT9NnkWOZ843uWenIBbNiwIWvHMsZt0R9S0epekzoLJoZWJ5jMmZK9YLJt27asHcsYt0WDSbSUblJnwcTQMtj4XpzjnBiTGxZMMmfBxLhSzWXMeDKjrJAiv4+OviC9gVCuszMuWTAxrlRzNTY2Zu1YxowkG9eaiDDHKZm3dlrpJB0WTIwrJRPrzWW8kq1rLVrVZQMX02PBxLjSZrJu3bqsHcuYkWzdunX0REmIXv8tndajKx0WTPJcMBTmVFc/IjB7sjXAm/Fnx44dWTlOtJq31Rrh02LBJM+1dfWjGgkkhX67HMz4s3bt2qwcx3p0Zca+PfKcW92CH3744awez5jh7Ny5MyvHGazmsoGLafE0mIjIDSJyUEQOicimBNtFRL7ubN8rIqud9ZeIyJ6Yvw5nFkZEZLOINMZsu8nLcxrvoo2N2br1fJSNgDfjzTwrmWQkrWl70yEifuABIhNtNQC7RGS7qr4Wk+xGYKnzdxXwIHCVqh4EVsUcpxF4Kma/r6nqV1w/iQko2g1yTpaDiYhkfPM9Y7wU/Qy0WgN8WrwsmVwJHFLVw6oaAJ4A1selWQ98RyNeBKaLSEVcmvcDb6nqMfezPPEN3uQxi2NMjBmPotVcJ9v77IdQGrwMJlXAiZjlBmddqmlu5cLb3t/jVIs9IiIzEr24M8FXnYjUtbW1pZ77Cepke+RX2Lxp1pPL5LfJxQWUOdP3dvbb9L2p8jKYSIJ18eF/xDQiUgSsA74bs/1B4CIi1WDNwFcTvbiqblXVWlWtLS8vTyHbE5tb1Vxr1qzJ6vGMGU62rjURibl7sLWbpMrLYNIAzI9ZrgaaUkxzI7BbVVuiK1S1RVVDqhoGtmFz1KfErQb4bPX9N2Y02bzW5kyxHl3p8jKY7AKWiship4RxK7A9Ls124BNOr66rgXZVjZ1M/Dbiqrji2lQ+COzLftYnLjcmxoLs9f03ZjTZvNaitxSyHl2p86w3l6oGReQe4FnADzyiqvtF5G5n+0NEZm+8CTgE9ACfiu4vImVEeoJtjDv0l0RkFZHqsKMJtpth9AZCdPQFKfL7mFFWmNVjZ6vvvzGjyWY39PMDF61kkirPggmAqj5N3HS/ThCJPlfgM8Ps2wPMSrD+9ixnM29Ef33NmVqMSKLmKmPGvpqamqwd63w1l5VMUmUj4POYW1Vcxnipqiq+w2f67JYq6bNgkscGbz3vQjCxfvpmPLJgkj4LJnms1akXnuPCdL3Zui24MV6y+3Olz4JJHmtxsWSycaP1gzDjz+A4k04bBZ8qCyZ57KS1mRgzREmhn2mlhQyElLM9A7nOzrhiwSSPuVnNZcx4db6qy9pNUmHBJI+52QC/fXv8eFRj3JHta80a4dNjwSRPqaqrXYOz2fffmJFk+1qLTt9rwSQ1FkzyVEdvkP5gmCnFBUwqzv7Y1Wz2/TdmJLW1tVk9nvXoSo8Fkzx1Mmb0uzHjWVNT/P1iM2PVXOmxYJKnbPS7mSg2b96c1eNZA3x6LJjkKTfHmADcddddrhzXmHhbtmzJ6vHsZo/psWCSp87f5NGdYGIj4M14ZdVc6bFgkqeiv7rmutRmYr25zHhVPqUYETjV1U8wFM51dsYNCyZ5ys0xJgC7d+925bjGuK3Q72PWpGLCCqe6ArnOzrhhwSRPtbpczWXMeDZvWqTEftKqupLmaTARkRtE5KCIHBKRTQm2i4h83dm+V0RWx2w7KiKvisgeEamLWT9TRH4qIm86jzO8Op/xzO1qroqKitETGZMFblxr0RL7yXYLJsnyLJiIiB94ALgRWA7cJiLL45LdCCx1/jYAD8Ztv05VV6lq7CilTcBzqroUeM5ZNiMIhZW2Lue+XFPcKZlku++/McNx41qbE3P3YJMcL0smVwKHVPWwqgaAJ4D1cWnWA9/RiBeB6SIy2s+O9cCjzvNHgZuzmOcJ6XRXP6GwMmtSEUUF7lwC2e77b8xw3LjWrGSSOi+DSRVwIma5wVmXbBoFfiIi9SKyISbNXFVtBnAe52Q11xPQ+Sou99pLst333xgv2S1VUpf9mzINTxKsi599ZqQ071TVJhGZA/xURF5X1V8l/eKRALQBYMGCBcnuNiGdn8fEbqVixj83SiY21iR1XpZMGoD5McvVQHxl57BpVDX62Ao8RaTaDKAlWhXmPLYmenFV3aqqtapaW15enuGpjG92KxUzkVRWVmb9mBZMUudlMNkFLBWRxSJSBNwKxE9EsB34hNOr62qgXVWbRWSSiEwBEJFJwAeAfTH73OE8vwP4kdsnMt61ehBM6urqRk9kTBY0Nzdn/ZiDbSYWTJLmWTWXqgZF5B7gWcAPPKKq+0Xkbmf7Q8DTwE3AIaAH+JSz+1zgKRGJ5vkxVf2xs+1+4EkRuRM4Dtzi0SmNWzZdrzEjm15WSFGBj86+ID2BIGVFXrYIjE+evkOq+jSRgBG77qGY5wp8JsF+h4GVwxzzNPD+7OZ0YnN7jAlE5piI/DuNGX9EhLlTizlxppeWjn4Wz7ZgMhobAZ+HrM3EmNFZ9+DUWDDJQxZMjBmdDVxMjQWTPNMfDHG2Z4ACnzBrUpFrr3Pfffe5dmxjYrl1rVnJJDUWTPJMa0f0NirF+HyJhvVkh42AN15x61qzgYupsWCSZ9yeFCvKjb7/xiTi1rVmY01SY8Ekz3jRkwvc6ftvTCJujWmaa2NNUmLBJM+4PSmWMV6rr6935bjzrGSSEgsmecarSbFWr149eiJjsmDdunWuHDdaMmnt6LcxU0mwYJJnvBr97tavRWO8UlrkZ2pJAYFQmLM9A7nOzphnwSTPNLd7U821YcOG0RMZM8bNte7BSbNgkmea23sBqJzubjDZtm2bq8c3xgvzpjntJjZwcVQWTPJIOKyDv7Aqp5fmODfGjH2D3YOtZDIqCyZ55FRXPwMhZeakIkoK/bnOjjFjXoVTMmmyYDIqCyZ5pPGcN1VcAI2Nja6/hjHg7rUWLcE3OZ8dMzwLJnkk2vheOc39Ki7rzWW84ua1VuUEk8azFkxGY8EkjzQNlkzcDyZu9f03Jt7WrVtdO/ZgyaTdgsloLJjkkaZz0cZ3G/1uJo4dO3a4duxoyaT5XB/hsA1cHImnwUREbhCRgyJySEQ2JdguIvJ1Z/teEVntrJ8vIr8QkQMisl9E/ixmn80i0igie5y/m7w8p/EkWjKp8KCayxivrF271rVjlxb5mTmpiEAozKkuu3vwSDwLJiLiBx4AbgSWA7eJyPK4ZDcCS52/DcCDzvog8BeqeilwNfCZuH2/pqqrnL8h0wKb886PMXE/mDz88MOuv4YxADt37nT1+NHSSYM1wo/Iy5LJlcAhVT2sqgHgCWB9XJr1wHc04kVguohUqGqzqu4GUNVO4ABQ5WHeJ4RGD6u5bAS8mSisET45XgaTKuBEzHIDFwaEUdOIyCLgcuClmNX3ONVij4jIjEQvLiIbRKROROra2trSPIXxqz8Y4lRXP36fMGeK+8FExL2Jt4zxknUPTo6XwSTRt0t8i9aIaURkMvB94F5V7XBWPwhcBKwCmoGvJnpxVd2qqrWqWlteXp5i1se/kzH35PK7OMOiMRNN1QynZGLBZEReBpMGYH7McjXQlGwaESkkEkj+XVV/EE2gqi2qGlLVMLCNSHWaidNw1rsBi8ZMJFXOZ8ZKJiPzMpjsApaKyGIRKQJuBbbHpdkOfMLp1XU10K6qzRKpM/kmcEBV/yl2BxGpiFn8ILDPvVMYv46d7gFgwcxJnrzemjVrPHkdY9y+1qqmlwHnf5CZxAq8eiFVDYrIPcCzgB94RFX3i8jdzvaHgKeBm4BDQA/wKWf3dwK3A6+KyB5n3RednltfEpFVRKrDjgIbPTmhcebYmW4AFs4q8+T13Oz7b0wst681q+ZKjmfBBMD58n86bt1DMc8V+EyC/X5N4vYUVPX2LGdzQjrulEy8CiZr1661gGI84fa1NqOskMnFBXT2BTnbHWDGpCLXXms8sxHweeJ8NZc3wcTtvv/GRLndDV1EWDQ78rk5crrb1dcazyyY5AFV5fiZaMnEmzYTY7xSU1Pj+msscj43R09ZMBmOBZM8cKY7QFd/kCnFBcwoK8x1dozJqqoq98cvL55twWQ0FkzywDGnVLJgVplngwkjzV/GTAzRkskRp7rYXMiCSR7wuvEd3L0tuDFeW2Qlk1FZMMkDh9u6gPO/rrywcaP10DYTR7Sa68ipbit1D8OCSR54oyUSTC6ZNyXHOTFmfJpRVsjUkgK6+oOc6grkOjtjkgWTPPBGSycAS+dYMDEmHSLCkvLJABxq7cpxbsYmCybj1EAozBefepXrvvI82351eNh0fQMhjp7uxiewpNy7aq7t2+PvlGOMO7y61i6tiPwYO9DcMWyaN1o6+fCDv2X9A79hX2O7J/kaKyyYjFP/+MzrPPbScY6c6uZ/PX2AJ+tOJEx3uK2bsEbaS0oK/Z7lz4u+/8aAd9fasnlTgeGDSXvvAP/lGy9Rf+wsvz9xjk9+axfnevKnSsyCyTjUcLaHR393FBG47crITZb/8ZnX6Q2ELkj7+snIhb907mRP8+hF339jAGpraz15nUsrnGByMnEweeiXb9Ha2c9lVdNYUT2NU139fOOFI57kbSywYDIOPfbScQZCytoVlfzDBy9jZfU0TncHeOqVxgvS7jlxDoCV86d7m0ljPNLUFD+ThTuWOdVcb7R0EQyFh2zrGwjxxMvHAdi87u3ctzYyq/hjLx8nEByadqKyYDLOBENhvr+7AYCPX7UAEeFP37UYgG//9sgF3RZfOX4OgMvnJ5yA0phxb/PmzZ68ztSSQhbOKiMQDPNaXFXX0682c7ZngHdUTWX1gumsXjCDZfOmcKY7wHMHWjzJX65ZMBlnXjh0ipaOfhbNKuPKxTMBuPEdFcyZUswbLV28ePjMYNreQIgDzR34BFZUT/M0n3fddZenr2fy15YtWzx7rSsXRT5zL8V8zgD+7cVjANx+9UJEBBHhltpIFfR36xs8y18uWTAZo06c6WHn3iY6+gaGrP/B7khV1kdqqgdvjVJU4OPWKxcA8P+cixoiVVzBsHLx3ClMKvZ0tgEbAW8mpKuXzALgpSOnB9fta2znlePnmFJSwLqV59sK16+qpMAn/PKNNto6+wfXqyq/e+s0L7zZNqEGQFowGYP2N7Vz/dd+yT2PvcLaf/k1LR2R+ds7+wb4yf6TANx8+dAG7o9duQC/T3h2/8nB9L842ArAu94228PcR1hvLjMRXbXEKZkcOcOA027y7d8eBeCWmvmUFp3vMTl7cjHvvaScUFjZ/vvz7Tr//Nyb3LbtRW7/5sv89Q/3TZiA4mkwEZEbROSgiBwSkU0JtouIfN3ZvldEVo+2r4jMFJGfisibzuO4ahzY19jOX/1gL0/WnUBV6Q+G+Isnf0/fQORCPXa6hy98fy+qyo7fN9MfDHP1kplUzxh6n61500r4wPK5BMPK4y8fJxxWfrwvEnjef+lcz89r9+7dnr+mMW6rnlHG0jmT6ewL8sKbkRLH9j1NiMAd1y68IP2HVlcD8F3n8/3i4dP883NvDm5/7KXj/HBPpLbh4MlOvvjUq3zrN0cIh8dfgPGs7kNE/MADwPVAA7BLRLar6msxyW4Eljp/VwEPAleNsu8m4DlVvd8JMpuAL7hxDj2BIANBZWppgbMcosAvDISUgWAYn0842x1ABMIKexvOUbtoJj97rYXT3QE+srqaf/n5mzz/RhtrVlRw/fK5bPhOPV39QR5/+QS7j51lamkhr5/sZNGsMr71qStZ/39/zfMH23j4V4cHe4vc5lRpxbv96oU8s+8kj/z6CLMmF3P8TA9V00u5YtG4iq/GjGk3X17Fl589yDdeOMKCmWUEQmGuXz434VxB71s2h9mTi3n9ZCf/sesE//zcm6jCPde9jfkzS/nC91/lb3+0n7KiAv77d39PR18QgNeaOrjxsnn8z50HKCv2s3nt22lq72N/Yzu31M4nFFZOdvRxxaIZnDjTy/yZpQSCYUqL/HT1BZleVkRYlUK/d+UF8aqIJSLXAJtV9Y+c5b8CUNX/HZPmYeB5VX3cWT4IvBdYNNy+0TSq2iwiFc7+l4yUl9raWq2rq0v5HH60p5E/e2IPpYV++oIhsvXWrayexsGWzsHSiE/gyY3XULtoJk+/2sx//ffzv/KXlE/iJ/e+h4IEF4mqcse3dvGrN9oG1/3tmuWDvb28VFlZ6VmXTZPfvL7WznYH+IMv/2Lwi7/AJzzzZ+9m6dzEtyv65q+P8Pc7z/9mvnzBdJ7ceA0FPuGu79Tzs5jeXiurp/FGSxe9AxeOGUtEhGG/hwr9wmSnrTQQDKNAod/HwlllPLnxmrQGMYtIvaomHNjjZTVXFRA7TLvBWZdMmpH2nauqzQDO45xELy4iG0SkTkTq2traEiUZVWtHP0UFPnoHIoGkwBdpAPdJ5E8EphQXUFLoo8AnrJw/Hb9PuGTuFN75tkjD3bJ5U/jaR1eyyLkd/C011Xz/09fyyB1XMLWkgEK/8D9vvoxap9fITZdV8PkPXEyBT5g9uZgHPrY6YSBxzpGv3rKSt1dGBlf98YoK7rh2UVrnmikLJMYrXl9rMyYV8bWPrmJKSQHFBT6++icrhw0kAJ+6dhHrVlYCsLxiKg/fXkOh34eIcP+HL+Mi5zZHH15dzXfvvpZvfrKWqSUFFPl9fO79S/nktYsQgXlTS7h++Vz8PmFqSQEXz52MKkwrHTrhXakTJAZCytmeAc72DNAdCNETCNHeO8DproA7d8NQVU/+gFuAb8Qs3w78S1ya/wTeFbP8HFAz0r7AubhjnB0tLzU1NZquvoGgnu3u13M9AQ2FwhoMhbVvIKh9A0HtDQQ1EAxp30BQu/oGVFU1GAoP7nuuO6DhcGQ5FArrue7AkGP3Bs7vF6+7f0ADwVBSeQyFwnqmqz+d08ua++67L6evb/JHrq613kBQu/sTf14TOdvdr6GY74OoUCis7b1Dvwtiv0NUVdt7AzrgfP47+wa0fyCk4XB4ME30+6jDOU7fQFDP9QS0rbNPWzp69VxPQNt7A9ra0acHT3akfK5RQJ0O873qZX/RBmB+zHI1EP+TYrg0RSPs2yIiFXq+mqs1q7mOU1zgp7hgaFT3+y6M8tGeuH7f+ZkNp8VMmevzyZBlYMRfC2VFyf+rfD5hxqSipNO7YcuWLZ4NJjMmF1L9dT+9LPFn0ucTppYM/S6IfM+cX47dPjlmQ7TLf7R0MsVJl+h7KpJpKJ9SnFK+k+VlNdcuYKmILBaRIuBWIP52n9uBTzi9uq4G2jVSdTXSvtuBO5zndwA/cvtEjDFjh/1oGRs8K5moalBE7gGeBfzAI6q6X0TudrY/BDwN3AQcAnqAT420r3Po+4EnReRO4DiRKjFjTJ6wzh5jg2e9ucaSdHtzmeTV19fbwEXjCRGZMAP/xrqx0pvLGGPMBGXBxLjCqzkmjDFjgwUTY4wxGbNgYowxJmN52QAvIm3AsVETjj+zgVO5zsQ4Yu9Xauz9Ss1EfL8Wqmp5og15GUwmKhGpG66nhbmQvV+psfcrNfn2flk1lzHGmIxZMDHGGJMxCyYTi82Vmxp7v1Jj71dq8ur9sjYTY4wxGbOSiTHGmIxZMDHGGJMxCybGGGMyZsFkghKRJSLyTRH5Xq7zMh6IyKUi8pCIfE9EPp3r/Ix1IvJeEXnBec/em+v8jAci8m7n/fqGiPw21/nJNgsm44iIPCIirSKyL279DSJyUEQOicgmAFU9rKp35ianY0OK79cBVb0b+BMgbwaaxUrl/QIU6AJKiMyQmpdSvMZecK6xncCjucivmyyYjC/fBm6IXSEifuAB4EZgOXCbiCz3Pmtj0rdJ4f0SkXXAr4HnvM3mmPFtkn+/XlDVG4EvAFs8zudY8m1S/0x+DHjcqwx6xYLJOKKqvwLOxK2+EjjklEQCwBPAes8zNwal+n6p6nZVvRb4uLc5HRtSeb9UNexsPwu4M6n4OJDqNSYiC4hMR97hbU7dZ8Fk/KsCTsQsNwBVIjJLRB4CLheRv8pN1sak4d6v94rI10XkYSLTR5uI4d6vDznv1b8B/zcnORu7Er5nzvM7gW95niMPeDYHvHGNJFinqnoauNvrzIwDw71fzwPPe5uVcWG49+sHwA+8zsw4kfA9A1DV+zzOi2esZDL+NQDzY5argaYc5WU8sPcrNfZ+pS4v3zMLJuPfLmCpiCwWkSLgVmB7jvM0ltn7lRp7v1KXl++ZBZNxREQeB34HXCIiDSJyp6oGgXuAZ4EDwJOquj+X+Rwr7P1Kjb1fqbP37Dy70aMxxpiMWcnEGGNMxiyYGGOMyZgFE2OMMRmzYGKMMSZjFkyMMcZkzIKJMcaYjFkwMSaGiIREZE/M36bR93KfRPxcRKY6y3NF5DEROSwi9SLyOxH54CjHOCIil8St+z8i8pcicpmIfNvFUzATnN2by5ihelV1VTYPKCIFzkC2TNwE/F5VO0REgB8Cj6rqx5zXWAisG+UYTxAZjb3F2ccHfAR4p6oeE5FqEVmgqsczzKvJQ1YyMSYJInJURLaIyG4ReVVEljnrJzkTJO0SkVdEJHqr8U+KyHdFZAfwExEpE5EnRWSviPyHiLwkIrUicqeIfC3mde4SkX9KkIWPAz9ynr8PCKjqQ9GNqnpMVf/FOYZfRL7s5GmviGx0kj1OJJhEvQc4qqrHnOUdcduNSZoFE2OGKo2r5vpozLZTqroaeBD4vLPur4Gfq+oVwHXAl0VkkrPtGuAOVX0f8F+Bs6q6Avh7oMZJ8wSwTkQKneVPkfgW5e8E6p3nbwd2j3AOdxKZM+MK4ArgLhFZrKp7gbCIrHTS3crQSZrqgHePcFxjhmXVXMYMNVI1V/SW6/XAh5znHyASDKLBpQRY4Dz/qapGJ056F/DPAKq6T0T2Os+7ReTnwBoROQAUquqrCV57pqp2JsqUiDzgHD/gBJAPACtE5CNOkmnAUuAITulERPYTmbDpb2MO1QpUDnPuxozIgokxyet3HkOc/+wI8GFVPRibUESuArpjV41w3G8AXwReZ/iJk4Ii4nNmONwPfDi6QVU/IyKziZQsoq/1WVV9NsFxHgd+AvwS2KuqrTHbSoDeEfJpzLCsmsuYzDwLfNZpFEdELh8m3a+BP3HSLAcui25Q1ZeIzH8x0tzgB4ElzvOfAyUi8umY7WVxefp0tOpMRC6OVr2p6lvAaeD+BK91MbBv2DM1ZgQWTIwZKr7N5P5R0v89UAjsFZF9znIi/wqUO9VbXwD2Au0x258EfqOqZ4fZ/z+B90JkmkPgZuAPnO6+LwOPOseFSEnnNWC3k6eHGVoL8TiwDHgq7jWuc17HmJTZLeiN8YCI+Im0h/SJyEXAc8DFqhpwtu8Evqaqzw2zfwXwHVW93qX8FROp+npXFroxmzxkbSbGeKMM+IVT9STAp1U1ICLTgZeJjCFJGEgAVLVZRLaJyFRV7XAhfwuATRZITLqsZGKMMSZj1mZijDEmYxZMjDHGZMyCiTHGmIxZMDHGGJMxCybGGGMyZsHEGGNMxv4/XkoT0JrOdrsAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.semilogx(EGeV, pa[0], lw=2)\n", "\n", "# the critical energy, onset of strong mixing regime\n", "plt.axvline(EminGeV(m_neV=ml.alp.m, g11=ml.alp.g, BmuG=ml.modules[0].B, n_cm3=ml.modules[0].nel),\n", " lw=1., ls='--', color='k')\n", "\n", "# maximum energy, end of strong mixing regime\n", "plt.axvline(EmaxGeV(g11=ml.alp.g, BmuG=ml.modules[0].B),\n", " lw=1., ls='-.', color='k')\n", "\n", "plt.xlabel(\"Energy (GeV)\")\n", "plt.ylabel(\"$P_{a\\gamma}$\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Sanity check\n", "\n", "We would expect to get the same result as above if we split up the calculations over multiple cells or domains that have the same magnetic field strength and orientation. Let's see if that's the case. For this purpose, we define a second module list where we set `L0` to some smaller value that fits $n$ times into our total legth. " ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "\u001b[0;36menvirons.py:\u001b[0;35m 257\u001b[0;0m --- \u001b[1;36m\u001b[1;36mINFO\u001b[1;0m\u001b[1;0m: Using inputted chi\n", "\u001b[0;36menvirons.py:\u001b[0;35m 270\u001b[0;0m --- \u001b[1;31m\u001b[1;31mWARNING\u001b[1;0m\u001b[1;0m: r_abell <= L0: assuming one domain from 0. to L0\n" ] } ], "source": [ "ml_many_cells = ModuleList(alp, src, pin=pin, EGeV=EGeV, seed=0)\n", "\n", "ml_many_cells.add_propagation(environ='ICMCell',\n", " order=0, # order of the module\n", " B0=1., # B field strength\n", " L0=0.001, # cell size\n", " nsim=1, # one single realization\n", " n0=1e-3, # electron density\n", " r_abell=10.001, # full path, chosen that we only have a single cell\n", " beta=0., \n", " eta=0.\n", " )\n", "\n", "ml_many_cells.modules[0].psin = np.ones_like(ml_many_cells.modules[0].psin) * np.pi / 2." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's check how many cells we have:" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(1,)\n", "(1, 10000)\n" ] } ], "source": [ "print(ml.modules[0].psin.shape) # first case\n", "print(ml_many_cells.modules[0].psin.shape) # many cell case" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we calculate the oscillation probability and compare:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "\u001b[0;36m core.py:\u001b[0;35m 652\u001b[0;0m --- \u001b[1;36m\u001b[1;36mINFO\u001b[1;0m\u001b[1;0m: Running Module 0: \n", "/Users/manuelmeyer/Python/gammaALPs/gammaALPs/base/transfer.py:799: UserWarning: Not all values of linear polarization are real values!\n", " warnings.warn(\"Not all values of linear polarization are real values!\")\n", "/Users/manuelmeyer/Python/gammaALPs/gammaALPs/base/transfer.py:802: UserWarning: Not all values of circular polarization are real values!\n", " warnings.warn(\"Not all values of circular polarization are real values!\")\n" ] } ], "source": [ "px_many_cells, py_many_cells, pa_many_cells = ml_many_cells.run()" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Absolute maximum difference: 6.081801728896608e-13\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEKCAYAAADXdbjqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAA41klEQVR4nO3deXxcZb348c93ZrLvadIkbbrTFQqlhFIUkEWwBS+VyyLIZbsoohTh/vRq9fdTweXKBZcrWqmALKKCuCCVnQsoa6EtlC6Utmlpm7Rp9n2bTOb7+2NOyjRM2qxzMpnv+/UKM+ec5znzPYd0vjnPeZ7ziKpijDHGDIXH7QCMMcbEPksmxhhjhsySiTHGmCGzZGKMMWbILJkYY4wZMksmxhhjhszndgBuyMvL06lTp7odhjHGxJT169fXqGp+pG1xmUymTp3KunXr3A7DGGNiiojs6WtbVJu5RGSJiGwTkVIRWRFh++UistH5eV1EjjtSXRHJFZHnRWSH85oTreMxxhgTErVkIiJeYCWwFJgHXCYi83oV+wD4hKoeC3wfuLsfdVcAL6jqTOAFZ9kYY0wURfPKZBFQqqq7VNUPPAIsCy+gqq+rar2zuAYo7kfdZcCDzvsHgc+M3CEYY4yJJJrJZCJQFrZc7qzry7XA0/2oW6CqFQDO6/hIOxOR60RknYisq66uHkT4xhhj+hLNZCIR1kV8yqSInEEomXxjoHX7oqp3q2qJqpbk50fsjGCMMWaQoplMyoFJYcvFwP7ehUTkWOBeYJmq1vajbqWIFDl1i4CqYY7bGGPMEUSza/BaYKaITAP2AZcCnwsvICKTgb8CV6jq9n7WXQ1cBdzmvD4+kgdhRhcNBmlrbaK5oRZ/Rwtt3mxavRm0+4NQv5uU6g0EA340GECD3dDtvAYDbJp4Kd0I3UFl9r7HSOs8gKiiPRe9zvQMlWmz2THuDFQhzV/DwgN/DLsuDgIgTtl1BRfRlFSEqjK77h8Ut2wk0kV0S0IeawovR5xr7jPLVx7cBwKKIAiKUJp7Gvsz5iNAYetWZtb9AyFUBhHnsl1QEdZM+jyIF4BjKh8nw1/l7AtUPrzAr02fzZ680wBI8dcx98Dfwvbp+bApwJtA2YQlBNIKSfB6yGveSmbHfjy+xNBPQiK+hGSS0rNJycwjLa+YtEQfHk+kxgQzlkUtmahqQESWA88CXuA+Vd0iItc721cB3wHGAb+S0C9+wGmailjX2fVtwKMici2wF7g4WsdkRl5TQy17du3g/eBEyuvbOVDbwJW7vkZqoJG0YBPZ2kSaBEhzyv+/rmv4XffZAFzifYnbE+7pc9//tnE+AeefwGOJj3C8pzRiuUcDn2BlYDIAc2UPNyT9ts99/veeWbyj7QB8z/cii33PRyy3NTiZe0s/dnD5G0l/wCfBiGWf2uvhoe5UAC72vs6VCQ/0+flX7Dit38d0R6Dg4DE9nbSy72PanMk7OtM5pvu5so9j2hGcyGL/HYhAeqKP1Z6v0uHLoD0pH39qIWRNJKVoLvnTj6Nw0lF4vN4+P9PEHonHybFKSkrUBi2OPhoMUla6kYp3niFh76sUtm1jglZRGpzAJ/0/PlhuY9LnyZS2g8vtmkizpNMpyfwx5bO8ln42KQleju3ezDktqwl6ElDxoR4viBf1+EA8/O+km/D4EvB4hOOrHyfDXw3itPyG/RVfkzbr4F/xqV11zD3Qc/HrlHHqCLCj8FzakkJ9QCbWv8m4lu0HrxwO/rmv0JGQzfaiT4cWFY7d+1tEg2HXMHrwyqgsZzHVGXNRVfJatjGl9lUOXu1oz39Cy2uKrwXxoMD8A4+R0VnpbPlwf6BUps5m+7gzUUJXWydU/NHZFnalpUHQbl7Ju4QaTwFd3UFOrH2c2a1r8QS78GgX3mAXvmAnycE29mgB13f/J23+bnwEKE2+ss//1z8IXs27Ez7LiVNzOXmilxNnFJGcmt5neTM6iMh6VS2JuM2SiXHbzuoW3nv+QRbu+DkTtfKQbZ2awI6EWdx71EqKc1KZlJvC7K6tZKRnkJ5bQFZugX0JjTKB7iCtHV00VWynqXIv7bVlBBr3I41lpDXvorBzNzf5v8RrwfkA3OD9Gzf4Huf9jMUw/2KOPfOz+BISXT4KE4klk14smbhPg0HWllaw8tV9/HN7NWd51vObxJ9QRyYfZJTQPfU0CuadxsSj5tsXyxhU39LJ22UNvLW7jvkbvs+nO588uO0Aeeye/jmOu/DrpKRluBil6c2SSS+WTNxVVrqJuj//Bx+0JnBz13JSErwsOyaPf5u4n7mLz8Xri8tHxsW1ij3b2PPKw0zc+UcmaaijZhW57Fz4TRb/yxcQsRv6o4Elk14smbhDg0HefOSHLNz2PyRKgAZN56GTn+TyU+eRm2ZXHwaC3d1sevmvpL36I47q3sly/400zvgX7rjoOAqzkt0OL+5ZMunFkkn0tTTVs2PV5zi+7XUA1mYvZdqld5BXOOkINU08CnZ38+oLf+MrazJoaOsiJzWBey+dzQmzprgdWlw7XDKxybHMiKur2kfFnWdzfNvrNJHGOx9byYk3P2KJxPTJ4/Vy2jkX8tzNp3HqzDyK27cx4/cns+7Jvrt6G3dZMjEjqq7Vz/p7vszMwA72SQHNVz7P8ef8m9thmRgxPjOZ+68+ka9O2Um2tHL8W//J+ifvdTssE4ElEzNi2vwB/v2BtdzUfCV/S/w0SV94nonTj3Y7LBNjfF4Pp1//P7wx+Yt4RTn2ra+z6Z9/dTss04slEzMiNBjk6396lw1lDeRk53Dy8t+QN8Hau83gLb76NtYUXk6CdDP5peXs27XlyJVM1FgyMSPirT//hDPe/w75iV08+O+LKMi0njhmaMTjYdEXfsGG1JPJopXO332Ozo62I1c0UWHJxAy7vds3sGDLf3Oh91V+ubiRo8bbCHUzPDxeL9Ov+z3lUsTk7r08/vif3A7JOCyZmGGlwSDNf76RJOlibdYSTjr3KrdDMmNMZvY4ms67iwu7bmXFhjzeLWtwOySDJRMzzNY+vpKj/RupI5OZV/zc7XDMGDWv5AxO/NgnCSrc8vctxON4udHGkokZNh1tLUx996cA7Dz+m2TnFbockRnLbj57FnlpiWSXv8SaZx92O5y4Z8nEDJsNf7md8dRR6p3BCZ/+otvhmDEuPcnHTxcc4P7EO5i25jv4OzvcDimuWTIxw6K5o4u9O98DoO3U/2sTH5mo+PjSy9jtmUQh1bz75K/dDieuRTWZiMgSEdkmIqUisiLC9jki8oaIdIrI18LWzxaRDWE/TSJys7PtFhHZF7bt3CgeknH8cW0ZX++4hpvz7mH+aRe4HY6JE16fj5oFywEo3LSK7kDA5YjiV9SSiYh4gZXAUmAecJmIzOtVrA74CvDj8JWquk1VF6jqAuAEoA14LKzIz3q2q+pTI3UMJrJAd5D7X9sNwKfP/ATisQteEz0Llv47+6WASbqfDc896HY4cSua/+oXAaWquktV/cAjwLLwAqpapaprga7D7OcsYKeq7hm5UM1AvPXPJ5jUtJ7p41I5c854t8MxccaXkEjZnGsBSNlwv8vRxK9oJpOJQFnYcrmzbqAuBXp33VguIhtF5D4RyYlUSUSuE5F1IrKuurp6EB9r+pK95nYeSfwBt854H4/HJjEy0TdvyRdo0yTm+Tex5/233Q4nLkUzmUT6lhlQ53ARSQTOB8KHvd4FzAAWABXATyLVVdW7VbVEVUvy8/MH8rHmMPbt2sI8/ybaNInjP3mp2+GYOJWRlcsbBZfx466LeXSr9epyQzSTSTkQPoFFMbB/gPtYCrytqpU9K1S1UlW7VTUI3EOoOc1Eyd4XQ48D35J9OumZES8KjYmK8ed/j192X8DvNrbQGeh2O5y4E81kshaYKSLTnCuMS4HVA9zHZfRq4hKRorDFC4DNQ4rS9Ft3IMC08tD/wtRF9tgU4675xVnMLcqksb2Ll7fXuB1O3PFF64NUNSAiy4FnAS9wn6puEZHrne2rRKQQWAdkAkGn++88VW0SkVTgbKD3aLjbRWQBoSaz3RG2mxHy/lvPcTQ17JfxzF28xO1wjOHiOUnUVv2N4AtPw7xfuh1OXIlaMgFwuu0+1WvdqrD3Bwg1f0Wq2waMi7D+imEO0/RTyzt/AWBP4RIm2CBFMwp8amYGE994lNaaZNpbf0RKWobbIcUNGxBgBkVVeacukX06jnEnXuh2OMYAMHH6XLb7ZpEmHWx9+c9uhxNXLJmYQdlY3shtrZ/moqS7Oeq409wOx5iD6qaGHoIR3PqEy5HEF0smZlCe3nwAgE8dU4THa79GZvSYsCj0OJ+jmtbY41WiyL4FzKB0bnqcfOo5Z16B26EYc4hJRx1LuRSSTQs73n7J7XDihiUTM2DV+3fz3bb/4p9J/4cTJtmUvGZ0EY+HfXmnAFD/7pMuRxM/LJmYAdv9VqgtekfqApKSkl2OxpiPSpr/Gf7afQpPNU51O5S4EdWuwWZskF2hpoP2yZ9wORJjIptz8rlc+ryPjpogNzZ3MD7D/ugZaXZlYgYk2N3N9Ka3AChaeJ7L0RgTWXKClxOn5gLw1gd1LkcTHyyZmAH5YMub5NLEAfKYPPNYt8Mxpk+nTkpkqedNWtY/6nYoccGSiRmQ6i2hJq6yrIU2CZYZ1U7LquauxJ+zuOxet0OJC/ZtYAaksaqcgHrQSSe7HYoxhzV9wWm0ayJTg2XUVpa7Hc6YZ8nE9Juq8u2WC5nfeS/jTr7c7XCMOazEpGR2JodmBt/z9vMuRzP2WTIx/VZW105VcyfJqRlMn2DT85rRr7lwMQBdO19xOZKxz5KJ6bd3dpYhBDlhSi4iNj2vGf0yZn4cgJyGTS5HMvZZMjH9VrjmB7yddD2XpNoc2yY2TJl/CkEVpnbtpKO91e1wxjRLJqbfxjdsJEdamDx5qtuhGNMvGVm5fOCdwg4tpnTXTrfDGdOimkxEZImIbBORUhFZEWH7HBF5Q0Q6ReRrvbbtFpFNIrJBRNaFrc8VkedFZIfzahORj4D21mYmd++hW4Upx1hPLhM7fj3nfs7z/4i1DTZR1kiKWjIRES+wElgKzAMuE5F5vYrVAV8BftzHbs5Q1QWqWhK2bgXwgqrOBF5wls0w2/Pem3hF2eudYrPXmZhy3JTQBK0byhrcDWSMi+aVySKgVFV3qaofeARYFl5AVatUdS3QNYD9LgMedN4/CHxmGGI1vTSUhh6hUp3ZO/8bM7otmJQNKBV7trsdypgWzWQyESgLWy531vWXAs+JyHoRuS5sfYGqVgA4r9ZndQR4D2wAQIsWuBqHMQM1e3w6ryXdxKPt19FQc8DtcMasaCaTSH1JdQD1P66qCwk1k90gIgOaK1ZErhORdSKyrrq6eiBVDZDfvBWA7BknuhyJMQPj83lpTsgHoHzrmy5HM3ZFM5mUA5PClouB/f2trKr7ndcq4DFCzWYAlSJSBOC8VvVR/25VLVHVkvz8/EGEH7/a/AH+o+Pz3Bq4minzFh25gjGjTGPWbABa9m5wN5AxLJrJZC0wU0SmiUgicCmwuj8VRSRNRDJ63gPnAJudzauBq5z3VwGPD2vUhvf2N/FO8CjezL+I5FSbWdHEoIJjAPBWbT5CQTNYUZscS1UDIrIceBbwAvep6hYRud7ZvkpECoF1QCYQFJGbCfX8ygMec0Zd+4A/qOozzq5vAx4VkWuBvcDF0TqmeLG1ogmAeRMyXY7EmMHJmnY8vAfjmu0m/EiJ6kyLqvoU8FSvdavC3h8g1PzVWxNwXB/7rAXOGsYwTS+ZWx7iP33bGZ9xtduhGDMok+aUEHxCKO4uw9/ZQaJNNz3sbAS8OaLZlU9xg281c1Ib3Q7FmEFJy8hmn6eIROmmbPsGt8MZk2wOeHNYGgwyoWs3AEUzT3A3GGOG4G9FN/HSB+1c0Z7PDLeDGYPsysQcVuW+XWTSRj0ZjCuYdOQKxoxWM87ibZ3F5qpOtyMZkyyZmMOq3BF6QvD+xGk2Ta+JaXOKQh1Itlc2uxzJ2GTNXOaw2spD80C0ZM1yORJjhmZmjocVvocp3tcC/M3tcMYcSybmsBJq3gNACo52ORJjhqY4L4trvM+QFOyiubGOjKxct0MaU6zdwhzWLn827wWnkDF1oduhGDMkvoRE9ntDjwPcv9NmXhxulkxMn4JB5TutF3Gu/0dMOPrjbodjzJDVpU4FoGmvJZPhZsnE9OlAUwcdXUHy0hPJSklwOxxjhsyfE7r3F6jc5nIkY48lE9OnveX7KKSW6ePS3A7FmGGRWDQHgOTGHS5HMvZYMjF9ki1/Zk3yjfyfrlVHLmxMDMidEnrgY177HpcjGXssmZi+1Th/veVMdTUMY4bLhBnz2RScxluBGXT4A26HM6ZY12DTp5TmDwBILpzjciTGDI+k5FRuzPgZu2vbOLqujTmF9iTs4WJXJqZP+R17AcidYvO+m7FjWl7oHuDumjaXIxlbLJmYiDraWijQGrrUS9HUuW6HY8ywmZqbwkSqqdlX6nYoY4olExNRxQdb8IhywFNAQmKS2+EYM2w+1f4kryXfxMxtv3Y7lDElqslERJaIyDYRKRWRFRG2zxGRN0SkU0S+FrZ+koi8JCJbRWSLiNwUtu0WEdknIhucn3OjdTxjWf3e0GNUalOmuByJMcMrpXAmAOkt1qNrOEXtBryIeIGVwNlAObBWRFar6nthxeqArwCf6VU9AHxVVd925oJfLyLPh9X9mar+eGSPIL6s9R7PDzu/y3lzJ7PA7WCMGUZ5k0PNtuP8+1yOZGyJ5pXJIqBUVXepqh94BFgWXkBVq1R1LdDVa32Fqr7tvG8GtgIToxN2fNrWAOt1NinTFrkdijHDqmDSUXSpl0Jq6GhrcTucMSOayWQiUBa2XM4gEoKITAWOB94MW71cRDaKyH0ikjOkKA0AZXWhni5TclNdjsSY4eVLSOSApwCAit3vuxzN2BHNZCIR1umAdiCSDvwFuFlVm5zVdwEzgAVABfCTPupeJyLrRGRddXX1QD42Ll1RdQff8v2eyRluR2LM8KtLLgagodySyXCJ5qDFciB83tdiYH9/K4tIAqFE8ntV/WvPelWtDCtzD/BEpPqqejdwN0BJScmAkli86Whr4fzgi3R5vUiuDeoyY097+hRof4vOKntG13CJ5pXJWmCmiEwTkUTgUmB1fyqKiAC/Abaq6k97bSsKW7wA2DxM8catyrJQ//sqTx6+hESXozFm+O2feTmf6fwezyQtcTuUMSNqVyaqGhCR5cCzgBe4T1W3iMj1zvZVIlIIrAMygaCI3AzMA44FrgA2icgGZ5ffUtWngNtFZAGhJrPdwBejdUxjVcO+7UwB6hInWC8HMyblTDmaDdpGWmOk1nczGFF9Npfz5f9Ur3Wrwt4fINT81durRL7ngqpeMZwxGuio2glAW9qkI5Q0JjZNygl1LCmvb3c5krHDHvRoPiJYHxrMFcya7HIkxoyM4uxkvu17iEnNNXQHTsXrs6/CobLHqZiPSG4OPeAxMX+6y5EYMzKSE30s873BOZ611BzY7XY4Y4IlE/MR27sL2RCcTsZEe/S8GbtqfaGxJrXl9sDH4WDXduYQqsr32i+h1X8hG2ae6HY4xoyYlpSJ0LyN1spdbocyJtiViTlEXaufVn83GUk+slIS3A7HmBHjzwj19QnU2QMfh4MlE3OIfZVVFFHLlJwkQsN7jBmbPDmhJ2J7G/e6HMnYYM1c5hD+95/hjeSvsb7rNOB0t8MxZsSk5E8DILWt3w/iMIdhVybmEF01oXnfu9JtuKIZ27KKZ7MuOIvNQesCPxzsysQcwtPkzPGQbQMWzdhWMGUun/Dfgi8gXNQdxOe1v62Hws6eOUSyc8mfNM5mWDRjW3KCl/EZSQSCyoGmDrfDiXmWTMwhMjurAMgomOpuIMZEwZTsBIqliorKyiMXNodlycQcYlwwlExyi6a5HIkxI+8bHXfyatLN8P6TbocS8yyZmIOaG+vIpI0OTSAnr+jIFYyJcV2ZzliT2t3uBjIG2A14c1BFK1zVeQtzsgL8l8f+zjBjnyerGPaBt3mf26HEPEsm5qB9zQHe1lmk5uW5HYoxUZGcF+q1mNxh90yGyv78NAftbwjN7VCUlexyJMZER8b40L3BLH+Vy5HEvkElExGZPdyBGPdllq7m//keYqHXnqJq4sO4CaFkktddjQaDLkcT2wZ7ZXKZiKwQkQE1k4nIEhHZJiKlIrIiwvY5IvKGiHSKyNf6U1dEckXkeRHZ4bzmDPKY4l5h1St83vc0R1HudijGREVmdh5tmkSadNDcVOd2ODFtUMlEVW8B/g7cJSJnSD+eCCgiXmAlsJTQvO6Xici8XsXqgK8APx5A3RXAC6o6E3jBWTaDkNpeAUByng1YNPFBPB5uSf0m53d+n4pWa/UfisE2c90K/AehednPBH7dj2qLgFJV3aWqfuARYFl4AVWtUtW1QNcA6i4DHnTePwh8ZuBHZACyu0I3IbMKp7obiDFRtG/cyWzUGexvCbgdSkwbbG+ul1T1HwOsMxEoC1suB04ahroFqloBoKoVIjI+0g5E5DrgOoDJk+3Bbr0Fu7vJD9aCQP7EGW6HY0zUFGWlAHCg0R6pMhSDva47XUSeFpF7ReSGftaJ1BSmUagbKqx6t6qWqGpJfn7+QKrGhbrqfSRKgHoySEnLcDscY6KmRN7nFt8DZO74m9uhxLTBJpNsYA3wQ6C/PbvKgfBH0RYD/Z1I4HB1K0WkCMB5tT5+g1C3PzR1aa034oWdMWPWNPZxte85xle/7nYoMW2wyaQO8BL64u5vF4i1wEwRmSYiicClwOphqLsauMp5fxXweD/3acJUN3fwdvAoKlJnuR2KMVGVPM4ZuNh+wOVIYtug7pmo6vdEZAJwJ7C5n3UCIrIceJZQIrpPVbeIyPXO9lUiUgisAzKBoIjcDMxT1aZIdZ1d3wY8KiLXAnuBiwdzTPFuq2cWP/B/j6tnTuVUt4MxJooynSdkZ3VZo8ZQHDGZiMhVwE8IXcU8Adygqs2quh+4diAfpqpPAU/1Wrcq7P0BQk1Y/arrrK8FzhpIHOaj9jeEbj7a6HcTb8ZNmA5AXncNGgwi9ly6QenPWfs2cDYwB9gD/NeIRmRc0VRfjZduCi2ZmDiTkZlDqyaTKp0019e4HU7M6k8yaVLVd5wxIN8mNObDjDHXlq1gR9KVHNWx5ciFjRlDxOOhxhvq4VlT8YHL0cSu/twzKXLGaGwF3gcSRjYk44asQA0eUbLyJrgdijFRtzd5DpUtqUhzG9PdDiZG9SeZfBc4FrgcmA+ki8hTwLvARlV9eATjM1GgwSB5wToQGFdoAzpN/Pn79G/z6LpyfuiZzoluBxOjjphMVPXu8GURKSaUXOYD5wKWTGJcQ20lORKgiVQy0zPdDseYqOsZBV/RYKPgB2vAXYNVtZzQIMKP9Kwysam+ci85QJ1nHJZKTDyakJ1MIl001FXS/3HYJpzNtGhort4LQEuCzbBo4tOctvVsT76aTbuPB/7hdjgxyTpUGzrrQvNftycXuByJMe7IGBfqeJLRZV2DB8uSiWFz8gl82f8Vthf/q9uhGOOKnIJQx5PcoE2QNViWTAyl/myeCi6mu7i/MwIYM7ZkjyvArz4yaaW9tdntcGKSJRNDVVOoB8v4DBv9buKTeDzUeHIBqDuwx+VoYpMlE8PJlX/gC94nmJjU7nYoxrimyRfqgNJYVXaEkiYS681lWNb2V/ISGqhM+arboRjjmrakPOiCtrpyt0OJSXZlEucC/k5ytZGgCrnjJx25gjFj1IYJl/JF/3+wPWm+26HEJEsmca62qhyPKPWSRUJiktvhGOOazgmLeTZ4Ih902tDdwbBkEud62ofrveNcjsQYdxVkhv6YqmzqdDmS2BTVZCIiS0Rkm4iUisiKCNtFRO50tm8UkYXO+tkisiHsp8mZhRERuUVE9oVtOzeaxxTrWqpDyaQ1Md/lSIxx16SEZr7kXc3CikfcDiUmRe0GvIh4gZWEJtoqB9aKyGpVfS+s2FJgpvNzEnAXcJKqbgMWhO1nH/BYWL2fqeqPR/wgxqCuhtDo945UG/1u4tv4xE6+kfAIZc0TgP92O5yYE80rk0VAqaruUlU/8AiwrFeZZcBvNWQNkC0iRb3KnAXsVFXrDD4Mmjq62a+5dKdPdDsUY1yVW9gzCr4WDQZdjib2RDOZTATCO3CXO+sGWuZSPvrY++VOs9h9IpIT6cNF5DoRWSci66qrqwce/Rj1TMqn+VjnLyk75ktuh2KMq9IzsmnTJNKkk+amerfDiTnRTCYSYZ0OpIyIJALnA38K234XMINQM1gF8JNIH66qd6tqiaqW5Ofb/YEeVc3O6PdMG/1u4pt4PNR6Qh1R6m0U/IBFM5mUA+EDGYqB/QMssxR4W1Ure1aoaqWqdqtqELgHm6N+QCob2gAotGRiDE0JoWTSXG2j4AcqmslkLTBTRKY5VxiXAqt7lVkNXOn06loMNKpqRdj2y+jVxNXrnsoFwObhD33s+mPzFbyZ9GUKE607pDEdyeMBaHemZTD9F7XeXKoaEJHlwLOAF7hPVbeIyPXO9lWEZm88FygF2oBreuqLSCqhnmBf7LXr20VkAaHmsN0Rtps+tLc2k0Mzfnwk5Ng4E2M60oopa8inqb3L7VBiTlSfzaWqT9Frul8nifS8V+CGPuq2AR/5xlPVK4Y5zLhRW7GHYqDGk8sEj41fNWbrvJu5/INPcXXqVM5yO5gYY98gcazJaRfueVqqMfGuwLl3WOlMy2D6z5JJHGurDT0dtS3JercZAx8mk+rGVpcjiT32CPo4FnBGv/tt9LsxAEzw1LIu6Xr8NcnADrfDiSmWTOJZ84HQa3qhu3EYM0rk5ReRLE34g61oMIjYvcR+s2QSx15L/DjPdfk4ddJpbodizKiQnJpOI2lkSSv1tZXk5Pd+mpPpi6XdOPZm11Hc172UpMnHux2KMaNGfc8o+Mq9LkcSWyyZxLEDTo8VG/1uzIeaE0K9G1tsFPyAWDNXnNJgkCXNf6HCk0VBxjluh2PMqNGRPB46oaPeRsEPhCWTONXUUMu3vA/R4kkhLfmHbodjzKgRSCuARuhu7P3oQHM4lkziVN2B3WQBtd5xpLsdjDGjSH3xWXx3D+QnLuJkt4OJIXbPJE71tAf3tA8bY0I8kxfxYPen2BCY7HYoMcWSSZzqqAuNfm93npJqjAn58JEq9iTtgbBmrjgVcNqDA6k2YNGYcAVpHi7wvMLEej9witvhxAxLJnHK44x+l0xLJsaEy89I5icJq9BuCHTdhi8h0e2QYoI1c8Wprs52AuohKXuC26EYM6okJCZRJ1l4Ramrsu7B/WXJJE7dkXwjszp/S3DOeW6HYsyo0+ANjYJvqLS54PsrqslERJaIyDYRKRWRFRG2i4jc6WzfKCILw7btFpFNIrJBRNaFrc8VkedFZIfzmhOt44lllU2dBPFQkG0dg43prTUxNC1DS3W5y5HEjqglExHxAiuBpcA84DIRmder2FJgpvNzHXBXr+1nqOoCVS0JW7cCeEFVZwIvOMvmMLqDSnVLqKfK+Ax7lIoxvXU40zJ0NVgy6a9oXpksAkpVdZeq+oFHgGW9yiwDfqsha4BsETnSYzuXAQ867x8EPjOMMY9J9ZVlvJbwZR5O/m8SfdbSaUxvwfTQ107QRsH3WzS/SSYC4U9OK3fW9beMAs+JyHoRuS6sTIGqVgA4rzZw4gjqK/dQKPWM9za5HYoxo5I3s4hO9eHvaHc7lJgRza7BEmGdDqDMx1V1v4iMB54XkfdV9eV+f3goAV0HMHlyfI9sbakJXbq3JNp0vcZE0j7vEmavPYpTUvP5hNvBxIhoXpmUA5PClouB3teQfZZR1Z7XKuAxQs1mAJU9TWHOa1WkD1fVu1W1RFVL8vPj+0vU74x+77TR78ZEND47HRAqnWkazJFFM5msBWaKyDQRSQQuBVb3KrMauNLp1bUYaFTVChFJE5EMABFJA84BNofVucp5fxXw+EgfSKzTpgoAgjZdrzER9czxc8CSSb9FrZlLVQMishx4FvAC96nqFhG53tm+CngKOBcoBdqAa5zqBcBjItIT8x9U9Rln223AoyJyLbAXuDhKhxSzPK09o99twKIxkWSnJvBY0neZpJW0tWwiNT3b7ZBGvag+TkVVnyKUMMLXrQp7r8ANEertAo7rY5+1wFnDG+nYltweaglMyu3d/8EYAyAiFHiayNMmyg6UkXpUttshjXrWLzQOPe05jbsD55EyofcwH2NMjyZfaBR8k80F3y/2oMc49EjHYuoDJ3BB8Wy3QzFm1GpNGg9d79FeZwMX+8OuTOJMZ6Cb+rYufB5hXJo9DdWYvnQ5o+ADDfawx/6wK5M4U11dxb94XqctbRIeT6RhPcYYADKLQgMNnOkazOFZMokzLXs38YvEX7JNZwPXux2OMaOWLzvUQSWxzZJJf1gyiTOttaFL9rak+B64acyR+CYcx68C59Pqnc/CIxePe3bPJM70PAXVn1LgciTGjG7Zk4/h9sClPO4vOXJhY8kk7vSMfs+w0e/GHE6BMwq+qqmT0BA4cziWTOKMr60y9Jplo9+NOZyURC+fSC7lXH2Z+oZ6t8MZ9eyeSZxJ7Qglk+TcSUcoaYy51XMPUxPL2Fm+jNycxW6HM6rZlUmcyeiqASC7aKq7gRgTA5oTnOl7a2wU/JHYlUkcCQaVT/p/THZ3Pf+cPMftcIwZ9TpSxkMndNbaKPgjsSuTOFLT0klnt9CVVkhyUpLb4Rgz6gXSQ/cWuxttFPyRWDKJI/saQlOQTshOdjkSY2KDNzt0b9HXbMnkSKyZK44E3n+GxxNvp1TPBE51OxxjRr2U/Kmh17YKdwOJAXZlEke6q7dznGcXxd4Gt0MxJiZkFk4HIN3puGL6Zlcm8cRp99WsYpcDMSY2jJ8yl5KOu2jyZvN+UO3hqIcR1SsTEVkiIttEpFREVkTYLiJyp7N9o4gsdNZPEpGXRGSriGwRkZvC6twiIvtEZIPzc240jymWJLbuByAhx8aYGNMfKclJBNPy8XcrNS2dboczqkUtmYiIF1gJLAXmAZeJSO+p/pYCM52f64C7nPUB4KuqOhdYDNzQq+7PVHWB83PItMDmQxmdoaefpo2f6m4gxsSQidkpAJQ7HVhMZNG8MlkElKrqLlX1A48Ay3qVWQb8VkPWANkiUqSqFar6NoCqNgNbAZvAfIByAtWh16JpLkdiTOz4d32MpxK/SeC9J90OZVSLZjKZCJSFLZfz0YRwxDIiMhU4HngzbPVyp1nsPhHJifThInKdiKwTkXXV1dWDPITY1dnRRh4NBNRDXuEUt8MxJmYU+ZqZ59lDsHq726GMatFMJpHuXPV+FOdhy4hIOvAX4GZVbXJW3wXMABYAFcBPIn24qt6tqiWqWpKfH39zeVTWN/PrwHms9p2D12f9LozpN2esCY02Cv5wovmtUg6E3/ktBvb3t4yIJBBKJL9X1b/2FFDVyp73InIP8MTwhj02lLV6+VHgck4szuFf3Q7GmBiSNG4yAMltvb+uTLhoXpmsBWaKyDQRSQQuBVb3KrMauNLp1bUYaFTVChER4DfAVlX9aXgFESkKW7wA2DxyhxC79tS2ATA5N83lSIyJLekFobEmmR02cPFwonZloqoBEVkOPAt4gftUdYuIXO9sXwU8BZwLlAJtwDVO9Y8DVwCbRGSDs+5bTs+t20VkAaHmsN3AF6NyQDGms2w9J8ouZto8JsYMSN7Eo0KvwSqXIxndotp47nz5P9Vr3aqw9wrcEKHeq0S+n4KqXjHMYY5JCz74DdckvcI6fzqh/gvGmP7IHldAi6aQKW001BwgO89mKY3E7sTGicyO0Oj3jKKZLkdiTGwRj4dnUs7lQEuAU+tayc5zO6LRyZJJHNBgkIJABUjo8RDGmIH5x+TlPLGxggktSRzndjCjlD3oMQ7U11SQLu00awrZ4wrcDseYmDMtL9RxZXdNq8uRjF6WTOJA1d73Aaj0TUA89r/cmIGakQWLPe+RsPcVt0MZtayZKw607N8BQFOKPYHGmMGYI2U8kvgDSitmAJ93O5xRyf5MjQOd1btCr5nTXY7EmNhUMC30XNnCwH40GHQ5mtHJkkkceCjhEhZ1rKTx2GuOXNgY8xHZ4wppIo10aae2yqbwjcSSSRzYXtVCFTlMnmxXJsYMhng8VPhCk8pV7tzocjSjkyWTGNXl7+TNX1xJ2a1zWfO7W/os19HVze7aVjwC0/PtUSrGDFZjZmiMVvOed/oss3vrOt7/4cncd8dX2byvMVqhjQqWTGLU+t/cxEm1jzNJ91Oy4+c8/7+R5wTbt20dTyd8ndvSHyU5wRvlKI0ZO3T8MQB4qrZE3N5YX0PaHy9isn8nd9Uez9X3r6WhzR/NEF1lySQGVezZxsIDjxJU4cXxV7HEfxsr3kig3d/9kbJ1O9cz21POzMQaFyI1ZuzImBp6DFFia+R7Ju/96VbyqWd3wnQmTZxITUsn977yQTRDdJUlkxi0+9mVJEo3b2eeyRlf+jlpE4+mttXPY+989Jc8WLYWgI7xC6IcpTFjS/Exp3Bix0ouaV9BoPvQHl0dbS3M3R+aGcPzqR/wf88/jjwakTUr8Xd2uBFu1FkyiTGBLj8z9v8dgNSPfQER4d9PmQYob7385Ee6LebWh24WZh71sWiHasyYkpmeTuq4ifgDynsVTYds2/z8g2TTQql3BrNPOJOFk3N4OO3HfFUfZPNLj7gUcXRZMokx6zesx6tdlEsRc0/6FABLjyni4ZQ7+J/WFbz3xtMHy7a3NjM18AHdKkw99uNuhWzMmLFoai4Ab5VWHrI+Y9ODANTNuxLxeBARamdcAID33T9EN0iXWDIZpcrq2nhi436aOroOWf+70mRO6lzJiyV3HXw0SqLPA8UlAHS+8euDZXdteJkE6WaPdwppGdlRi92YsWppzj6eS/xPSt666eC67VveZnZgG02kMn/JtQfXzzzrGrrUy9Fta6k5UHZwvQaDbHntSTa9/NiYGgBpyWQU2rK/kbN/9k++/Yd/su32s6jevxuA5o4unttygAA+zvzYSYfUmbFkOQH1ML/51Q/Lb3oSgKr8k6MZvjFj1tyZM5jl2cf0tnfp8ncCcPd7Xj7ZeTvPTP82KWkZB8uOKyhmS9oifBKk9MUHDq5f88A3OPr5zzH/xat5a+XVYyahRDWZiMgSEdkmIqUisiLCdhGRO53tG0Vk4ZHqikiuiDwvIjuc15xoHc9wKH33Nd66899467E70WCQzkA3X330XTq6gnzV9ydODL7L/t9+Hg0GefXl/yUl0Mji6bkU56Qesp/xE6exMf0UEqSb0mdWEgwq9zYs5FeB80lbeIlLR2fM2FI0ZTa7PZPIpI33Xn2M6uZOVm/Yz06KOem8qz5SPjj/swCML/0zGgzy3utPctKeew5uP6n2cdY/EWpN2HagmZd+tZw1D/+QYPdHe2aOdlF70KOIeIGVwNlAObBWRFar6nthxZYCM52fk4C7gJOOUHcF8IKq3uYkmRXAN0biGNr8AboCSmZyaLxGW2sTvoREuvydBPydiNdHfXcKIuCp/4DK91+n+Ngz2PP6nwm21jD5jM+zb/X3+GvzXFIWXMgF+fuZ/vQVHCXtUPd33tr7Jl3J41hS00xX7uc47ZLbaXrgNI7rWMua332X+R/8ideS6lg76/cR40s8+Qvwvy9z9N7f85dXLuV/G4rYmn0N15WcPhKnw5i4VDHlfKZ+sBLvm3fx6D7F353B2fMKmTLuo4OC551+CTVv3sL04G5efvIh7t/YwXd1PBXF5+HNncKiTd9l1tvf459Fi7nxySr+rauTr1fdzlu/2EzSscvIfeVWOj0pdH3qdnYkzWPLvkYuLpmEt2YrzVV7Oarkk1Tu3U7B5Fl0dXaQnJ5FS5eQnZpIUJUEb/SuFyQ0U24UPkjkZOAWVf2Us/xNAFX9UViZXwP/UNWHneVtwOnA1L7q9pRR1QoRKXLqzz5cLCUlJbpu3boBH8PjG/bx1KN387OEu0jGj0cOPXd1mk5J5yqCePis9yX+O+GeiPu5yf9lHg+ewgzZxx8Tv09dQiGTunaTIqEBTt0qbPuXx5hXcgZvP30/C9+8+WDdvZ6JTPjmBnwJiR/ZrwaDbLr9bI7tWMfdgfP4r8DlfOfT85zeXsaY4dBQcwDPLxeSSWhuk9eCxzD+S08ysyg7Yvk1f/g+Se//jS/6/4Mqcji5OJHffvET+Lw+Nvz4PI5ve50LO7/Lep3NkqIWflp3I6nSecg+NgRn8Bn/9wH4tOcNfpn4CwCCKod8D33FfwOrg6HONqf6tnCn704AErULQemURK4b9wAPXX/GoAYxi8h6VS2JtC2azVwTgbKw5XJnXX/KHK5ugapWADiv4yN9uIhcJyLrRGRddXX1oA6gqqmTRb4dpEonHlG6NPQ/o1uFbhVeCh5PWlIiyQkeVutpvJ1wAgH18IFnCpuTFgCwyzOVf1myhKnjUtmpE7ln1t1M/8Zr7Dz7NzSRhl+9rDvm28wrOQOAhUuv4Y2pX6JLvdSQTeBf74+YSCD0/KAJVz9AqXcGWbRy3rFFXPWxqYM6VmNMZNl5hew69ac0kUqHJpB04pV9JhKARZ/9Fg/N/hVV5DCvKJOfX3UqCQmJiMfD5Kvu4Q+JF7FeZ3PhwmLuvOFidp3T813g441Jn+fN/It4qPtsCjOTOXvueG7wPU4Taez2TMYjSiMfXhHN9ewlxUkSZ7KOHJrJoZk06SBVOmkIprK/1TMiT8OI5pXJxcCnVPXzzvIVwCJVvTGszJPAj1T1VWf5BeDrwPS+6opIg6pmh+2jXlUPe99ksFcmAJ2d7XS0NIHHS0ZmDqpKIBC6olBPIl6vh6AqgW4lLclHdyCA1xdqTWysqyYzexzi8RAMKs0dAbJSEw7uu6O9le5AV8SeV20tjSQkJpOQmHTEGIPd3TQ21pOTa5NVGzNSOtpaCAa7SU3P6lf5hjY/mckJeDxyyPpgUGnxB8hM/vC7oLOjjUCX/+B3QVNHF6kJXnxeDy1N9SQmpZCQkEhbaxNpGdk01teQkZlDa2szGZnZdAa66WiuJ+DvQDVIQlIqIoI/EKSuO5lZBRkMxuGuTKI5OVY5MClsuRjY388yiYepWykiRWHNXFXDGnUvSUkpJCWlHLKuJ1kcUs730W1ZufkH33s8ckgiAUhO6ftBjP39hQXweL2WSIwZYcmp6QMqn50auUXB45FDEglAUnIqSckfdrIJ356e+eHfyj3JJisn9O89IzO0nOTzkpQT+TtgpL4ZotnMtRaYKSLTRCQRuBRY3avMauBKp1fXYqDRabo6XN3VQE83iquAx0f6QIwxxhwqalcmqhoQkeXAs4AXuE9Vt4jI9c72VcBTwLlAKdAGXHO4us6ubwMeFZFrgb3AxdE6JmOMMSFRu2cymgzlnokxxsSr0dKbyxhjzBhlycQYY8yQWTIxxhgzZJZMjDHGDFlc3oAXkWpgj9txjIA8wObn7T87XwNj52tgxuL5mqKq+ZE2xGUyGatEZF1fPS3MR9n5Ghg7XwMTb+fLmrmMMcYMmSUTY4wxQ2bJZGy52+0AYoydr4Gx8zUwcXW+7J6JMcaYIbMrE2OMMUNmycQYY8yQWTIxxhgzZJZMxigRmS4ivxGRP7sdSywQkbkiskpE/iwiX3I7ntFORE4XkVecc3a62/HEAhE51Tlf94rI627HM9wsmcQQEblPRKpEZHOv9UtEZJuIlIrICgBV3aWq17oT6egwwPO1VVWvBy4B4magWbiBnC9AgRYgmdAMqXFpgL9jrzi/Y08AD7oR70iyZBJbHgCWhK8QES+wElgKzAMuE5F50Q9tVHqAAZwvETkfeBV4IbphjhoP0P/z9YqqLgW+Adwa5ThHkwcY+L/JzwEPRyvAaLFkEkNU9WWgrtfqRUCpcyXiBx4BlkU9uFFooOdLVVer6seAy6Mb6egwkPOlqkFnez2QFMUwR5WB/o6JyGRC05E3RTfSkWfJJPZNBMrClsuBiSIyTkRWAceLyDfdCW1U6ut8nS4id4rIrwlNH21C+jpf/+qcq4eAX7oS2egV8Zw5768F7o96RFEQtTngzYiRCOtUVWuB66MdTAzo63z9A/hHdEOJCX2dr78Cf412MDEi4jkDUNXvRjmWqLErk9hXDkwKWy4G9rsUSyyw8zUwdr4GLi7PmSWT2LcWmCki00QkEbgUWO1yTKOZna+BsfM1cHF5ziyZxBAReRh4A5gtIuUicq2qBoDlwLPAVuBRVd3iZpyjhZ2vgbHzNXB2zj5kD3o0xhgzZHZlYowxZsgsmRhjjBkySybGGGOGzJKJMcaYIbNkYowxZsgsmRhjjBkySybGhBGRbhHZEPaz4si1Rp6EvCgimc5ygYj8QUR2ich6EXlDRC44wj4+EJHZvdb9j4h8XUTmi8gDI3gIZoyzZ3MZc6h2VV0wnDsUEZ8zkG0ozgXeVdUmERHgb8CDqvo55zOmAOcfYR+PEBqNfatTxwNcBHxcVfeISLGITFbVvUOM1cQhuzIxph9EZLeI3Coib4vIJhGZ46xPcyZIWisi74hIz6PGrxaRP4nI34HnRCRVRB4VkY0i8kcReVNESkTkWhH5WdjnfEFEfhohhMuBx533ZwJ+VV3Vs1FV96jqL5x9eEXkDiemjSLyRafYw4SSSY/TgN2qusdZ/nuv7cb0myUTYw6V0quZ67Nh22pUdSFwF/A1Z93/BV5U1ROBM4A7RCTN2XYycJWqngl8GahX1WOB7wMnOGUeAc4XkQRn+RoiP6L848B65/3RwNuHOYZrCc2ZcSJwIvAFEZmmqhuBoIgc55S7lEMnaVoHnHqY/RrTJ2vmMuZQh2vm6nnk+nrgX5335xBKBj3JJRmY7Lx/XlV7Jk46Bfg5gKpuFpGNzvtWEXkR+LSIbAUSVHVThM/OVdXmSEGJyEpn/34ngZwDHCsiFzlFsoCZwAc4VycisoXQhE3fCdtVFTChj2M35rAsmRjTf53Oazcf/tsR4EJV3RZeUEROAlrDVx1mv/cC3wLep++JkwIi4nFmONwCXNizQVVvEJE8QlcWPZ91o6o+G2E/DwPPAf8ENqpqVdi2ZKD9MHEa0ydr5jJmaJ4FbnRuiiMix/dR7lXgEqfMPGB+zwZVfZPQ/BeHmxt8GzDdef8ikCwiXwrbntorpi/1NJ2JyKyepjdV3QnUArdF+KxZwOY+j9SYw7BkYsyhet8zue0I5b8PJAAbRWSzsxzJr4B8p3nrG8BGoDFs+6PAa6pa30f9J4HTITTNIfAZ4BNOd9+3gAed/ULoSuc94G0npl9zaCvEw8Ac4LFen3GG8znGDJg9gt6YKBARL6H7IR0iMgN4AZilqn5n+xPAz1T1hT7qFwG/VdWzRyi+JEJNX6cMQzdmE4fsnokx0ZEKvOQ0PQnwJVX1i0g28BahMSQREwmAqlaIyD0ikqmqTSMQ32RghSUSM1h2ZWKMMWbI7J6JMcaYIbNkYowxZsgsmRhjjBkySybGGGOGzJKJMcaYIbNkYowxZsj+P+BYfivUOl0YAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.semilogx(EGeV, pa[0], lw=2)\n", "plt.semilogx(EGeV, pa_many_cells[0], lw=2, ls='--')\n", "\n", "plt.xlabel(\"Energy (GeV)\")\n", "plt.ylabel(\"$P_{a\\gamma}$\")\n", "\n", "print(\"Absolute maximum difference:\", np.max(np.abs(pa[0] - pa_many_cells[0])))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As expected the difference is very close to zero and can be explained by numerical accuracy. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Mixing in the ALPS II experiment\n", "\n", "We can even use `gammaALPs` to calculate the mixing in case of the ALPS II experiment! ALPS II is a light-shining-through-a-wall experiment where ALPs are produced by immersing a strong laser in a magnetic field. Some of the photons should convert to ALPs, pass a light-tight barrier, and reconvert in a magnetic field behind this wall back to photons. \n", "\n", "The foreseen laser in the ALPS II experiment has a wavelength of 1064 nm, and the production cavity has a length of roughly 100m. To be precise, twelve refurbished HERA dipole magnets with a length of 8.8m will be used in the experiment which have a field strength of 5.3 T. The laser beam will be polarized parallel to the magnetic field to maximize the mixing. \n", "\n", "We provide these specs below:" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "from astropy import units as u" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Laser energy in eV: 1.1652650228684234 eV\n", "53000000000.0 3.42226292965325e-18\n" ] } ], "source": [ "# center energy around 1064 nm\n", "wavelength = 1064 * u.nm\n", "energy = ((c.c / wavelength.to('m')) * c.h).to('eV') \n", "print(\"Laser energy in eV:\", energy)\n", "\n", "# define an array centered around 1064 nm\n", "EGeV = np.logspace(np.log10(energy.value) - 4,np.log10(energy.value) + 14.5, 2000) * u.eV.to(\"GeV\")\n", "\n", "# B field and length of production cavity\n", "B = 5.3 * u.T.to('1e-6G')\n", "L = 8.8 * 12 * u.m.to('kpc')\n", "print(B, L)\n", "\n", "# fully polarized beam\n", "pin = np.diag((1., 0., 0.))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We set up the module list as before and run the calculation:" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "\u001b[0;36menvirons.py:\u001b[0;35m 257\u001b[0;0m --- \u001b[1;36m\u001b[1;36m\u001b[1;36mINFO\u001b[1;0m\u001b[1;0m\u001b[1;0m: Using inputted chi\n", "\u001b[0;36menvirons.py:\u001b[0;35m 270\u001b[0;0m --- \u001b[1;31m\u001b[1;31m\u001b[1;31mWARNING\u001b[1;0m\u001b[1;0m\u001b[1;0m: r_abell <= L0: assuming one domain from 0. to L0\n", "\u001b[0;36m core.py:\u001b[0;35m 652\u001b[0;0m --- \u001b[1;36m\u001b[1;36m\u001b[1;36mINFO\u001b[1;0m\u001b[1;0m\u001b[1;0m: Running Module 0: \n", "/Users/manuelmeyer/Python/gammaALPs/gammaALPs/base/transfer.py:802: UserWarning: Not all values of circular polarization are real values!\n", " warnings.warn(\"Not all values of circular polarization are real values!\")\n" ] } ], "source": [ "ml_alps = ModuleList(alp, src, pin=pin, EGeV=EGeV, seed=0)\n", "\n", "ml_alps.add_propagation(environ='ICMCell',\n", " order=0, # order of the module\n", " B0=B, # B field strength\n", " L0=L, # cell size\n", " nsim=1, # one single realization\n", " n0=1e-100, # electron density\n", " r_abell=L * 1.1, # full path, chosen that we only have a single cell\n", " beta=0., \n", " eta=0.\n", " )\n", "\n", "ml_alps.modules[0].psin = np.ones_like(ml_alps.modules[0].psin) * np.pi / 2.\n", "\n", "# set the ALP parameters\n", "ml_alps.alp.m = 1e4\n", "ml_alps.alp.g = 2.\n", "\n", "# run the calculation\n", "px_alps, py_alps, pa_alps = ml_alps.run()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We want to compare our results against the simple formula for a two-state mixing system given by \n", "\n", "$$P_{a\\to\\gamma} = \\Delta_{a\\gamma}^2 \\frac{\\sin\\left(\\Delta_\\mathrm{osc} L /2 \\right)^2}{\\left(\\Delta_\\mathrm{osc} L /2 \\right)^2} $$\n", "\n", "For this we first calculate the values of the momentum differences $\\Delta$ and then calculate $P_{a\\to\\gamma}$." ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [], "source": [ "from gammaALPs.base.transfer import Delta_a, Delta_ag, Delta_pl, Delta_QED" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [], "source": [ "# calculate the Deltas\n", "DQED = Delta_QED(ml_alps.modules[0].B[0], EGeV)\n", "Dpar = Delta_pl(ml_alps.modules[0].nel[0], EGeV) + 3.5*DQED\n", "Dpl = Delta_pl(ml_alps.modules[0].nel[0], EGeV)\n", "Dag = Delta_ag(ml_alps.alp.g, ml_alps.modules[0].B[0]) * np.ones(EGeV.size)\n", "Da = Delta_a(ml_alps.alp.m, EGeV)\n", "\n", "# calculate the Dosc\n", "Dosc = np.sqrt((Dpar - Da)**2. + 4.*Dag**2.)\n", "\n", "# calculate the analytical expression\n", "Pag_analytical = Dag ** 2. * (np.sin(Dosc * L / 2.))**2. / (Dosc / 2.)**2." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And the comparison is plotted below. We find very good agreement between the code and the simple analytical expression." ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0, 0.5, '$P_{a\\\\gamma}$')" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEVCAYAAADkckIIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAxgklEQVR4nO3dd5xcZb3H8c9v+szWZHfTE0ggBEILELpoEPACUlS6gBdEIyAqigUQyxUVvYoiIE3ACHpD04uAQYVALiUECRBaQkIggTSS3WSzffrv/jEDbku2ZPc8s7O/9+u1r+zMc/ac78nOzm+eU55HVBVjjDHmAz7XAYwxxhQWKwzGGGM6sMJgjDGmAysMxhhjOrDCYIwxpgMrDMYYYzoomsIgIneKyCYReX2A1vd3EdkqIo90ev5pEVmS/1ovIg8OxPaMMaZQFE1hAOYAxw7g+n4BnNv5SVU9QlVnqOoM4DngLwO4TWOMca5oCoOqPgVsaf+ciOyS/+T/Yv6T/u59WN98oGlb7SJSBnwceLCfkY0xpiAFXAcYZLcBF6rqWyJyMHATuTfzgfBpYL6qNg7Q+owxpiAUbWEQkVLgMOB+Efng6XC+7TPAj7r5sXWq+h+93MRZwO07mtMYYwpN0RYGcofJtubPBXSgqn9hB84NiEgVcBC5XoMxxhSVojnH0Fn+EM8qETkNQHL2HaDVnwY8oqrxAVqfMcYUjKIpDCIyl9xVQtNEZK2IXACcDVwgIq8AbwAn92F9TwP3A0fl19f+ENOZwNyBS2+MMYVDbNhtY4wx7RVNj8EYY8zAsMJgjDGmg6K4Kqm6ulp33nln1zGGrOXLlwMwbdo0x0mMMV558cUX61S1pru2oigMO++8M4sXL3YdY8iaNWsWAAsWLHCawxjjHRF5d1ttRVEYzI656qqrXEcwxhQQKwyGo48+2nUEY0wBsZPPhiVLlrBkyRLXMYwxBcJ6DIZLL70UsHMMxpgc6zEYY4zpwAqDMcaYDuxQUhGKtzbTsGUjyXgr6UQr6WScrZV7kspAuHEVoeZ1qCpoFlRp3lpLOlTJ/GUbKalfRqRlfYf1qc9P3dhZAJRveY1w26YO7Rl/mC1jPgJAxeYlhOKbO7SnAzHqRx8KQGXtYoLJhg7tqVA5W2sOBGDEpucJpJo7tCfDI2io3h+AkRsX4k+3dWhPRKpprMqNj1j1/tP4MsmO/x/R0TSN3AuA6vVPIprt0N5WMo7myj1AlZr1T3T5/2wtnURLxVQkm6J6w1Nd2lvKJtNaPgVfOk7Vxme7tDdXTKWtdBL+VDMjNz3fpb2pcnfiJeMJJBsYUdv1suvGEXuSiI0hGN9C5eaXu7Q3VO1LMlJNqG0TFVte69K+tXp/UuERhFvfp7z+jS7t9TUHkQ6VEWlZS9nW5V3at4w6hEywhGjTu5Q2ruzSvnnMR8j6w8Qa36akaXWX9rqxH0N9AUoaVhBrXtOlvXb8UQCU1i8l2rqhQ1ufX3uJLYhPQPyI+NBgCc2jZyICpfVvEkg34xMf+Hy59lApqZHT8PuEimw9pbEYpeUj8AeG91vj8N77IS4Rb2XFvx6jefkT/CFwOss3pzi96S4ukj8T6bTs7vHfEyfM9wJ3c0Hg0Q5tpVtbWJTdgwv+sJifBW7jzMCCDu2NGuXTiTsAuCF4PSf6F3Vof19HcEritwDcGfxvPu5f0qH97exYTk1eC8C9oR8xw/dmh/ZXslP4YvLHAPwtdBV7+jpeXr0wM50vpnKX1D4Z+g6TfRs7tD+WOYAvpi4D4IXwN6iRjoXnwcxhXJq6BICl4a8Sk0SH9j+lj+K76QsQsqyKXERnt6Y/yTXpsymjlde6af9V6lSuz3yG0Wzh+cglXdqvTp3DHZnj2UXWMT/8rS7tl6e+wD2Zj7OPvM1D4e91af9K8hIezh7Gob43mBv6SZf285Pf4snsfhzjW8zvQr/q0n564nv8S/fgZN8z/CZ0U5f2TyZ+yhu6M2f7H+cnwTu7tM9KXMtqHcts/8NcGew6duSB8ZuopZKvB+7na4H/7dK+R/xO2oh0+9rLqjAl8ScArgn8jrMCT3Zo39HX3srsOE5O/hLIvfYO7vTaW5Kdwqc+fO1dwfj8a6+BEjYFxvF2+cGs3e8yDt2liuljy2k3t0tRK4pB9GbOnKnD6Qa3+toNLP/z1Ux//0HKaSGlfi4su57I2OkcJEuZllyKlFThD8fwBSP4QlFaJnwMfzBEtOk9QvGNCD4QARFefnUpiZHT2G/mwYSa1xBIbO34ByA+4lV7AhBsfA9/p0/86guSGJmbNTXUsBpfqqlTe4jEyNxd1aGt7+BLt3RozwaiJCt3BSBc/xaSiXdqLyFZOSXXvmU5ku3YI8iGykiW75xr37wM0XSH9kyoglT5JAAim9/I9ZTat0dGkiodD6pENr/e5f87HakmXToWshkiW5Z2bY+NIh0bjWSShOu7fuJOxcaQidUg6TjhrW91bS8dTyYyEkm1Em54u5v2iWQilfiSzYQaV3VpT5bvRDZUji/RQKjpvW7aJ5MNleKPbyXYzSf2ROWuaCCKv20zwU69xVz7VDQQwd9aS7D1/S7t8ZG7gy9IoOV9Am213bRPB5+fQPMGAvG6ru3VewMQbFqLP1HfsbGPrz1JNqLZLJrNoJol6wvRPHIvsqpE617Dn2hANZtbRrOkAyXUVx1AKpOlYtU8fE3r0HgjvrbNRJvf5c3UWL7d8lmELHNLr0OmHc/MT3+tKHoUIvKiqs7sts2rwiAiEeApcrOoBYAHVPUHnZYR4DfA8UArcJ6qvtTTuodTYXj2mSfY8/FzKdVWXin7KP79zmK3g48jVlrhOpoxRWlDQxuLXl3GXv83m6npt1jp34Xo2XczfsqerqPtkO0VBi9PPieAj6vqvsAM4FgROaTTMscBU/Nfs4GbPcxX8OY8u4pzH2lhfvRY1p75GAd886/MOOrMHS4KCxcuZOHChQOU0pjiMrYiyqeP2J9dr/wXLx50HdWZjZTedQxvv1q8fzOe9Yc01zX54KxiMP/VubtyMnBXftlFIlIpImNVdQPD3AvzH+C3j8U5evrufPKs24gE/QO27iuvvBKw+xiM2R7x+Tjg+PNZt/tBBO46kcRfLmHDpKcZW1niOtqA8/RyVRHxi8gSYBPwmKp2vkRjPND+IOja/HPdrWu2iCwWkcW1tV2PaxaTuvXvMu3pr3Jd+VyuP2u/AS0Kxpi+GT9lT9rO/isXZK/ksvtfI5sd+udpO/O0MKhqRlVnABOAg0Rkr06LdHfKv9v/dVW9TVVnqurMmppuR44tGu/O/RohTTLptJ9aUTCmAOw8dW++8smD+NfbG3miCA/DOrnBTVW3AguAYzs1rQUmtns8Aeh6mcQwsmrpCxzQ9CQvTzyXiVP3dR3HGJN35oET+X3ZrUyffx6pZKLnHxhCPCsMIlIjIpX576PA0cCbnRZ7CPic5BwCNAz38wubH72GFo0w/dOXu45ijGnH5xPKDjqHcbqJVx//o+s4A8rLi3HHAn8QET+5gnSfqj4iIhcCqOotwDxyl6quJHe56vke5is4W5oTrG1IkBh7GodXjR607Vx33XWDtm5jitneR57OuoU/JPbyHXD8Ba7jDBgvr0p6Fdivm+dvafe9Al/2KlOhe+iV9fwweTGPfuqIQd3OjBkzBnX9xhQrfyDAe5PP4NB3rmf96uWM27k4pse1QfQK2FMvvc70seXsMbZ8ULfz+OOP8/jjjw/qNowpVjsdcTYAq565z3GSgTP07+suUls2reP2unP5v10uAwa3x/DjH+fGirGZ3Izpu3GTd+cbldfxTv0uHO46zACxHkOBenvhg/hEmbDPLNdRjDE92Gmvj/DK+ia2tiZ7XngIsMJQoOTt+dRRyS57H+Y6ijGmBx+dGOQH/jmsWPiw6ygDwgpDgZrQ9Arvls7A57cb2owpdHtNHstp/qfQZQ+5jjIgrDAUoPfXrGQMdaTGH+Q6ijGmF4KhMCujezO6vsfBoIcEO/lcgF7amOWfyYu5cMZJnmzv1ltv9WQ7xhSz1poZ7PXe7bQ0baWkrNJ1nB1iPYYC9MKGFP/wf4xdd/NmvPdp06YxbVpxXH9tjCuxyTPxi/LuG4t6XrjAWWEoQJF3HuPomnoCfm9+PQ8//DAPP1wcJ82McWXCnoezXkeydv0611F2mB1KKjCazXJR/X+zrOpo4CxPtnnttbn5mE888URPtmdMMaoaPZHDo3dwQGIEn3AdZgdZj6HAbFz3DuW0wKihPW2gMcPRtDFlrNjY1POCBc4KQ4F5f0Vu7uqKnbsMK2WMKXCfkQXcWH8h6dTQvtHNCkOBaVv7KgDjph3gOIkxpq9GlUfYVdax4d3lrqPsECsMBcZfv4paRlBeWeU6ijGmjyp32huAundecZxkx9jJ5wJza+hcRlSdwC893Obdd9/t4daMKV7jds3NshjfsNRxkh1jhaHAvFofYtZunafCHlwTJ07seSFjTI9Ky0ewiZEEtqx0HWWH2KGkAtLS1MAZrfewX3Sjp9u99957uffeez3dpjHFanHsCFZkxruOsUOsMBSQ91cv45vB+5nm8/YGmZtvvpmbb77Z020aU6we3/kb3Jg6wXWMHWKFoYA0rHsTgIrxuzlOYozpr0kjY7zf2EYilXYdpd+sMBSQxMa3ABgzebrjJMaY/jo4sYilofPZtPoN11H6zQpDAfE1rqGeMkrLR7iOYozpp8rq0UQlyZZ1Q/cEtBWGAhJu3chm/yjXMYwxO6B6Ym6k4vimtx0n6T+7XLWAfDd8BbvWZPiNx9t94IEHPN6iMcWravREkhogW/+u6yj95lmPQUQmisiTIrJMRN4Qka91s8wsEWkQkSX5r+97la8QrG9MUjbS+x5DdXU11dXVnm/XmGLk8/up840k2OLtZecDycseQxq4TFVfEpEy4EUReUxVO98i+LSqDu1rvfoh3trMZclbiPJZYG9Ptz1nzhwAzjvvPE+3a0yxWhA7liZfNTNdB+knz3oMqrpBVV/Kf98ELAOG9l0gA6hu/SrOCcxnEt5/ypgzZ86HxcEYs+OeHX8+96Y/5jpGvzk5+SwiOwP7Ac9303yoiLwiIo+KyDYnJRCR2SKyWEQW19bWDlZUzzRszB2PjFbb8BTGDHWjS8O0Ndai2azrKP3ieWEQkVLgz8ClqtrYqfklYCdV3Re4AXhwW+tR1dtUdaaqzqypqRm0vF5prXsPgIrROzlOYozZUcc0PsBzvi/Q1FjvOkq/eFoYRCRIrij8SVX/0rldVRtVtTn//TwgKCLD4qxoeutaAKrHTXGcxBizo4KVYwGo37DabZB+8vKqJAHuAJap6q+2scyY/HKIyEH5fJu9yuhSqqWBzVQQLSlzHcUYs4NiVblDwo217zlO0j9eXpV0OHAu8JqILMk/dyUwCUBVbwFOBS4SkTTQBpypquphRmfmln+elW2n8ZiDbc+bN8/BVo0pXhWjJwHQtnmt4yT941lhUNVnAOlhmRuBG71JVFg2NyepKgs72XYsFnOyXWOK1ch8Ycg0bXKcpH9sSIwCcfbm33By+h9Otn3TTTdx0003Odm2McUoWlLGDdnTWBEcmgNiWmEoELNSTzEl6+Z45H333cd9993nZNvGFKv7Sj/LEt8ermP0ixWGApBMxKmghWzJ0L/s1hiTs1M0gTTYyWfTTw11G6gBfKVWGIwpFt9svY6yrZuAk1xH6TPrMRSAxs3rAQiWj3acxBgzUJLhkZRltrqO0S9WGApAY1Mza7Wa6IhxrqMYYwZIJlZNpTYMyWEx7FBSAVgd25NPJ67nySmHONn+ggULnGzXmGImJdWEJENDwxYqRgytARysx1AANjcnAagqDTlOYowZKIGy3NwqjXXrHCfpOysMBWDSyj/xu9CvKQu76cD98pe/5Je//KWTbRtTrHTCgVye+gK1mVLXUfrMCkMBqGxYyj6+d8gPE+W5Rx55hEceecTJto0pViVjduOezMfZlB56IwtYYSgAocQWmv0VrmMYYwZQZUTYU1aR2DL0xkuywlAAIulG2gLlrmMYYwZQZSjL38LfZfTqh1xH6TMrDAUgmmkkGbTCYEwxicbKSGoA4ltdR+kzu1y1ALyno0nGdnG2/Wg06mzbxhQr8floklLECoPpK1Xl88lvcsHkKRzlKMOjjz7qaMvGFLdmXxnB5FbXMfrMDiU51pbKkMoolbGg6yjGmAHW6i8jlOo8tX3hs8LgWOOmNTwa+g7TGxc6y3D11Vdz9dVXO9u+McXqb1Xnc1fwTNcx+swKg2Mt9RvZw7eG0qC78VTmz5/P/PnznW3fmGK1YeTBPJue5jpGn1lhcCzeuBmAUNlIx0mMMQNtJ18t+7Q97zpGn1lhcCzZXAdAtNzmYjCm2BzY9Bg3y89IJROuo/SJFQbHUs1bACipHFqjLxpjeuaLjQCgsb7WcZK+scLgWL2W8VxmOmWV7noMVVVVVFVVOdu+McXKX5I7RNzSUOc4Sd/YfQyOLSn9CLdnx7Gi1N2dz3/+85+dbduYYhbKF4bWIVYYPOsxiMhEEXlSRJaJyBsi8rVulhERuV5EVorIqyKyv1f5XNnamqIiGnI2sqoxZvBEKnKHiBNNVhi2JQ1cpqp7AIcAXxaR6Z2WOQ6Ymv+aDdzsYT4nTlh9DTfpT5xmuOKKK7jiiiucZjCmGIXH7MHZySt4L7qX6yh94tmhJFXdAGzIf98kIsuA8cDSdoudDNylqgosEpFKERmb/9miVJFYR0DcXrHw3HPPOd2+McWqrGIEz2b35phsiesofeLk5LOI7AzsB3S+wHc8sKbd47X557pbx2wRWSwii2trh9YZ//bCmRYSgaE3w5MxpmelIT/H+xYRrnvDdZQ+8bwwiEgp8GfgUlXtPIhIdwfatbv1qOptqjpTVWfW1AzdewAimRbSQSsMxhSjUNDPtcFbmLxhnusofeJpYRCRILmi8CdV/Us3i6wFJrZ7PAFY70U2V6LaSsYKgzFFq0ViSLLJdYw+8fKqJAHuAJap6q+2sdhDwOfyVycdAjQU8/kFgP/LzqC2Ym+nGSZMmMCECROcZjCmWLVJjEBqaBUGL+9jOBw4F3hNRJbkn7sSmASgqrcA84DjgZVAK3C+h/k8l8pkuSw5m29M2M1pjj/+8Y9Ot29MMWvzlxJMNbuO0SdeXpX0DN2fQ2i/jAJf9iaRey2JNABlEbvP0JhilfTHCGVaXMfoExsSw6HWTat4Pfx5dt/8mNMcl156KZdeeqnTDMYUq/tHXcrPwl91HaNP7KOqQ21N9YyTOJFw2GmOJUuWON2+McWspXwKKzYNrTufrTA4lGiuByAYq3QbxBgzaKZlVlAeXwTOZnXvOzuU5FCytQGAUEml2yDGmEGzV8siruIOspmM6yi9ZoXBoXS+MERKKxwnMcYMFgmX4ROlpbnBdZRes0NJDtUFxnBPehZHVY52mmO33dxeLmtMMfNFcx/8Whq3UFYxNKbwtcLg0KroXvw8PZtllaOc5rjtttucbt+YYuaPVgK5i02GCjuU5FBLvA2/DyJB+zUYU6yCsdwkXInmrW6D9IG9Izl0xMpreSF0ofNJembPns3s2bOdZjCmWGUnHsIRiV/zfunurqP0mh1KcsiXaiYuEdcxWLFihesIxhSt0rJy1uhoGlN+11F6zXoMDgVSzbTJ0JrAwxjTN6W+FBf6HyL8/suuo/SaFQaHQulmEn4rDMYUs7Kwj8uD91CxqfO8ZIXLCoND4UwLSZu9zZiiFo2VkVVBk0NnID07x+DQ3/0fo7JiFDMc55gxw3UCY4qXz++nhTCSHDpDb1thcOj3meM5dvQY1zG47rrrXEcwpqi1ShRfauj0GOxQkkOBxBbKg91OaW2MKSJxieIfQoXBegyOpFNJng/M5rlNFwL7OM1yzjnnADaTmzGD5VsVv2REeTkHuA7SS1YYHGltaaIckHDMdRTWrl3rOoIxRS0bHUl92nWK3utXYRCRaaq6fKDDDCfxlsZcYQjZVUnGFLujkgvwx+uAQ11H6ZX+nmM4S0QuFxHrcfRTvCU3BK8/YoXBmGJ3QPJ5PtH2qOsYvdavwqCqPwQeBm4WkSPF9WA/Q1CytQkAf9hucDOm2GUDMSLZNtcxeq2/h5L+Cxiff/hx4CzARmHrg8bASH6ROp2P17gfWOvQQ4dG99aYoSobLCFKkRcG4ElVXTCQQYabhkANv818iqOqdnEdhWuuucZ1BGOKmgZLiGkczWYRX+HfJdDfhLNE5FERuV1EvtzbHxKRO0Vkk4i8vo32WSLSICJL8l/f72e+gpdsaWAcdZTYWRpjil+4DL8oifjQuJehv4WhElgE/ASY1oefmwMc28MyT6vqjPzXj/oXr/CNeO9RFka+Sllyk+sonHLKKZxyyimuYxhTtFbu8jmmxP9ISzbkOkqv9Pfz6pb8z27Kf98rqvqUiOzcz20WlWw8N25KtKTccRLYvHmz6wjGFLVoJEoWHy2JDFVD4ELE/l6V9CPgFuB6oGFAE8GhIvJK/lDVngO87oLxwUiL0VL3hcEYM7jGxFfyk8AdJDavdh2lV3osDCLynyJSJyJbROQuESkDUNX1qnqBqv56APO8BOykqvsCNwAPbifXbBFZLCKLa2trBzCCR5ItZFQIh6OukxhjBllFegtnB+aT2To0RhnoTY/he8AxwO7Au8BPByuMqjaqanP++3lAUESqt7Hsbao6U1Vn1tTUDFakQSOpVlqJDIkrFIwxOyYYyx0Z+OD+pULXm3MMjar6wZx03xORQZuGSETGABtVVUXkIHKFqygPgC8pPYInNldypesgwFFHHeU6gjFFLRwrAyDT1ug4Se/0pjCMFZHZwDLgTSDY342JyFxgFlAtImuBH3ywPlW9BTgVuEhE0kAbcKaqFuW41G8E9+L16MSCKAzf+973XEcwpqiFYxUApBNDY7Ke3hSGH5AbF/psYG+gVETmAa8Ar6rq3N5uTFXP6qH9RuDG3q5vKCttfpfJ/pTrGMYYD0RLy4lrkFQy6TpKr/RYGFT1tvaPRWQCuUKxN3A80OvCYP7tnNpfIZoBTnMdheOOOw6ARx8dOoN8GTOURCtHs3viD3x71DQOcx2mF/p8H4OqrgXWAvMGPs7wEcy2EQ9UuI4BQFvb0BnDxZihKBzw4fcJLYmhMSmDXRLjSCjbRjpgl6oaMxyICD8N/Z7d1j3oOkqv2Eg9jkSycTJ+97O3GWO8MYsXebdxaLzlWo/BkQhtZINWGIwZLuK+KIH00BhEb2iUryL0g/QF7D96Xw52HQQ44YQTXEcwpuglfDH86VbXMXrFCoMDqUyWv6YPZpeRu7mOAsA3v/lN1xGMKXopX5RQZmgUBjuU5EBrWxsHyzKqtN51FGOMR5qC1TRrxHWMXrHC4EB86/vcG76aqVufcR0FgFmzZjFr1izXMYwpanMnfo9vhYfGKANWGByIt+QG0vKFSxwnMcZ4JRby05rMuI7RK3aOwYFkW64wBCJljpMYY7xyaOPfOTr5D+Bo11F6ZIXBgQ+G3g1Eh8BUTsaYATEqvYGDeBHNZgt+uP3CTlek0vHc0LvBqPUYjBk2QqX4RIm3Ff69DNZjcGBj6XS+mPwG36mZ6joKAKeffrrrCMYUPcmfU2xrbiBaUtgfCq0wOLBVRvBYdiY/KBvpOgoAF198sesIxhQ9Xzh36DjeWvhzMtihJAf89W8zy/cysQIpy62trbS2Do0bb4wZskpqeCs7nrZkwnWSHllhcGDchn8yJ/QLYoHCmJzu+OOP5/jjj3cdw5ii1rrTxzkm+Qu2Ria5jtIjKwwuJFtIq49w2IbdNma4iIX8ALQmC39OBisMDkiylVYiBX/JmjFm4FTE13J/6IdE1j7rOkqP7J3JAV+6hbgMjTFTjDEDIxbwcaBvBdK0wXWUHllhcMCXbiNhhcGYYSVSWg5AJl74VyUVyHUxw8v/lp1NXOq51nWQvPPOO891BGOKXiSWu3dBk1YYTDfeZjzZ2HjXMT5khcGYwRfNFwYShX/ns6eHkkTkThHZJCKvb6NdROR6EVkpIq+KyP5e5vPK9Obn2FeXuY7xobq6Ourq6lzHMKao+QMBFuvubPVVuI7SI6/PMcwBjt1O+3HA1PzXbOBmDzJ57j+b7+C4lgddx/jQqaeeyqmnnuo6hjFF78Lgj3m68lOuY/TI08Kgqk8BW7azyMnAXZqzCKgUkbHepPNOOBsnE4i5jmGM8Vh0iMzJUGjnGMYDa9o9Xpt/rvCv7+qDCHGyVhiMGXZ+nLqW7LqRwN2uo2xXoRUG6ea5bseNEJHZ5A43MWlS4d9i3l5M42jIZm8zZrgZo3Wk7ORzn60FJrZ7PAFY392Cqnqbqs5U1Zk1NTWehBsIyUSckKTRoPUYjBluUv4ogUzcdYweFVqP4SHgEhG5BzgYaFDVojqM1JYWTkr8jP+ctD+Hug6Td9FFF7mOYMywkPLHKEnVu47RI08Lg4jMBWYB1SKyFvgBEARQ1VuAecDxwEqgFTjfy3xeaElleVMnIWWFc079jDPOcB3BmGEhE4gR1jbXMXrkaWFQ1bN6aFfgyx7FcSLeUMu5/n9SlaoCCuPcyJo1ufP9EydO7GFJY8yO2Fi6O5ua4oxzHaQHhXaOoehltqzm6uAcatpWuY7yoXPPPZdzzz3XdQxjit4rE87msvQlrmP0yAqDx5JtuXFSAhG7KsmY4SYa9NOWypDJFsYkXdtihcFj6XgjAKFYueMkxhiv7b/5YRaFv0xbS4PrKNtlhcFj6XyPIRwtc5zEGOO1iC/LGKkn3tzoOsp2WWHw2AdjsYdipY6TGGO85o/k/u7jLYVdGArtPoai92bVMXwzUc5DVYVzBdBll13mOoIxw4IvnDtSkGi1wmDaadIw7+loopGw6ygfOvHEE11HMGZYCEZzF50k25ocJ9k+O5TksapNz/HFwKOE/IXzX798+XKWL1/uOoYxRc9XMYFHMgfTTGEPiVM4707DxM51/8dXAn9BpLvxAt340pe+xJe+9CXXMYwpev7R07gk9TXqYru6jrJdVhg85ku10EbEdQxjjAMlodzR+5YCn5PBCoPH/OlW4r6o6xjGGAdi6a28Gv4CE97+H9dRtstOPnsskGklKdZjMGY4ipWUEZVWNG4nn007wUwbSb/1GIwZjiLRErIqkGx2HWW7rMfgse+Xfp+xJX5ucB2knauuusp1BGOGBfH5aCWMpFpdR9kuKwwe25IKMSZaWOMkHX300a4jGDNstEkEX6qwp/e0Q0keO7N1LgfGn3Mdo4MlS5awZMkS1zGMGRb+GTySFaHprmNsl/UYPHZm5iGWthXWpWqXXnopAAsWLHCaw5jh4I+lFzAuGuEc10G2w3oMHtJslpjG0ZDNxWDMcFUS9tOWSLqOsV3WY/BQItFGRLJghcGYYeu7W39ILF0PvOA6yjZZj8FD8ZbctcsSsiG3jRmuMv4woWyb6xjbZYXBQ235Mdh9ESsMxgxX2UCMcDbuOsZ22aEkDzVFxnJ4/I/cMHlf11E6+OlPf+o6gjHDRjYYI4IVBpPXkkiTxUcsUlhDYhx22GGuIxgzbGSDJcS0sAuDHUry0qZlXB24kxGJda6TdLBw4UIWLlzoOoYxw8KmkTP5feZYUunCumy9PU8Lg4gcKyLLRWSliFzeTfssEWkQkSX5r+97mW+wSf0qzg08TimFddfjlVdeyZVXXuk6hjHDQu2Yj/Kz9Fm0prKuo2yTZ4eSRMQP/BY4BlgLvCAiD6nq0k6LPq2qJ3iVy0uZ/IiK4aidfDZmuCoNQjnNtMYTVESDruN0y8sew0HASlV9R1WTwD3AyR5u37l0PNdTCJdUOE5ijHFl19rHeDUym1TtStdRtsnLwjAeWNPu8dr8c50dKiKviMijIrLntlYmIrNFZLGILK6trR3orINC80PtRkoKaxA9Y4x3ApEyABKthTsng5eFobtJjrXT45eAnVR1X+AG4MFtrUxVb1PVmao6s6amZuBSDqJ0Ok1Cg8RKylxHMcY4EsgfSk4WcGHw8nLVtcDEdo8nAOvbL6Cqje2+nyciN4lItarWeZRxUD096rNc8NZhLA+GXEfp4LrrrnMdwZhhIxTLHTFIxxt7WNIdLwvDC8BUEZkMrAPOBD7bfgERGQNsVFUVkYPI9Wg2e5hxULUmMpSEC+/WkRkzZriOYMywEcr3GFJthXV1YnuevUupalpELgH+AfiBO1X1DRG5MN9+C3AqcJGIpIE24ExV7Xy4acjaf8M97CHvkbswq3A8/vjjgE3YY4wXgpVj+VXqVKZHJ7uOsk2efnxV1XnAvE7P3dLu+xuBG73M5KXJTS9Srut7XtBjP/7xjwErDMZ4IVZew/WZz/CDcOEWBrvz2UPBdDNxn93DYMxwFgv5GMNmMi1bXEfZJisMHgpnWkgFbC4GY4azkN/H0+FL2XP1Xa6jbJMVBg9Fsi2kgnapqjHDmfh8tEkESRXuyWcrDB5q1iht4WrXMYwxjsUJ4yvgwlB4104WsZNS1/D5KZM5wnWQTm699VbXEYwZVuK+KP50q+sY22SFwSOJdIZkJktZpPD+y6dNm+Y6gjHDSkKiBNLWYxj2mrfWcVfwGpJNFwG7uo7TwcMPPwzAiSee6DiJMcPDo+Wn0aIhCmsux3+zwuCRtq21fNT/Gi9og+soXVx77bWAFQZjvLJ05DG8XdvMFa6DbIOdfPZIW/NWAIKxSqc5jDHujQs2M6r1LdcxtskKg0cSLfUABGM2F4Mxw90ntt7L7alC7S9YYfBMqiV3CClcWuk2iDHGvUglUUkSL9CB9KwweKQl7eft7Fgi5VWuoxhjHJP8IeXm+sKcUcBOPnvkzbKDOSd5La+M2cV1lC7uvvtu1xGMGVb8JSMAaGmoo3rcTo7TdGWFwSP1rUn8PqG8AO9jmDhxYs8LGWMGTKhkJABtjYU53YwdSvLIjFW/43fhXyPS3Qynbt17773ce++9rmMYM2z4xu7N15MXURvqbtp79wrv42uRqmpaQZVscB2jWzfffDMAZ5xxhuMkxgwPpVXj+d/sERyhhXmVovUYBkE2q1ww5wW+8cdn0WwWgHCqnhZ/Yb4IjDHeqgj7mClvktmyGoC69e/y3o+ms2juT9wGy7PCMAheXlNP3fKFfP+t03jn9UUAxNKNJIJWGIwxUB4NcF/oaia8+78ArHxyDpOy6zjwzV+QzWQcp7PCMChqX/knx/n/RaW0UPv6EwCUZRtIhUc4TmaMKQT+QJAtUoG/+X0Agu8vAeCx7ExWr9/oMFmOFYZBMPKdv3JG4Cnep4bQ+n+RzSqvZXemqbywBs8zxriz1V9FOL4JgKrmFTwdOJQLU1/nxY3uewx28nkQjGhawbrwFJKBcka3vEldS4Lzk9/mv6bt6Tpatx544AHXEYwZdlpCNZQmN1LfkuQTbT/h8iPHEnlmCyvWbwHcXkJuPYYBlk4lmZR+l+bKPUiM2ofxupG3Vq8BYFxl1HG67lVXV1NdbTPLGeOleGw0IzKbeWXtVpIE2X2XXbgr+itOfuOrrqNZYRgoi+ZcyZL/Ppb3lr9MWFL4x+5NZNrR3JI+kU3P38+Toa+zE++7jtmtOXPmMGfOHNcxjBlW3pp0Fhcnv0rrC3/iysCf2Gt8OZSMYnzibdasfI3lPz6YJfPvcZLN08IgIseKyHIRWSkil3fTLiJyfb79VRHZ38t8/fH+mpU0bt3MIat/y4zW59j06M/JqlC160ym7PMRfpY+i9XrNjDZt5FRo0a7jtstKwzGeC8yYW8WZadTueoR/iP0CuXRENlRezKCJp74+4NMS79J2cKfU9uUYGNj3NNsnp1jEBE/8FvgGGAt8IKIPKSqS9stdhwwNf91MHBz/t+CkM1kaG5t4e0taRJP/wb8YQ5Z/nNeyu6KMpUDfG8Rb6zjYLmTRbsfgD8QYOpIPx9tfo5a3whqqse43gVjTIHYc3SEE30LOSyzmH+NOIGdgPKd94Pl8OR6H+/KuXyfuznyp3ezKTieG6csonzUJMrG7MqEvQ4nFhq8t28vTz4fBKxU1XcAROQe4GSgfWE4GbhLVRVYJCKVIjJWVQflluFr/uefHLPhFhBBNHclgGQzPFHxKV4L7M3Y5Co+t+UGStL11GRqiUmC5zIzmZP5D+aGrvtwPfdnPsaTJZ/kVzUPc8jaPzK1ugR/IPdfe03gdxzge4uXSz5CzWDshDFmSJo6upwbQjcC4NvlYwBM2OMgav9eThZh6pHnwDN38z+hn/Cb9Gc4cvXtsDr3sxc8cBnzswdw0r7j+Pkp+xAN+Qc0m5eFYTywpt3jtXTtDXS3zHigS2EQkdnAbIBJkyb1K9B7G+sY07wUIUuW3H9sRvw06mbe87UyMbuG0vQWFKFFYmyWKh4If4Z3fNO4Ifotjmr7O5umf55Vm/bgh4dNZkrVNGb/+RN8/YQDP9zGyE98m6UPbSQ66+v9ymiMKU6BYIhFUy8jvO459j76HADKK6t4bvLZfKF+GYcdeTmL3joLX6qZk0t8vLFxHwLZBOFsK4dH3mV1yUepa04QCQ78GQHJfTgffCJyGvAfqvqF/ONzgYNU9SvtlvkbcI2qPpN/PB/4tqq+uL11z5w5UxcvXjx44YvcrFmzAFiwYIHTHMaYvlHVfg/MKSIvqurM7tq87DGspePFuROA9f1YxgywefPmuY5gjOmHwRqt2curkl4AporIZBEJAWcCD3Va5iHgc/mrkw4BGgbr/IL5t1gsRiwWcx3DGFMgPOsxqGpaRC4B/gH4gTtV9Q0RuTDffgswDzgeWAm0Aud7lW84u+mmmwC4+OKLHScxxhQCz84xDCY7x7Bj7ByDMcPP9s4x2J3PxhhjOrDCYIwxpgMrDMYYYzqwwmCMMaaDojj5LCK1wLuuc/RTNVDnOsQOsn0oHMWwH7YP3thJVbsdqacoCsNQJiKLt3VlwFBh+1A4imE/bB/cs0NJxhhjOrDCYIwxpgMrDO7d5jrAALB9KBzFsB+2D47ZOQZjjDEdWI/BGGNMB1YYjDHGdGCFwRhjTAdWGAqUiHxKRH4nIn8VkU+4ztMXIjJFRO4QkQdcZ9kRIlIiIi+KyAmus/SHiBwhIreIyO0istB1nr7o7jWU/338If93cbbLfL21jf3YI/97eUBELnKZb1usMAwCEblTRDaJyOudnj9WRJaLyEoRuXx761DVB1X1i8B5wBmDGLdX+rJPqvqOql7gJum29eP38h3gPm9Tbl8ffw9Pq+qFwCPAH1zkbW8AXkOfAR7I/12c5FHsLnZ0P1R1Wf73cjpQkDfBWWEYHHOAY9s/ISJ+4LfAccB04CwRmS4ie4vII52+RrX70avyP+faHHq5T95H67U59P73cjSwFNjodcgezKHvv4fPAnO9Crgdc9ix19AEYE3++8wgZeyNOezg34KInAQ8A8wfvJj9Z4VhEKjqU8CWTk8fBKzMf4JIAvcAJ6vqa6p6QqevTfnpTX8OPKqqL3m9D531ZZ88D9dLfdyHI4FDyL2pflFECuJvpa+/BxGZRG6K3EZvk3Y1AK+hteSKAzh87xqIvwVVfUhVDwMK8pBYQbzYh4nx/PvTDuRe5OO3s/xXgKOBUz+Y/rQAdbtPIlIlIrcA+4nIFW6i9Vq3+6Cq31XVS4H/AX6nqlkX4Xppe6+tC4Dfe56o9/ryGvoLcIqI3Aw87HHOnvR6P0RklohcLyK3kpvOuOB4NuezQbp5bpt3F6rq9cD1gxdnQHS7T6q6GSjUYtbZdn8vqjrHuyj9ts19UNUfeJylr3r9GlLVFgp3Hvi+7McCYIEHmfrNegzeWQtMbPd4ArDeUZaBUgz7ZPvg1lDO3l6x7AdghcFLLwBTRWSyiISAM4GHHGfaUcWwT7YPbg3l7O0Vy34AVhgGhYjMBZ4DponIWhG5QFXTwCXAP4BlwH2q+obLnH1RDPtk++DWUM7eXrHsx/bYIHrGGGM6sB6DMcaYDqwwGGOM6cAKgzHGmA6sMBhjjOnACoMxxpgOrDAYY4zpwAqDGVZEJCMiS9p9bXf4c6/kB018QkTKt7PMHBH5UqfnPiUi80QkJCJPiYgNc2N2mBUGM9y0qeqMdl8/29EVDtCb8fHAKz2MgjqX3B217Z0JzM2P6DmfApi7wwx9VhiMAURktYj8l4i8JCKvicju+edL8hOzvCAiL4vIB8NZnyci94vIw8A/RSQmIveJyKsicq+IPC8iM0XkAhH5dbvtfFFEftVNhLOBv7Zb7hwR+Ve+V3Nrfrz/x4HdRWRsfpkYuRF4H8z/2IMU6DDOZmixwmCGm2inQ0ntP2HXqer+wM3AN/PPfRd4QlUPJDdHwy9EpCTfdijwn6r6ceBioF5V9wGuBg7IL3MPcJKIBPOPz6f7YbAPB16E3NSP5D75H66qM8hNSnO2qmbIDT19ev5nTgKeVNWm/OPXgQP7/l9iTEd2PNIMN235N9vu/CX/74vkppEE+AS5N/YPCkUEmJT//jFV/WDClo8AvwFQ1ddF5NX89y0i8gRwgogsA4Kq+lo32x7Z7g3+KHKF5QURAYgCm/Jtc4Ff5Ld1JnDXBytQ1YyIJEWkrN26jOkzKwzG/Fsi/2+Gf/9tCHCKqi5vv6CIHAy0tH9qO+u9HbgSeJNtT5qTFhFffkIgAf6gqt1NcvQsMFZE9gUOo+s5hzAQ304WY3pkh5KM2b5/AF+R/Ed3EdlvG8s9Q/4Qj+Tm+t37gwZVfZ7cWP3bm3t5OTAl//18cjP3jcqvb6SI7JRflwL3AX8A5qnqh0VARKqAWlVN9WM/jfmQFQYz3HQ+x9DTVUlXA0HgVRF5Pf+4OzcBNflDSN8BXgUa2rXfBzyrqvXb+Pm/AbMAVHUpcBW5k9qvAo8BY9stOxfYl9z5i/aOpECnijRDiw27bcwAyF81FFTVuIjsQu5T/275y0gRkUeAX6vq/G38/FjgLlU9Zgcy/AW4ovNhL2P6ys4xGDMwYsCT+auPBLhIVZMiUgn8i9w9Ct0WBQBV3SAivxOR8h7uZehWftawB60omIFgPQZjjDEd2DkGY4wxHVhhMMYY04EVBmOMMR1YYTDGGNOBFQZjjDEdWGEwxhjTwf8DT8g4CeG8j0gAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.semilogx(EGeV * 1e9, pa_alps[0])\n", "plt.semilogx(EGeV * 1e9, Pag_analytical, ls='--')\n", "plt.axvline(energy.value, ls='--', color='k')\n", "\n", "plt.xlabel(\"Energy (eV)\")\n", "plt.ylabel(\"$P_{a\\gamma}$\")" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.18" } }, "nbformat": 4, "nbformat_minor": 4 }