{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mean Free Path and Lorentz Invariance Violation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebooko demonstrates the calculation of the mean free path and optical depth with and without Lorentz invariance violation from a given EBL photon density"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from ebltable.ebl_from_model import EBL\n",
    "from ebltable.tau_from_model import OptDepth\n",
    "import astropy.units as u"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initiate the class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "ebl =  EBL.readmodel('gilmore')\n",
    "tau = OptDepth.readmodel('gilmore')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define some redshifts and energies for the interpolation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "z0 = 0.031\n",
    "steps_z = 50\n",
    "ETeV = np.logspace(-2,2,50)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate and plot mean free path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "gam = ebl.mean_free_path(z0, ETeV)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Mean free path (Mpc)')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAELCAYAAADOeWEXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAwoUlEQVR4nO3deXxU1fnH8c+ThCQkQFjCvgVJQBZRARcUV6yCiHtbcam2VGpbtFXbqrX9WdvaahetWtTiXm1BtLbu+wK4ExbZl8giYUvYQgghIcnz+2MGTFOSzCSZTCb5vl+vvMg9M/feh/NK5sk5555zzN0REREJVVy0AxARkdiixCEiImFR4hARkbAocYiISFiUOEREJCxKHCIiEhYlDhERCUtCtAOoj/T0dM/IyIh2GCIiMWXevHnb3L1zXc9vMonDzPoAfwW2Aavc/c7azsnIyCA7OzvisYmINCdmtr4+50e0q8rMHjOzPDNbUqV8rJmtNLMcM7s5WDwAeMXdvwMMjmRcIiJSd5Ee43gCGFu5wMziganAOAIJYqKZDQYWAJeY2bvAexGOS0RE6iiiicPdZwM7qhQfC+S4+xp3LwVmAOcB3wZuc/fTgfGRjEtEROouGk9V9QQ2VDrODZa9DlxnZg8B66o72cwmm1m2mWXn5+dHNFAREflf0Rgct0OUubsvAS6u7WR3n2Zmm4EJiYmJIxo8OhERqVE0Why5QO9Kx72ATVGIQ0RE6iAaiWMukGVm/cwsEbgEeDGcC7j7S+4+ucgTIxKgiIhUL9KP404HPgYGmlmumU1y9zJgCvAGsByY6e5Lw7zuBDObtr2olA9ztjV84CIiUi2L5R0Ak7pn+agbpvH6j08muVV8tMMREYkJZjbP3UfW9fyYXKvqQIujlcG67Xt54L2caIckItJixGTiODDG0Ts9FYAHZ31BTl5hlKMSEWkZYjJxHGhxlBUXMfHY3uwvd37+/BIqKmK3201EJFbEZOI40OJIS0vjprGHk94mkc/W7eC5ebnRDk1EpNmLycRxoMVRUFBA+5REfnlOYE3EO15dzrY9JVGOTkSkeYvJxFG5xQFw7pE9OCkrnYLi/fzuleVRjk5EpHmLycRRlZnx2/OHkpQQx/MLNmpuh4hIBDWLxAHQt1Mq143JAuDWfy9m3/7yKEckItI8xWTiqDzGUdnVJx1GVpc2rNu+l6ma2yEiEhExmTiqjnEckJgQx+8uPAKAB97/goUbdkUhOhGR5i0mE0dNjsnoyHdH96O8wvnRjAUUlZRFOyQRkWal2SUOgJ+OHcjh3dqyfvtebn8prPUTRUSkFs0ycSQlxHPfxKNJSohjZnYury7eHO2QRESajZhMHNUNjlc2oGtbbh0/CIBbnl/M5oLixgpPRKRZi8nEUd3geFVXHN+X0wZ2pqB4Pzc887nWshIRaQAxmThCZWb84eIjSW+TyMdrtvPwnDXRDklEJOY168QB0LltEn+8+EgA/vTmSpZsrL57S0REatdkEoeZnWRmD5nZI2b2UUNe+7TDu3DlqL7sL3eum7GA4lLNKhcRqatI7zn+mJnlmdmSKuVjzWylmeWY2c0A7j7H3a8BXgaebOhYbjl7EFld2rAmv0iP6IqI1EOkWxxPAGMrF5hZPDAVGAcMBiaa2eBKb7kUmN7QgSS3iufeSwKP6M6Yu4EZn33Z0LcQEWkRIpo43H02sKNK8bFAjruvcfdSYAZwHoCZ9QEK3H13JOIZ3KMdd1wQWJLk/15YqiVJRETqIBpjHD2BDZWOc4NlAJOAx2s62cwmm1m2mWXn5+eHffOLR/TiiuP7UlpewfefnqeNn0REwhSNxGGHKHMAd7/N3WscGHf3acDtwPzExMQ6BfDLcwYzom8HNhfsY8o/51NWXlGn64iItETRSBy5QO9Kx72ATeFcINQJgNVJTIjjgcuG07ltEp+s2cFdr6+o03VERFqiaCSOuUCWmfUzs0TgEuDFcC4QypIjtenaLpkHLhtOQpzx8Jy1vPR5WLlLRKTFivTjuNOBj4GBZpZrZpPcvQyYArwBLAdmuntUno89JqMjvzwn8EDXz55bxMothdEIQ0Qkpph77K7fNHLkSM/Ozq7XNdydG2d+zvMLNpLRKYUXpowmrXWrBopQRKTpMbN57j6yruc3mZnj4WiIrqpK1+KOC45gcPd2rNu+l2unL9BguYhIDWIycdR3cLyq1onx/O2KEXRMTWT2qnx+8/KyBrmuiEhzFJOJoyFbHAf07pjCtCtGkBgfx5Mfr+fJj9Y12LVFRJqTmEwcDd3iOGBkRkfuujgws/z2l5Yya1X4EwxFRJq7mEwckXTB0b249vRMKhym/GM+q7fqSSsRkcpiMnFEoquqsuvPGMD4I7pTWFLGd56cy3YtSyIiclBMJo5IdVUdEBdn/OnrR3JkrzQ27Cjme0/No6RMe3iIiECMJo7G0Doxnoe/NZLuaclkr9/JLf9aTCzPeRERaSgxmTgi3VV1QJd2yTxy5UhSEuN5fsFG7VkuIkKMJo5Id1VVNqRHGnd/4ygA7np9JZ+s2R7xe4qINGUxmTga29ih3bjmlP6UVzhT/rmArbv3RTskEZGoUeII0U/OHMCowzqxbU8JU/45n/1alkREWigljhAlxMdx38Sj6douibnrdnLna9rDQ0RapphMHI01OF5V57ZJB/fwePSDtby8SHt4iEjLE5OJozEHx6sa0bcjvxg/CICbnltETp5mlotIyxKTiSParjwhg3OP7EFRaTnXPD2fopKyaIckItJolDjqwMz4/YVHkNWlDTl5e/jZvxZpcqCItBhKHHWUmpTAg5ePIDUxnlcWbWZm9oZohyQi0ihqTRxmFmdmR5vZeDM73cy6RiKQ4H3uMLP7zezKSNyjoWV2acMdFwSWYf/1S8v4cvveKEckIhJ51SYOM+tvZtOAHOBOYCLwA+AtM/vEzL5tZjUmHjN7zMzyzGxJlfKxZrbSzHLM7OZg8XlAT2A/kFuP/1OjOu+oHow/ojtFpeXc+OxCyivUZSUizVtNH/y/BZ4G+rv7We5+ubtf7O7DgHOBNOCKWq7/BDC2coGZxQNTgXHAYGCimQ0GBgIfu/sNwPfr8p+JBjPjt+cPpUvbwPwOrWclIs1dtYnD3Se6+2w/xKivu+e5+1/c/cmaLu7us4EdVYqPBXLcfY27lwIzCLQ2coGdwfdUu4a5mU02s2wzy87Pbxo79HVITeSui4cBcPebq1i+eXeUIxIRiZxQxjh+aGbtKx13MLMf1OOePYHKI8m5wbLngbPM7H5gdnUnu/s04HZgfmJiYj3CaFinDezCZcf1obS8guufWaj9O0Sk2Qrlqaqr3X3XgQN33wlcXY972iHK3N33uvskd7/W3afWdIFoTgCsya3jB5HRKYUVWwq5+61V0Q5HRCQiQkkccWZ28MM+OEZRnz/1c4HelY57AWGt3RGtJUdqk5KYwN3fPIo4g2mz1/DZ2qq9dCIisS+UxPEGMNPMxpjZ6cB04PV63HMukGVm/cwsEbgEeLEe12tShvfpwA9OzcQdbnx2IXs0q1xEmplQEsdNwLsEnnT6IfAO8LNQLm5m04GPgYFmlmtmk9y9DJhCICEtB2a6+9Jwgm6qXVUHXDcmiyE92rFhRzG/fims/5qISJNnoSyVEWwZDAIqgJXBp6GixswmABMyMzOvXr16dTRDqdbqrYWMv/8DSssqeOjyEYwd2i3aIYmIAGBm89x9ZF3PD+WpqvHAF8C9wF+BHDMbV9cbNoSm3uIAyOrallvGHQ7Azc8v0q6BItJshNJV9WfgNHc/1d1PAU4D7olsWDVrqoPjVV11QganDOjMrr37uXHm51RoVrmINAOhJI48d8+pdLwGyItQPCGJhRYHBGaV//Hrw+iYmsgHOdt47MO10Q5JRKTeQkkcS83sVTO7Krj44EvAXDO70MwujHB8Ma9L22T+cFFgVvkfXl/J0k1Nu5UkIlKbUBJHMrAVOAU4FcgHOgITgHMiFlkNYqWr6oAzBnfl8uMDs8p/NGMhxaWaVS4isSukp6qaqpEjR3p2dna0wwhJcWk5E/76ATl5e7ji+L785vyh0Q5JRFqo+j5VlVDDhe+r6UR3v66uN22JWifGc+8lR3H+1A956pP1nDqwM2MGRWRrExGRiKqpq+oaYDSB5UCygXlVvqIm1rqqDhjSI42fnRV4RPdnzy0ir1CP6IpI7KkpcXQHpgFnEdh3oxXwors/Wdty6pEWK09VHcqk0f04MbMT24tKuXHm59r4SURiTk37cWx394fc/TTgKqA9gSesatu8SWoQF2f8+etH0SGlFXNWb+Put1ZGOyQRkbCEMnN8OPBj4HLgNaLcTdUcdEtLZuqlw4kzmPreF7y2eHO0QxIRCVlNe47fbmbzgBuAWcDI4H4ZyxotumbshMx0fn72IABufPZzVm0tjHJEIiKhqanF8UsC+4ofCfwemG9mi8xssZktapToqhGrg+NVTRrdj/OO6sHe0nIm/z2bguL90Q5JRKRW1c7jMLO+NZ3o7usjElEYYmkeR3WKS8u56MGPWLZ5N6cN7MyjVx5DXNyhNkkUEWkYkVwd90t3X1/dV/Dm+oSrp9aJ8fztihG0T2nFeyvzuedtbTkrIk1bTYnjPTO71sz6VC40s0QzO93MngSujGx4LUPvjin8dWJgsPz+d3N4fcmWaIckIlKtmhLHWKAcmG5mm8xsmZmtAVYDE4F73P2JRoixRRidlc7Nwf07bpy5kNUaLBeRJqqmeRz73P0Bdz8R6AuMAYa7e193v9rdFzZkIGZ2qpnNMbOHzOzUhrx2rLj6pMOYcGQPikrLuerxuWwp0MxyEWl6QlkdF3ff7+6b3X1XOBc3s8fMLM/MllQpH2tmK80sx8xuPnAbYA+B1Xhzw7lPc2Fm3HXRERzVuz0bdxXzrcc+ZdfeqO7SKyLyP0JKHPXwBIEur4PMLB6YCowDBgMTzWwwMMfdxwE3AbdHOK4mKyUxgcevOoasLm1YtXUP335iLntLy6IdlojIQRFNHO4+G9hRpfhYIMfd17h7KTADOM/dK4Kv7wSSqrummU02s2wzy87Pz49I3NHWITWRpyYdR8/2rVnw5S6+99Q8Sssqaj9RRKQRRLrFcSg9gQ2VjnOBnsEdBf8GPAX8tbqT3X2au49095GdO3eOcKjR0y0tmae/exydUhOZs3obN8xcqAURRaRJCGWtqgvNbLWZFZjZbjMrNLPd9bjnoeZ+uLs/7+7fc/dvuvv7tcTULGaO16ZfeipPfudY2iYl8PKizfzfC0uI5Y23RKR5CKXF8QfgXHdPc/d27t7W3dvV4565QO9Kx70I7PkhhzC0ZxoPXzmSxIQ4/vHpl9zzliYIikh0hZI4trr78ga851wgy8z6mVkicAnwYjgXiOX9OOri+MM6MfXS4cTHGfe9m8OD738R7ZBEpAWraXXcC83sQiDbzJ4xs4kHyoLltTKz6cDHwEAzyzWzSe5eBkwB3gCWAzPdfWk4QbeUrqrKvja4K3ddNAwzuOv1FfzxjRXqthKRqKhpkcPHazjP3f07kQkpdM1hkcNw/WfBRm58NrBz4JWj+nLbhCFaFFFEwlLfRQ4TqnvB3b8dvMGJ7v5hlZueWNcbNgQzmwBMyMzMjGYYUXH+0T1JTUrgh/+cz5Mfr6ewpIw/XDSMhPhoPCAnIi1RKJ8294dY1mha2hhHVV8b3JXHrzqGlMR4np+/kR/+cz4lZeXRDktEWoiaxjhGmdmNQGczu6HS16+A+EaL8NCxtbgxjqpOzEznqUnH0S45gTeWbuW7T2ZrhrmINIqaWhyJQBsC3VltK33tBi6OfGjVa+ktjgNG9O3AjMmjSG8TmCR4xaOfaW0rEYm4agfHD77BrG9T2O2vskpjHFevXr062uFE3Zr8PVz+yKdsKthHv/RUHrlyJP07t4l2WCLSRNV3cDyUxNEZ+BkwhMDKtQC4++l1vWlDaYlPVVVn065iJj2ZzfLNu2mbnMADlw3npKzmuySLiNRdJLeOPeAfwAqgH4FVa9cRmMQnTUiP9q157ppRnDWkK4X7yrjq8bk8+dE6zfUQkQYXSuLo5O6PAvvdfVZw/sbxEY5L6iA1KYEHLxvBlNMyKa9wbntxKb/4zxL2l2tlXRFpOKEkjv3Bfzeb2XgzO5rA+lJRo6eqqhcXZ/zkrIH85ZtHHVzf6srHNGguIg0nlMTxWzNLA24EfgI8Alwf0ahqoaeqanf+0T15ZvLxpLdJ4qMvtnP+1A9Zvrk+ixqLiATUOjjelGlwvHabdhVz9d+zWbppN0kJcdx+7hC+eUxvzLRMiUhLFfHBcTM7zMxeMrNtwf3DXzCzw+p6Q2lcPdq35l/fP4FLjulNSVkFNz+/mOufWUhRiSYLikjdhNJV9U9gJtAN6AE8C0yPZFDSsJJbxXPnRcO455tH0rpVPP9ZuIkJf/2AFVvUdSUi4QslcZi7P+XuZcGvp4Go9m9pcLxuLji6Fy9deyIDurZhTX4R50/9kJlzN+iRXREJSygTAO8EdgEzCCSMbwJJwFQAd98R2RCrpzGOuikuLef/XljCs/NyAbjw6J78+vyhtEmqdrFkEWlGGmPm+NoaXnZ3j9p4hxJH/Tw3L5df/Gcx+/ZX0KdjCvd88yhG9O0Q7bBEJMIinjiaMiWO+svJK+S66QtZtnk38XHGdadn8cPT+mt/D5FmrDGWHGk0ZpZqZvPM7Jxox9JSZHZpy79/eALfO/kwyiuce95exTenfcKX2/dGOzQRaaIimjjM7LHgI7xLqpSPNbOVZpZjZjdXeukmAk9wSSNKSojnlrMH8Y/vHke3dsnMW7+Ts++bw7/m5WrgXET+R6RbHE8AYysXmFk8gYH1ccBgYKKZDTazM4BlwNYIxyTVODEzndd/fBLjhnZjT0kZNz77OVOmL2BnkZYrEZGvhDIB0MzscjP7v+BxHzM7NpSLu/tsoOpTV8cCOe6+xt1LCTytdR5wGoHFEy8FrjazJtWN1lK0T0nkgcuG84eLh5GSGM8rizZz5l9m896KvGiHJiJNRCgfzg8Ao4CJweNCgo/i1lFPYEOl41ygp7vf6u4/JjDh8GF3P+SSrmY22cyyzSw7Pz+/HmFIdcyMb4zszWs/OoljMjqQX1jCt5+Yyy3PL2KPZpyLtHihJI7j3P2HwD4Ad99JYFvZujrUIkkHO9Ld/Ql3f7m6k919GoF9QeYnJtYnDKlN306pzJg8ip+ffTiJ8XFM/2wD4+6dzadrtkc7NBGJopCWVQ+OSzgc3BGwPhs85AK9Kx33AjbV43oSQfFxxuST+/PStaMZ0qMdG3YUc8nDn3DHK8vYt7882uGJSBSEkjjuA/4NdDWzO4APgN/V455zgSwz62dmicAlwIvhXEDLqje+gd3a8u8fnMi1p2cSZ8bDc9Zyzv0fsHDDrmiHJiKNLKQJgGZ2ODCGQDfTO+6+PKSLm00HTgXSCTwtdZu7P2pmZwN/AeKBx9z9jrCCNpsATMjMzLx69erV4ZwqDWDhhl3cMHMha/KLiDOYfHJ/fnxGFsmt4qMdmoiEoFFmjpvZaCDL3R8PdlW1cfealiJpFJo5Hj379pdz91ureHjOGtwhq0sb/vT1Izmyd/tohyYitWiM/ThuIzAx75ZgUSvg6bresCFoddzoS24Vz8/PHsRz15zAYemprM7bwwUPfMhdr6+gpExjHyLNWShjHBcA5wJFAO6+CWgbyaBqozGOpmNE3w68+qOTuPqkfjjw4PtfcM59GvsQac5CSRylHujPOvBUVWpkQ6qdWhxNS3KreG4dP5hnvzeKfsHWx4UPfMhvX15GcalaHyLNTSiJY6aZ/Q1ob2ZXA28DD0c2rJqpxdE0jczoyKvXncT3Tg6stP/IB2s56y+z+ShnW5QjE5GGVOPguJkZgXkWhwNnEniq6g13f6txwquZBsebrkW5u/jZc4tYsaUQgEuO6c0tZw8irXWrKEcmIo2xkdM8dx9R1xtEgh7HjQ2lZRX8bdYX3P9uDqXlFXRpm8Rvzx/KmUO6RTs0kRatMfbj+MTMjqnrDSJBXVWxITEhjmvHZPHqj0YzvE978gpLmPzUPL73VDabC4qjHZ6I1FEoLY5lwABgPYEnq4zAlrHDIh9ezdRVFTvKK5ynPl7HH99YSVFpOamJ8Vz/tQFcdUKGdhsUaWQR66oys37uvtbM+h7qdXdfX9eb1pe6qmLX5oJibn9xGa8v3QLA4O7t+N2FR3CUJg6KNJpIJo557j7CzN5x9zF1jjCC1OKIXe8s38r/vbCUjbuKMYPLj+vLT8cOpF2yBs9FIi2SiWMB8B/gu8A9VV9397vretOGosQR2/aWlnHvO6t5dM5ayiqc9DZJ3DR2IBcN70Vc3KFW3xeRhhDJwfFLCOzBkUBgpnjVL5F6SUlM4JZxg3j5utGM6NuBbXtK+Olzi7jgwY9Y8OXOaIcnItUIZXB8nLu/1kjxhEUtjuajosL5z8KN3PnaCvIKSwC4aHgvbho7kC7tkqMcnUjz0iir4zY1GhxvvvaUlDH1vRwenbOW0vIKUhPjuXZMFt8+MYOkBC3bLtIQWmTiOEAtjuZr3bYifvvKct5evhWAvp1S+MmZAxl/RHeNf4jUU2NMABRpdBnpqTxy5Uie/M6x9O+cyvrte7l2+gLOf+BDrX0lEmWhbuR0ApBBYKAcAHf/e+TCCo1aHC1DWXkFM7Nz+cvbqw6Of5w8oDM3jR3IkB5aPUAkXI2xVtVTQH9gIXBgjWx39+vqetNq7jMI+BGBbWbfcfcHaztHiaNl2VtaxuMfruOh97+gsKQMgPOP6sENXxtIn04pUY5OJHY0RuJYDgz2OgyGmNljwDlAnrsPrVQ+FriXwJ7jj7j7nZVeiwMedvdJtV1fiaNl2llUytT3cvj7x+spLa8gPs644Oie/PC0TPqlR327GJEmrzHGOJYAdV3O9AlgbOUCM4sHpgLjgMHARDMbHHztXOAD4J063k9agA6pifzinMG8+5NTuGh4LwCem5fLmD+/z49nLCAnrzDKEYo0b6G0ON4DjgI+A0oOlLv7uSHdwCwDePlAi8PMRgG/cvezgse3BK/3+0rnvOLu42u7tlocAvDl9r088H4Oz83LpazCMYOzh3ZnyumZDOreLtrhiTQ59W1xJNT+Fn5V14tXoyewodJxLnCcmZ0KXAgkAa9Wd7KZTQYmA/Tp06eBQ5NY1KdTCndeNIwpp2fy0KwvmDk3l1cWb+aVxZs5/fAufOfEfpyY2YnAvmQiUl+1Jg53n9XA9zzUb6+7+/vA+yHEM83MNgMTEhMTm9QGUxJdvTqk8Nvzj2DKaVk8NOsLpn/2Je+uyOPdFXkM7NqW74zO4LyjepLcShMJReqj1jEOMzvezOaa2R4zKzWzcjPbXY975gK9Kx33AjaFcwFt5CQ16ZaWzK/OHcJHN5/OT84cQJe2SazcWshN/1rMiXe+y91vriSvcF+0wxSJWaGMcWQTWPDwWWAk8C0gy91/HtIN/neMIwFYBYwBNgJzgUvdfWnIQWvJEQlDaVkFryzexKMfrGXJxsDfPK3ijTMHd+Mbx/RmdGY68ZqNLi1IY4xx4O45Zhbv7uXA42b2UYjBTQdOBdLNLBe4zd0fNbMpwBsEHsd9LJykIRKuxIQ4Lji6F+cf1ZO563by2AdreXPZloPjID3Skrl4ZG++PqIXvTtqPohIbUJpccwGzgAeAbYAm4Gr3P3IyIdXMz1VJXW1uaCY57JzmTlvAxt2fLX/+YmZnfjGyN6cMagrqUkh/V0lEnMaYwJgX2ArkAhcD6QBD7h7Tl1vWl/qqpKGUlHhfLJ2OzPnbuC1JVsoKasAICkhjlMHdubsI7ozZlBX2iiJSDPSKKvjmllroI+7r6zrjSJBLQ5pSAXF+3nx8038Z8FG5q3/aiOpxIQ4ThnQmbOP6MaYQV21va3EvMZocUwA/gQkuns/MzsK+HWoEwAjQS0OibQtBft4bclmXlu8hbnrd3Dg1yQhzhjetwOnDOjMyVmdGdKjnZZ5l5jTGIljHnA68L67Hx0sW+Tuw+p604aiFoc0hq279/HG0i28smgzc9ftoKLSr0yn1ERGZ6VzclZnRmel01W7FUoMaIynqsrcvUCzbqWl6toumW+NyuBbozIoKN7PRznbmL06n1kr89lUsI8XFm7ihYWBqUg927dmeN8OjOjTnuF9OzCoeztaxWvbG2leQkkcS8zsUiDezLKA64CQHseNlEpdVdEMQ1qgtNatGHdEd8Yd0R1354v8PcxatY3Zq/KZt34nG3cVs3FXMS99Hkgkya3iGNarPcN6pnF493Yc3q0tmV3aaPa6xLRQuqpSgFuBMwksF/IG8Bt3j/rUW3VVSVNSXuGszitk3vqdzF+/i/lf7mTttqL/eV+cQb/0VA7v1o6B3dpyWOdU+nZMpU+nFNJaa+BdIk97jitxSBO2o6iUBV/uZNmm3azYWsiKzbtZu63ov8ZJKuuQ0oo+nVLJ6JRC344pdG/fmu5pyXRPa023tGTaJSdosUapt4glDjN7saYT9VSVSN3s219OTt4eVmwpZNXWQtZuK+LL7XtZv6OIffsrajw3NTGebsFE0qVtEp0rf7X56vu01q2UYKRakUwc+QSWP58OfEqVVW0jsGpu2NTikObE3ckrLGH99r2s317Ehh172VSwjy0F+9hcUMzmgn3sLS2v/UIE1uLqlJpEettE0tskHfzq3DaJLsGvru2S6dIuiZRETW5saSL5VFU34GvAROBS4BVgutaVEokMM6Nru2S6tkvm2H4d/+d1d2f3vrKDiSS/sIT8PSXkF5aQVxj4d1vw38KSMrbs3seW3bUPRbZJSqBL2yS6pSXTo31reqQl071964Pf92jfWsuvyH+p9qchuKDh68DrZpZEIIG8b2a/dvf7GytAEQkwM9JatyKtdSsGdmtb43v37S9n254Stu0pZVthSfD7QFLZuruEvMJ95AUTzp6SMvaUlLHmEAP5B3RKTaRPpxQyOqXSp2MKGekp9O2USkanVDqmJjb0f1WauBr/jAgmjPEEkkYGcB/wfOTDEpH6SG4VT68OKfTqUPNqv+7O7uIythbuY3PBPjbvKmbTrmI2Fexj065A99jGXcVsLyple1EpC77c9T/X6JiaSGaXNgzo2oasLm3JCv6b3iZR4yzNVE1jHE8CQ4HXgBnuvqQxA6uJBsdFGk9FRWDsZd32wCD+uu1FrN8RGIdZt20ve0rKDnlex9REhvRox9CeaQzrmcbQnmn06tBayaQJiOTgeAVwoO1a+U1GYKvXdnW9aUPR4LhIdLk7W3bvY9XWPazeWsjqrXtYlVdIztY9FB4iobRPacURPdMY1iuNY/t1YmTfDho/iQLN41DiEGly3J2Nu4pZsrGAxRsLWLxxN0s2FrCjqPS/3hcfZwztmcbx/Tpy3GEdGZnRUasPNwIlDiUOkZjg7mwq2Mfi3AIWfLmTT9ZsZ8mm3ZRXmg0ZZzCsV3vOHNKVMwd3I7NLmyhG3HwpcShxiMSsPSVlZK/bwWdrd/Dp2h0syt3F/vKvPpMO65zKWUO6cebgrhzZq72WsG8gzSpxmNn5BJ7i6gJMdfc3a3q/EodI87K3tIw5q7fx5tKtvLNiK7v27j/4Wpe2SYwf1p1Lj+1DVteaH0eWmjX5xGFmjwHnAHnuPrRS+VjgXiAeeMTd76z0WgfgT+4+qaZrK3GINF9l5RXMXbeTN5Zu4a1lW9m466u94Uf27cDEY/swflh3rTRcB7GQOE4G9gB/P5A4zCweWEVgZnouMBeY6O7Lgq//GfiHu8+v6dpKHCItg7uzeGMBz8zdwAsLNx18BLhdcgIXDu/FJcf25vBuUX/QM2Y0+cQBYGYZwMuVEsco4Ffuflbw+JbgW+8Mfr3l7m9Xc63JwGSAPn36jFi/fn2EoxeRpqSopIyXF23in59t4PMNuw6Wn9C/E1NOz2TUYZ00V6QWsZo4LgbGuvt3g8dXAMcRaIVcSaAFstDdH6rpumpxiLRsSzcVMOOzDfx7wcaDrZCRfTsw5fRMThnQWQmkGrGaOL4OnFUlcRzr7teGeD3NHBeRgwqK9/PkR+t47MO1BwfUh/VKY8ppmXxtcFclkCrqmziitRlyLtC70nEvYFOUYhGRGJfWuhXXjcnig5tO5+Zxh5PeJpFFuQVMfmoe4+6dw7srttKUniCNddFqcSQQ6JYaA2wk0DV1abhLtqurSkQOpbi0nOmffcnfZn/B1t0lAIzOTOfW8YMY1F2D6E2+xWFm04GPgYFmlmtmk9y9DJhCYP/y5cDMcJKGmU0ws2kFBQWRCVpEYlrrxHi+M7ofs356Gr8YP4h2yQl8kLON8ffN4eZ/LSKvsPZ9SqR6TWoCYKg0xiEi4dhZVMq976zm6U/WU1bhpCbG84PTMpk0ul+LnAcSE4PjkaKuKhEJxxf5e/j9qyt4e/lWAHq2b80vzxnEWUO6tagB9CbfVSUi0lT079yGR64cyT+/exyDurdj465irnl6Pt95Yi4bduyNdngxIyYTh8Y4RKQ+TshM5+VrR/Ob84bQNimB91bmc8bds5j6Xg6lZRXRDq/JU1eViLRoeYX7uOOV5bywMDAjoH/nVH57/hGM6t8pypFFTovsqlKLQ0QaSpe2ydx7ydH847vHcVh6Kl/kFzHx4U+44ZmFbN9TEu3wmiS1OEREgkrKyvnbrDX8Ndhl1SGlFb8YP5gLh/dsVoPnLbLFISISCUkJ8Vw3Jos3f3wyJ2Z2Yufe/dz47Odc8ehnrN9eFO3wmgwlDhGRKjLSU3l60nH8+etH0j6lFR/kbOOsv8zmoVlfsL9cg+cxmTg0xiEikWZmXDSiF+/ccArnH9WDffsruPO1FZz71w9ZlLsr2uFFlcY4RERCMGtVPrf+ezG5O4uJM/jWqAx+ctZA2iQlRDu0sGmMQ0SkEZwyoDNvXn8yV5/UDzPjiY/WccafZ/HG0i3RDq3RKXGIiIQoJTGBW8cP5sUpJ3Jk7/Zs2b2P7z01j+8+mf1fe6I3dzGZODTGISLRNKRHGs9//wR+fd4Q2iQl8PbyrXzt7lk8MmcNZS1g8FxjHCIi9bClYB+/fnkpry4OdFkN7t6OX583hJEZHaMcWfU0xiEiEkXd0pJ54LIRPHrlSHq2b82yzbu5+KGPueGZheTtbp77fihxiIg0gDGDuvL2Dadw3ZgsEhPieH7BRk7/8ywenr2m2c39UOIQEWkgrRPjueFrA3j7+lM4Y1BX9pSUcceryxl37xw+WL0t2uE1mCaTOMzsMDN71Myei3YsIiL10adTCo9cOZLHv30M/dJTycnbw+WPfsrVf88mJ68w2uHVW0QTh5k9ZmZ5ZrakSvlYM1tpZjlmdjOAu69x90mRjEdEpDGdNrALr//4JH42diApifG8tWwrZ94zm1ueXxzT4x+RbnE8AYytXGBm8cBUYBwwGJhoZoMjHIeISFQkJcTzg1Mzef+np3LZcX0wM6Z/9iWn/PF97n5zJXtKyqIdYtgimjjcfTawo0rxsUBOsIVRCswAzotkHCIi0dalbTJ3XHAEb15/MmcN6Urx/nLuezeHU/7wHk9+tC6mdh6MxhhHT2BDpeNcoKeZdTKzh4CjzeyW6k42s8lmlm1m2fn5+ZGOVUSkQfXv3Ia/XTGS564ZxfA+7dleVMptLy7l1D8GEsi+/eXRDrFW0Ugch9oNxd19u7tf4+793f331Z3s7tOA24H5iYmJEQtSRCSSRmZ05F/fP4GHLh9BVpc2bCrYx20vLmX0Xe8xbfYXFDXhLqxoJI5coHel417ApijEISISVWbG2KHdeOPHJ/PQ5cMZ0qMd2/aU8LtXV3DiXe9y/zurKSjeH+0w/0fElxwxswzgZXcfGjxOAFYBY4CNwFzgUndfGu61teSIiDQn7s77K/O5/93VzP9yFwBtkhK4eEQvLj++D5ld2jbIfeq75EhEE4eZTQdOBdKBrcBt7v6omZ0N/AWIBx5z9zvCvO4EYEJmZubVq1evbtigRUSizN35eM127n8nh4/XbD9YfvxhHbni+AzOHNKVVvF17zBq0okj0tTiEJHmbummAp7+5EteWLiRvaWBgfMubZO45JjeXHJsH3q0bx32NVtk4lCLQ0Ramt379vPv+Rt5+pP1rM7bA4AZHJvRkQlH9mDc0G50apMU0rVaZOI4QC0OEWlp3J1P1+7g6U/W8+ayrQfnf8THGSf078SEYT04a0g30lJaVXuNFpk41OIQEYHCfft5a9lWXvp8E3NWb6OsIvB53ireOKF/OqcO7MwpAzrTLz0Vs69mQrTIxHGAWhwiIgE7i0p5Y+kWXlq0iY+/2E5FpY/23h1bc8qAzpw6oAuj+neiTXKrlpc41OIQEanetj0lzFqZz6xV+cxenc+uvV/NBWkVb+T8bnzLSxwHqMUhIlKz8gpn8cYC3l+Zx6xV+Xy+YRdr7zxHiUNEREKzs6iUjm2StOe4iIiEpkNq/df4i8nEYWYTzGxaQUFBtEMREWlxYjJxuPtL7j45LS0t2qGIiLQ4MZk4REQkepQ4REQkLEocIiISlphMHBocFxGJnphMHBocFxGJnpieAGhmhcDKCN8mDahr0ybUc0N5X3XvCae8alnV43RgWy1x1FdTqM+aXg+l3kIpa4y6rC6Ohj4vGvXZXH82Q3lvY/yuD3T3um8n6O4x+wVkN8I9pkX63FDeV917wimvWnaI4xZRnzW9Hkq9hVLWGHVZn/oM57xo1Gdz/dmsT302pd/1mOyqamQvNcK5obyvuveEU161rD7/t7pqCvVZ0+uh1Fs4ZZFW13uGc1406jOW6jLcc+tan03mdz3Wu6qyvR7rrch/U302HNVlw1J9Nqz61mestzimRTuAZkb12XBUlw1L9dmw6lWfMd3iEBGRxhfrLQ4REWlkShwiIhIWJQ4REQlLs00cZna+mT1sZi+Y2ZnRjieWmdlhZvaomT0X7VhilZmlmtmTwZ/Jy6IdT6zTz2TDCvfzskkmDjN7zMzyzGxJlfKxZrbSzHLM7OaaruHu/3H3q4GrgG9GMNwmrYHqco27T4pspLEnzLq9EHgu+DN5bqMHGwPCqU/9TNYuzPoM6/OySSYO4AlgbOUCM4sHpgLjgMHARDMbbGZHmNnLVb66VDr1F8HzWqonaLi6lP/2BCHWLdAL2BB8W3kjxhhLniD0+pTaPUH49RnS52VCw8XYcNx9tpllVCk+Fshx9zUAZjYDOM/dfw+cU/UaZmbAncBr7j4/wiE3WQ1Rl3Jo4dQtkEsgeSyk6f7BFlVh1ueyRg4v5oRTn2a2nDA+L2PpB7gnX/3FBoFfxJ41vP9a4AzgYjO7JpKBxaCw6tLMOpnZQ8DRZnZLpIOLcdXV7fPARWb2INFZTiNWHbI+9TNZZ9X9fIb1edkkWxzVsEOUVTt70d3vA+6LXDgxLdy63A4o+YbmkHXr7kXAtxs7mGaguvrUz2TdVFefYX1exlKLIxfoXem4F7ApSrHEOtVl5KhuG5bqs2E1SH3GUuKYC2SZWT8zSwQuAV6MckyxSnUZOarbhqX6bFgNUp9NMnGY2XTgY2CgmeWa2SR3LwOmAG8Ay4GZ7r40mnHGAtVl5KhuG5bqs2FFsj61yKGIiISlSbY4RESk6VLiEBGRsChxiIhIWJQ4REQkLEocIiISFiUOEREJixKHNGtmVm5mCyt91biEfGOxgHfNrG+l2LaY2cZKx4lVzrkq+Gx+5bJ0M8s3syQzm2FmWY37P5GWSPM4pFkzsz3u3qaBr5kQnEhVn2uMB85w9+srlf0K2OPuf6rmnHbAGqCPu+8Nll0DHOPuk8zsFODy4L4KIhGjFoe0SGa2zsxuN7P5ZrbYzA4PlqcGN8CZa2YLzOy8YPlVZvasmb0EvGlmKWY208wWmdkzZvapmY00s0lmdk+l+1xtZncfIoTLgBdqiG+Emc0ys3lm9oaZdXf33cBsYEKlt14CHGiFzAHOMLNYWrxUYpAShzR3rat0VVXe3Wybuw8HHgR+Eiy7FXjX3Y8BTgP+aGapwddGAVe6++nAD4Cd7j4M+A0wIvieGcC5ZtYqePxt4PFDxHUiMO9QAQfPvR+42N1HAI8BdwRfnk4gWWBmPYABwHsA7l4B5ABHhlAvInWmv0ykuSt296Oqee354L/zCGztCnAmgQ/+A4kkGegT/P4td98R/H40cC+Auy8xs0XB74vM7F3gnODmOK3cffEh7t3R3QuriWsgMBR4K7AfGfHA5uBrLwMPBLutvkFgO9rKOwrmAT2oJimJNAQlDmnJSoL/lvPV74IBF7n7yspvNLPjgKLKRTVc9xHg58AKDt3aACgzs7hgK6EqA5a6+6iqL7h7sZm9DlxAoOVxfZW3JAPFNcQmUm/qqhL5b28A11rwT30zO7qa931A4C9+LLBn8xEHXnD3TwnseXApX40/VLUSOKyG1zqb2ajg9VuZ2ZBKr08HbgC6Ap9UOXcAoNVjJaKUOKS5qzrGcWct7/8N0ApYZGZLgseH8gCBD/dFwE3AIqCg0uszgQ/dfWc1578CnHqoF9y9FLgYuMvMPiewT/kJld7yJoHuqGe80mORZtaVQNfcZkQiSI/jitSBmcUTGL/YZ2b9gXeAAcEPfczsZeAed3+nmvO7A3939681YEzXA7vd/dGGuqbIoWiMQ6RuUoD3gk9AGfB9dy81s/bAZ8Dn1SUNAHffbGYPm1m74GO2DWEX8FQDXUukWmpxiIhIWDTGISIiYVHiEBGRsChxiIhIWJQ4REQkLEocIiISFiUOEREJy/8DZIEjGfZ+y0cAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.loglog(ETeV,gam, lw = 2)\n",
    "plt.gca().set_xlim(1e-2,1e2)\n",
    "plt.gca().set_xlabel('Energy (TeV)')\n",
    "plt.gca().set_ylabel('Mean free path (Mpc)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate and plot optical depth w/ and w/o LIV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "from time import time  as t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "it took  0.4837489128112793 s\n"
     ]
    }
   ],
   "source": [
    "t0 = t()\n",
    "# with LIV\n",
    "tauLIV = ebl.optical_depth(z0,ETeV, LIV_scale=1e-7, nLIV=2)\n",
    "#w/o LIV\n",
    "tauCalc = ebl.optical_depth(z0,ETeV, LIV_scale=0.)\n",
    "t1 = t()\n",
    "print('it took ', t1 - t0, 's')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Do the plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'Attenuation')"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEKCAYAAAARnO4WAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAvm0lEQVR4nO3deXxU9b3/8dd3JpMMWSEhCUsIIQFkDQkEFRBkUUQQxF3EBWtd7r3W2l9tq7YutVZtr71tqdaliqhVUVGsdamKShG1SiKL7AQIEEIgC4Tsycx8f3+cSUgg22SZM8vn+XjkkZxlzvnMl+Gdk+8553uU1hohhBCBz2J2AUIIIbxDAl8IIYKEBL4QQgQJCXwhhAgSEvhCCBEkJPCFECJIhJhdQEuUUvOB+VFRUTcPHz7c7HKEEMJv5OTkFGut41tapnz5OvysrCydnZ1tdhlCCOE3lFI5WuuslpZJl44QQgQJCXwhhAgSEvhCCBEkJPCFECJISOALIUSQkMAXQoggIYEvhBBBQgJfCCGChAS+EEIECQl8IYQIEhL4QggRJCTwhRAiSEjgCyFEkJDAF0KIICGBL4QQQUICXwghgoQEvhBCBAkJfCGECBIS+EIIESQk8IUQIkhI4AshRJCQwBdCiCAhgS+EEEEixFs7UkpFAH8F6oA1WutXvLVvIYQQXTzCV0otU0odVUptOWX+HKXUTqVUrlLqbvfsS4GVWuubgQVd2a8QQgjPdfUIfznwBPBSwwyllBV4EjgfyAfWK6XeBZKA792rObu431YVFx4k9n9Ho/6vvNV19FXhcIbNmFhbg1pT2/KKEQr90+jGSfWnE3BCt7zu1DD0DLux/d31WFZUtbp/5x3REGNBo7C+UYHaVe/eQZOVFLjSQqm9KgaNgkpN+BPFxvYVoJSxvgKtFCcu7kvt0EhcStHr2xNEfnPcWM9iQVtAWyxoi8IRYePA9WNxKSsuSwhJK7djqXXhCgnBFWLFZbOhbTZcNhtlY9I4njkKZbMTVl5H7y37ICICFd0HFdsXW98EbPH9CEvoT8SQoVjDwlp9z0II83Up8LXWa5VSKafMPhPI1VrvBVBKrQAuxgj/JGAjbfxloZS6BbgFIDk5uVN1WVQrodywDzS0s87Jepqs19ZLlD65rtJtrhticZ0Md63B1fJ6VpeLcNXwy8gFtcZGVQvF9HGeAKqN2ZU1UFLX8kYjFaPrNp2c3lEOFS0XO6BuO8SsNiZ21cNr1a2/qTsjqYwO54SKotfHJ7AV1FIfEU59ZCT1vXvj6DcAnTQYS/p4ImacT5++/VEWOYUkhDcprTsWfK1uwAj897TWY9zTlwNztNY/dE9fB5wF/ALjr4EaYF1H+vCzsrJ0dna2R/VolwvtchlB2nrRxhcY67W1btNQcrWSzE3W1Q3ba2Fd7XL/YWO1nozr+np3zU6jbvd+tMuFVqDDwnC5XGinE33iBDgdaIcD7XLhcjrA4UA7nTiiI9GhNmNecTGUlKDr640vRx26zvhy4aJq+BBcTgcuZz2R/8lBVVRAbQ2qrgZVWwt1tajaWipS4ikfEoulvpqIvAL6/Xs71ro6rLV1hNTUEVLrwFLrwFLrhNujsPRyv6eXKmFfK3/EnRECV4dTp0MorY6h90tHqEqIpWbgQBypw7CNziD67GlEjM0Am63t9hZCnEYplaO1zmppWU+ctFUtzNNa60rgxh7YX/OdWyw9d+RotXb/NkM96AaJje/Yeqke7H/W1R6s3DqX08mJ8uOUlx6lLvM7nHl74GgBqugIISVHCSsuwn7sOK4BVkJQRKtK+h07CkeqsB+pgu/zgW+AvwOgLbDn5gkUjZ+Erf8Y+lniSEw+A9uoUT3z7yBEEOiJwM8HBjWZTgIKemA/wodYrFaie8cR3TsOUke2u351ZTkle7dRnbUGtm8iLC+X8MMFRBQfx15SjSpzMTR8J0MLdhufnreqYIsDV6iFiuR4KseMxTbzAuIWXIFKTj75F5sQolU90aUTAuwCZgGHgPXANVrrrR5scz4wf+jQoTfv3r27S/UJ/+N0ODi8NZujh7ZRe3gr9pLtDFv5FRE7T6BaOGl+LD2RHXfeQJ9Rs0hNP4fQ0DD5K0AErba6dLoU+Eqp14DpQF/gCPCA1vp5pdRc4E+AFVimtf5tZ7bfmT58Ebgqy4+Tv/ZD9Kf/IGZzDnF5+YQeqoGzQmG6cYVUzUELtlcrODFiCI5zZxF79Q1YJ06UXwAiaPRY4Pc0CXzRFu1yceRgLgXffUr9ofUkln5Lyn/2wIc1zdZzhNsoyxxD6GWLiLr9TjkZLAKaBL4IGsWFByj44FUiPn+Xflu3ErGvHI67P+Phityfj6NowLn0GTePYToaa2YmhIaaW7QQ3UgCXwStQ3u2Ufzu88SteZ/E6nxsk91dOzUa/b/luGwhlGeNo9fiGwm75DLo18/cgoXoIr8LfDlpK3pCTVUFu775F1Xb/sWQrZ+Q+FY+HG1+v0RFahIsvJTIX/wSEhJMqlSIzvO7wG8gR/iip2iXiwM7N1D84XIGrnmPxO35qDwHOIzle382liMj5pCQtZBUVzgqLQ2ioswtWogOkMAXoh0lR/LZ9/nr9PnoNQbt307oNPctKlrjWFqNpcJF5bgx9LpiESHzF8DIkXLtv/BJEvhCeKCmupKdX/+Tmi3vM7RgDXGvHoaDzYeKqI2LwXn+bMLv+RWkp5tUqRCnk8AXopNcTie5m76gbM2rJH/9MYm7CmCPAyqN/zeFS1LIO+tCws+YybDaXvTSFjjnHOn+EaaRwBeimxw9tI+8r1YRtW4VqVs2EjZZgdXo2nG9Wo1ldz3aoqgeOZzQCy8i5NzpcNZZEN/BcZCE6CK/C3y5Skf4g7raGnat/5jybZ8Qd/Q/DP30eyy766HAddrw2Ecvmonrr8+QODAV5XAYo6na7eYULgKa3wV+AznCF/7kxPES9uZ8Qv2mj0jasJb++/Ih3wkFTpgZBmeHcZRYyvP6kvpyDjXJA1GZEwg7axIqIwPGjZP7AESXSeALYYKyY8Uc2LyWil1fEln4HYOdO4mmEr6tg3/VtPiQnPre0Rz49H2SRo0nzB4OmzbBgAHSJRRM9n8NOS/AgicgxPO7wCXwhfABLqeTg7nfc3T7OnTed/TdmUP//Dx6HamGIy4odEKEgjuicGgLh6wD6P/4HkLLa6mPicI5fDi2zAlYx46FUaMgIwNiY81+W6K7HPwWPn8E9n5uTM9fChNu8Hgz3n4AihCiBRarlcFnZDD4jIzGedrlouDAbo7sWk/Nge8Iz99CrDrIAH2YwXUHIdoBtWArK8e2PgfW5zS+dt/VMzm2ZAn9hk+kf3kdasMG45fAyJEyPpA/yc+BNY9ArvtxoqFRMOm/YdTF3b4rOcIXwgfVVFdyaPcmjuVtor7ge2L2byExP5e4o6VQ5IJil3FeIMU4Zqv5twv7mgoAXCFW6tNSsU2bjuWcc2DKFEhLM/PtiJYUboHPfgO7/mVMh0bCWbfBpP+B8M7/5eZ3XTpylY4QLassP86h3Zs4vn8zriPbCT++g4E1ucRtL4Yt9VDogtLm4wPV9I1hzzuvkJoxjV4RUbBrFwwbJncKm+m7l+D9n4KzDmzhcOYtMPkOiIjr8qb9LvAbyBG+EO3TLhclhQc5tPNbqvdvwH5oI0l7vqfvwSLjDuEYC1xop15b2V83mLTfbcQR35eQBQtRF10Es2ZBZKTZbyM4OOrgo3tg/XPG9IQlMONXENl9J+Ul8IUIQkcP7ePgps9w7PuKuNINDHHsxXqoHt6ohvKT/+9dthD0tHOxXnklXHONhH9PqTgKb1wPB74Gayhc9EfIvLbbdyOBL4SgvKyUvTmrqdnyHkO3fkrc7hLY7TDuFQC0gp0rn2P4xUuwyCMhu1d+Drx+LZQXQNQAuOrvkDShR3YlgS+EaEa7XOz5/muKcv5B4q5/kbpjt9H3P9NOgUrkQPJljH/pK0IvnAc33ACJiWaX7L82vgr/vBOctZA8Ca54EaJ6rj0l8IUQbSo8mMu+1X9jyP636EcR5DngxSoAXDYb6rbbUPfcA/37m1ypn9m7Bl5yX16ZdRPMeaxTN1N5QgJfCNEhToeDLV+8g/ObZaRvWkvIhlrYaTwVxhUWhuX22+EXv5A7fzuirgqemgzH9sG0n8HMX3llt34X+HJZphDmKz58gN0fPcWob5YRs6YUdhjBXz+wP7b9B0H6+dv2yf3w5Z8hYRTc8u8eP7Jv0FbgW7xSgYe01v/UWt8SExNjdilCBK2+/ZOZtORR1GPb+equO6m6ORaGh2AbVcrGP11Cwe5NoLXxJZo7vAm+egJQsOAvXgv79vjkEX4D6dIRwneUHC1g58qHmFD4JmEWhzHez54skutjUc88A717m12ib3A64LlZcHijcefshb/z6u797ghfCOF74hIGMPm/n6b0h9/wn97zCKl3MvjtT1FvvIFzzBhYt87sEn3DN08ZYR8zyGv99h0lgS+E8Ej/5KGcfeer5Mx6nmM3JcIAC9ZDh9DnToMHHgCHw+wSzXMszxjxEmDe/0GYbz3qUgJfCNEpE2ZdQc3d/2HzrVNgSijKpeGhh9BTz4H8fLPL8z6t4b2fQH0VjLkMhs82u6LTSOALITqt/8DBjLr3M9b84C4c10VClEL95xvqFl1pdmnet/l12PMZ9OoDc7zbb99REvhCiC4JCQlh+g9+y5b/eoeC24bAyBCOZx2mvOiA2aV5T2UJ/Ose4+fZv+3WwdC6kwS+EKJbZEyahe3nX7Pj6nEkxJRy/Jl51JQVgdNpdmk9b9OrUF0KKVMh4xqzq2mVBL4QotvE900g6qZ32KuSGFS/n6orM3FdOAfq6swurWftdD/EJOtGn37OgE8GvlJqvlLq2bKyMrNLEUJ4aODAJPTiVRRU9SV2XQGWT1ajF10duEf6VaXGkMeWEBh6ntnVtMknA1/utBXCv6UNHU7x9W9TurgfhIF6exXc/MPAvCs3dzVoJwyeAnbfziyfDHwhhP9LH5vBzhtfp3xRHIQALyyH++4zu6zut/ND4/sZF5pbRwdI4AshesykSVP55pqXqL4yBhToRx+BjRvNLqv7OOog91Pj5+FzzK2lAyTwhRA96rzz5/LJFU/hmmjcnOW4/TazS+o+B76C2jKIHwmxQ8yupl0S+EKIHjf/4qv5ZP7lMMFG0bk1xgBjgaDh6pwzfP/oHiTwhRBeoJRiyA1/oGBeEv1t+yhc/RezS+o6rWHnB8bPw32//x4k8IUQXjJ8UD/+nfZzAGK+fAz9+nJzC+qqoh1wfD+Ex0FSi6MR+xwJfCGE18y94iY+0xPp9WIJ6uob4aOPzC6p8xquzhk+Byz+8fQvCXwhhNfE9LJRPvMRaofZAdA/XALV1eYW1VlNA99PSOALIbzqoqkTeXnWzZBgQeUXwq/vN7skz1UUQf56sIZC2kyzq+kwnwx8GVpBiMBltSgyr76HPXONyxj1H/4Ptm0zuSoP7f4Y0MZgaWGRZlfTYT4Z+DK0ghCBbcKQBP5x7sO4JoSiHC64eYnZJXmm4eocP7i7timfDHwhROBbfOlC3pw5D+zAV+th0yazS+qY+hrY87nxsx/134MxwoUQQnhdYrSdigseoHLTZ0SE1aJtxfjuwMJN5K2D+kpIHAu9B5ldjUfkCF8IYZprzx3DG1MugfGhlK5fbnY5HdPYneNfR/cggS+EMFFYiJXK9BtwakXvfe9DeaHZJbVNa9jlvnfAz/rvQQJfCGGyCyZN4N9lY7F+WInrsrlml9O2wu/hRD5EJkL/TLOr8ZgEvhDCVMMSo/g69iLIrkOt3gAH8swuqXUNR/fDLwCL/8Wn/1UshAg4SbOupHx4FEoDf/yV2eW07lCO8X3IuebW0UkS+EII0y3IGMg/xhl3rOrXVvnuoxALvze+9x9nbh2dJIEvhDBdn4hQcmbfiSvSgjpSBe++aHZJp6sqNfrvbeEQm2p2NZ0igS+E8Alzzx7J92OHA6Cf+L3J1bSg4eg+cbTfjI55Kgl8IYRPmH5GPM9kLTEm1u6AooOm1nOahsDvN9bcOrpAAl8I4RNsVguJ51/EnpkpqOvDYcdKs0tqrvEIf4y5dXSBBL4QwmdcPiGJ35x1Ewy0orOfB5fT7JJOajzCTze3ji6QwBdC+IzRA2I4En8Oea5EVNlB2PGB2SUZHLVQvBNQkDjK7Go6TQJfCOFTLstK5vO8MbC8Eu7+qdnlGIp2gMsBcUMhNMLsajrNJwNfHoAiRPC6OGMga3UG7HfCmlw4stPskgLihC34aODLA1CECF7xUWGomXOo7N0LyjW887TZJUngCyFET7lswiC+GzLSmPj4fXOLAQl8IYToKbNGJvDJkOkA6M37oOKoecVoDYVbjJ8l8IUQonvZbVYcM2YDoA44YOt75hVz/ADUlkFEvDEssh+TwBdC+KRhmcMp7tMb6oCPXjGvkKbdOcovHsLYKgl8IYRPmpgSy1/PvgLXgl5QtRFqy80pJED670ECXwjho0b2j+bNiZeRM2402B2Qu9qcQgLgDtsGEvhCCJ9ktSjGD+7Dx84sY8YOk67WkSN8IYToeRNT+lC6pxd8WANfvA+OOu8WUH0Myg5AiB1i07y77x4ggS+E8FkTU2KZtHUzfFsHO45B3hfeLeDIVuN7wiiwhnh33z1AAl8I4bPGDepN9mB3V8p+p/e7dQKoOwck8IUQPsxus3LirMnGxH4nbH8fXC7vFSCBL4QQ3pM8YQwFUX2hWsPeQ1CwwXs7L9xsfA+AK3RAAl8I4ePOTI3lm0Hup0ztd8IOL91166iDozvw9zHwm5LAF0L4tAnJsXyb3BD4Du/14xfvBFc9xKZCWJR39tnDJPCFED4tJtzGkfGT2JqYirN/uBHExbt7fseN/ff++wzbU0ngCyF83sCJ6cxbspTdl80zZnijWydARshsSgJfCOHzJg6JBeBj5wRjhje6dQLshC1I4Ash/MDElD6gNWt39EIfsED+eigv7Lkdah1wl2SCBL4Qwg/0j+nFFEcRK5/8H/TbNUYg7/yw53ZYlg81xyE8DqL699x+vEwCXwjhFxInZlAcHoOlrAZKXbDn057bWQCNgd+UBL4Qwi9MTI3j26TRxsR+J+xdC05Hz+wsALtzQAJfCOEnJqb04ZtkI4D14XDjsYOHsntmZ0fcgZ8ogS+EEF6XFh/JjuGZAOi8eqMfP7eHunUKNhrf+wfOFTrgxcBXSqUqpZ5XSq301j6FEIFDKUXviZkcs0dhKS6H47pnnoJVXghlByE0CvoO7/7tm6hDga+UWqaUOqqU2nLK/DlKqZ1KqVyl1N1tbUNrvVdrfVNXihVCBLeJaX1ZP2g0jhAbFCtjILXKku7dSb67m2jgeLBYu3fbJuvoEf5yYE7TGUopK/AkcCEwCliklBqllBqrlHrvlK+Ebq1aCBGUslJiuWvunVzw8IcwezqgYe/n3buThvMCSVndu10f0KHA11qvBUpPmX0mkOs+cq8DVgAXa62/11pfdMrX0W6uWwgRhEYPiKY+KoY9x2qoGHSuMXPPZ927k8Yj/CAN/FYMBA42mc53z2uRUipOKfU0kKmUuqeN9W5RSmUrpbKLioq6UJ4QItDYrBbGD+4NwEbbeGNm7qfGCdzu4HKeHG8/WI/wW9HS3QittrrWukRrfZvWOk1r/Wgb6z2rtc7SWmfFx8d3oTwhRCAa0y+KN175ORPmXQphCVBRePLZs11VtAPqKqB3MkQGXk90VwI/HxjUZDoJKOhaOUII0ba0xGjiqsroVX4c7MZlmt12123+euN7AHbnQNcCfz0wTCk1RCkVClwNvNs9ZQkhRMtS4yPIjXMfa9a7e5G763r8hv77pIndsz0f09HLMl8DvgbOUErlK6Vu0lo7gNuBj4DtwBta6275u0opNV8p9WxZWVl3bE4IEUDS4iMbA1+XWgEFB76Gusqub/xQjvE9APvvoeNX6SzSWvfXWtu01kla6+fd8z/QWg9398v/truK0lr/U2t9S0xMTHdtUggRIPpEhFIwYAgAtTv2woBMcNZB3rqubbjmBBzdDhZbQI2B35QMrSCE8Dt1w0cA4Nq6DYbOMmZ2tVunYAOgjQHTbPaubctHSeALIfxO6JiRuFCE5e2F5Ibr8bsY+AF8w1WDELMLEEIITyUnxbN0ytWkpQ9jfv9MCIuBklw4th/6DO7cRgP4hqsGPnmELydthRBtSY2P5E/nLOaNCXMhIhJSpxkLOnuUr3WTK3Qk8L1KTtoKIdqSFh8BwN4i95U5aV3sxz9+ACqPQq9YiE3thgp9k3TpCCH8zqDYcOJrTjDxy8+pe+UIofPcgb/33+CsB6vNsw027b8PoEcansonj/CFEKItNquFs1zH+NN7f8D56KPGUAhxw6Cu/OTdsp7Id19/H8D99yCBL4TwU5ZRowAIzd0NTufJyzN3feT5xhp+SSRN6KbqfJNPBr6ctBVCtGfA4H4cjozDWlsD+/fDyPnGgk0rjG6djnLUweFNxs8DJfC9Tk7aCiHa02xMnW3bYPAUo1unohB2/avjGzqyBZy1xmt79emZYn2ETwa+EEK0Jy0+kty+TQJfKci60ZjOXtbxDQXB5ZgNJPCFEH4prckRvt623Zg5bhFYw4ynYJXu69iGguAO2wYS+EIIv9Q7PJQjSanUWm1UV9UYM8NjYfQlxs/fvdixDQXBHbYNJPCFEH6rfMKZjPx/K/nut385ObOhW2fD340Tsm2pKoXSPRBih8TRPVeoj/DJwJerdIQQHTEkMQaXxcqeooqTMwedBQmjoLIIdrzX9gYaxr8fkOn5zVp+yCcDX67SEUJ0RFp8JAB7j5ZDvftSTKVggvsoP+eFtjfQ2J0T2JdjNvDJwBdCiI5IjY/gtv+s5BfXTIInnzy5IP1KCOkF+9ZCcW7rG2i84Srw++9BAl8I4cfS4iOptoURXl1pXJrZoFdvGHuZ8XNrR/lHtsLBb42fA/QZtqeSwBdC+K2kPr3IizcuzXQ2DXyACT8wvm98Feprmi/LXQ3PX2CMvZM2E6IHeqFa80ngCyH8VojVQu0w43GHestWY1z7BgPHG8+mrS6F7e+enJ/9ArxypRH2oy+Bq18L6BEym5LAF0L4tZi0ZE6ERRBSdhyOHj25oNmdty+AywWf3A/v3QnaCef8P7hsWcA+v7YlEvhCCL+WlhBFblySMXFqt87YKyA0Eg58BX+/BL78MygrzF8K5z0AluCKQJ98t3IdvhCio1LjI9kdl2xMbN/efGFYlBH6AHvXQFg0XLsSJtzg1Rp9hU8GvlyHL4ToqLT4CFaNmcGfr/o5zJ59+gpn3mzcSRszCG762DhJG6TkEYdCCL+WGh/Jf5LT2Wiz8KPUtNOPYhNHw49yIDwObL3MKNFn+OQRvhBCdFRMLxt9I8OoqXdRUFbdykpJQR/2IIEvhAgAqfERXLR9La5f3A3Hj5tdjs+SwBdC+L20+Ehu+fZtkv/2l9Ov1BGNJPCFEH4v7dTHHYoWSeALIfxeWnykBH4HSOALIfxeanwE+TGJxsTBg+YW48N8MvDlxishhCeS+oRTHNMXAKcEfqt8MvDlxishhCesFkXIYPeomQcPmVyN75Ibr4QQASEmdTDH7ZEQFUOo1kEzAqYnfPIIXwghPDV4YBwZP17B80vfkrBvhd8d4dfX15Ofn09NTU37KwvRDex2O0lJSdhsgf+Qa3+WHBcOwKFjrdxtK/wv8PPz84mKiiIlJQUlv8VFD9NaU1JSQn5+PkOGDDG7HNGGxGhjXPvCsmpwOsFqNbki3+N3XTo1NTXExcVJ2AuvUEoRFxcnf1H6gX7Rdn669mWW3T4DnnnG7HJ8kt8FPiBhL7xKPm/+oV+0nXprCPb6WjgkV+q0xC8DXwghThXdK4QS97X49fsPmFyNb5LA74T8/Hwuvvhihg0bRlpaGj/+8Y+pq6sztabIyMhOvW7NmjVcdNFFba6zceNGPvjgg8bpd999l8cee6xT+xOipyilqOvXH4D6A/kmV+ObJPA9pLXm0ksvZeHChezevZtdu3ZRUVHBL3/5S7NL6zGnBv6CBQu4++67TaxIiJa5+g80fjgkgd8SCXwPffbZZ9jtdm688UYArFYrf/zjH1m2bBlVVVXMnTuXzZs3A5CZmclDDz0EwH333cdzzz3HmjVrmD59OpdffjkjRoxg8eLFaK1P209ubi7nnXce48aNY/z48ezZs4eKigpmzZrF+PHjGTt2LP/4xz9arPH3v/89Y8eOZdy4cY3BPH36dLKzswEoLi4mJSXltNd9++23TJ48mczMTCZPnszOnTupq6vj/vvv5/XXXycjI4PXX3+d5cuXc/vttwOwf/9+Zs2aRXp6OrNmzeLAAeNP6SVLlnDHHXcwefJkUlNTWblyZRdaXYiOsSYbd9vaCg+bXIlv8snLMpVS84H5Q4cObXO9lLvf75H95z02r9VlW7duZcKECc3mRUdHk5ycTG5uLtOmTeOLL74gJSWFkJAQvvzySwDWrVvHtddey+HDh9mwYQNbt25lwIABTJkyhS+//JJzzjmn2TYXL17M3XffzSWXXEJNTQ0ul4vQ0FBWrVpFdHQ0xcXFnH322SxYsKDZScUPP/yQd955h2+++Ybw8HBKS0s7/L5HjBjB2rVrCQkJYfXq1dx777289dZbPPTQQ2RnZ/PEE08AsHz58sbX3H777Vx//fXccMMNLFu2jDvuuIN33nkHgMOHD7Nu3Tp27NjBggULuPzyyztcixCdEd2vL9UhYfSqqoQTJyA62uySfIpPHuH78lg6WusWr9pomD916lTWrl3LunXrmDdvHhUVFVRVVZGXl8cZZ5wBwJlnnklSUhIWi4WMjAzy8vKabau8vJxDhw5xySWXAMaNP+Hh4Wituffee0lPT+e8887j0KFDHDlypNlrV69ezY033kh4uHETSmxsbIffW1lZGVdccQVjxozhJz/5CVu3bm33NV9//TXXXHMNANdddx3r1q1rXLZw4UIsFgujRo06rU4hekK/3r14dPoS3vnRb+Q6/Bb45BF+R7V1JN5TRo8ezVtvvdVs3okTJzh48CBpaWmEhISQnZ1Namoq559/PsXFxfztb39r9ldBWFhY489WqxWHw9Fsey118QC88sorFBUVkZOTg81mIyUl5bTrw1v7hRQSEoLL5QJo9Zry++67jxkzZrBq1Sry8vKYPn166w3Riqb7bvo+W3tPQnSnxGg7D0+YT9GYfiyMiDC7HJ/jk0f4vmzWrFlUVVXx0ksvAeB0OvnpT3/KkiVLCA8PJzQ0lEGDBvHGG29w9tlnM3XqVB5//HGmTp3a4X1ER0eTlJTU2DVSW1tLVVUVZWVlJCQkYLPZ+Pzzz9m/f/9pr509e3bj+QSgsUsnJSWFnJwcgFb708vKyhg40Djp1bTbJioqivLy8hZfM3nyZFasWAEYv5BO7ZoSwpsa7rY9ckJulGuJBL6HlFKsWrWKN998k2HDhjF8+HDsdjuPPPJI4zpTp04lMTGR8PBwpk6dSn5+vkeBD/Dyyy+zdOlS0tPTmTx5MoWFhSxevJjs7GyysrJ45ZVXGDFixGmvmzNnDgsWLCArK4uMjAwef/xxAO666y6eeuopJk+eTHFxcYv7/PnPf84999zDlClTcDqdjfNnzJjBtm3bGk/aNrV06VJeeOEF0tPTefnll/nzn//s0fsUojv1i7aTVnyQSZ+shDVrzC7H5yhf/lM7KytLN1xZ0mD79u2MHDnSpIpEsJLPnX+oqXfy6Nz/4dern0Hfeivq6afNLsnrlFI5WuuslpbJEb4QImDYbVbK4xIAqN8vT746lQS+ECKgONw3Xznz5earU0ngCyECikoyAt96uMDkSnyPBL4QIqDYkwbgUBZCS4qhttbscnyKBL4QIqAk9ImgKKKPMXFYhlhoSgJfCBFQEqPtFEb1pToiCjwYWiQYSOB3ggyP7L3hkefOncvx48c5fvw4f/3rXxvnd6RuMAZxO/VGs7y8PMaMGUNlZSVxcXGUlZU1W75w4ULeeOON7nkDwusSo+1csfh3/PdfVsP48WaX41Mk8D0kwyN7d3jkDz74gN69e58W+N0hIiKC2bNnN97RDMbdxuvWrevQLxPhm/pF23FYQyg8If33p5LA95AMj9x9wyP//ve/Z+nSpQD85Cc/YebMmQB8+umnXHvttYAxJERxcTF33303e/bsISMjg5/97GcAVFRUtNuO7Vm0aFHj0BAAq1atYs6cOY2Dzwn/kxhjjOEkwyuczq8HT+PBHhpN88GyVhfJ8MjdNzzytGnT+MMf/sAdd9xBdnY2tbW11NfXs27dutOGonjsscfYsmULGzduBIwunY60Y3vmzJnDD3/4Q0pKSoiLi2PFihX86Ec/8mgbwrfERYQxbf9GHvzorzhzz8f64nKzS/IZcoTvIRkeubmuDI88YcIEcnJyKC8vJywsjEmTJpGdnc0XX3zRobGH2mvHjggNDWXBggWsXLmS4uJiNm7cyOzZsz3ejvAdVosiMtJO6rECHDt3mV2OT/HJI/yOPgClrSPxniLDI7fNk+GRG97DCy+8wOTJk0lPT+fzzz9nz549HRq3pr127KhFixbx8MMPo7Xm4osvxmazdWo7wne4BjQ86vCQuYX4GJ88wvflB6DI8MjNdXV45GnTpvH4448zbdo0pk6dytNPP01GRsZpv7TaqqGrZsyYwe7du3nyySdZtGhRj+xDeFfIoCTA/ahD94GO8NHA92UyPHL3Do88depUDh8+zKRJk0hMTMRut7fYVnFxcUyZMoUxY8Y0nrTtqFtvvZWkpCSSkpKYNGnSacstFguXXXYZJSUlTJs2zaNtC98UG9+b4/ZILI56aOXzHoxkeGQhOkA+d/7lyc9zmXXleYwo3g/ffQeZmWaX5DUyPLIQIqj0c99tC0g/fhM+edJWCCG6IjHazlujzuXI2PFc1d7FH0FEAl8IEXD6xYSxasxMNvaN4KoWznUFK+nSEUIEnIQmDzP35fOU3iaBL4QIOFFhISS4ajh729dUr2p5CJJgJIEvhAg4SinSHcdZ9tZDWO+9x+xyfIYEficopbjuuusapx0OB/Hx8R6PsNgwMFhX1xFCtKDhUYcF8qjDBhL4nRAREcGWLVuorq4G4JNPPmm8Q1UI4Rsi+ydQa7URUn4CKirMLscnSOB30oUXXsj7778PwGuvvdbslvzS0lIWLlxIeno6Z599duNwySUlJcyePZvMzExuvfXWZieT/v73v3PmmWeSkZHBrbfe2uxOVyGE5xJjenFYrsVvxv8DX6nWv5599uR6zz7b9roeuvrqq1mxYgU1NTVs3ryZs846q3HZAw88QGZmJps3b+aRRx7h+uuvB+DXv/4155xzDhs2bGDBggWNY8dv376d119/nS+//JKNGzditVp55ZVXutYuQgS5xGg7R6LijAkJfECuw++09PR08vLyeO2115g7d26zZevWrWscUXPmzJmUlJRQVlbG2rVrefvttwGYN28effoYD1r+9NNPycnJYeLEiQBUV1eTkJDgxXcjROBJjLZTGOkO/Px8c4vxEf4f+B29xvaWW4yvbrRgwQLuuusu1qxZQ0lJSZOSTq+pYfTH1sbSv+GGG3j00Ue7tT4hglm/mDDWNxzhHz1qbjE+wv+7dEz0gx/8gPvvv5+xY8c2mz9t2rTGLpk1a9bQt29foqOjm83/8MMPOXbsGGAMubxy5UqOuj+UpaWlLQ59LITouMRoO09OvooZD/4T7rrL7HJ8gv8f4ZsoKSmJH//4x6fNf/DBB7nxxhtJT08nPDycF198ETD69hctWsT48eM599xzSU5OBmDUqFE8/PDDzJ49G5fLhc1m48knn2Tw4MFefT9CBJKEKDvlYRFU1SqcLo3V4vm5ukAjwyML0QHyufNPE37zCSWVdXz7y1kkRNnNLscrZHhkIURQSg2p59XX7iVi1gyzS/EJ0qUjhAhYMX17c/aB71EHgfp6CPLnFfvlEb4vd0OJwCOfN//VNzaSosg+KK2hsNDsckznd4Fvt9spKSmR/4TCK7TWlJSUYLcHR/9voGl2Lb7cfOV/XTpJSUnk5+dTVFRkdikiSNjtdpKSkswuQ3RCvxj33baFu+XmK/ww8G02G0OGDDG7DCGEH0iMDiNfjvAbea1LRym1UCn1N6XUP5RSs721XyFE8EqMtlMo4+k06lDgK6WWKaWOKqW2nDJ/jlJqp1IqVyl1d1vb0Fq/o7W+GVgCXNXpioUQooMSo+1sGDCCN7PmwqRJZpdjuo526SwHngBeapihlLICTwLnA/nAeqXUu4AVOHVQmB9orRsGs/iV+3VCCNGjYsNDyU4dx9eD05l/0RyC/dR7hwJfa71WKZVyyuwzgVyt9V4ApdQK4GKt9aPAaY9+UsaoYY8BH2qtv2ttX0qpW4CGUc4qlFI73T/HAGWnrN503qnL+wI99aiolmrprte0tV5ry9prm9bmNZ2W9pL2Cuj26vU7z1/Tznq+2l6tj8mite7QF5ACbGkyfTnwXJPp64An2nj9HUAO8DRwW0f32+T1z7Y179TlQLan++hKLd31mrbWa21Ze23TRhs1bT9pL2kvaa8Ab6+uXKXT0khErV4cr7VeCiztwv7+2c68lpb3lM7sq6OvaWu91pa11zatzfNWm0l7eUbayzPSXh3U4cHT3F0672mtx7inJwEPaq0vcE/fA6CNLh3TKaWydSsDCInTSXt5RtrLM9Jenump9urKZZnrgWFKqSFKqVDgauDd7imrWzzb/iqiCWkvz0h7eUbayzM90l4dOsJXSr0GTMc4kXAEeEBr/bxSai7wJ4wrc5ZprX/bE0UKIYToOp8eD18IIUT38bvB04QQQnSOBL4QQgSJoAx8GdfHM0qpVKXU80qplWbX4quUUhFKqRfdn6vFZtfj6+Qz5Znuyiy/C3wZ18cz3dRee7XWN/Vspb7Hw7a7FFjp/lwt8HqxPsCT9grWz1RTHrZXt2SW3wU+xrg+c5rOaDKuz4XAKGCRUmqUUmqsUuq9U74Smrw0GMb1WU73tVewWU4H2w5IAg66V3N6sUZfspyOt5foXHt1KbP8bjx87cVxfQJBd7RXsPKk7TAGEEwCNuKfB1Jd5mF7bfNyeT7Hk/ZSSm2nGzIrUD6YAzl5dAXGf76Bbaz/I+A84HKl1G09WZiP8qi9lFJxSqmngcyGO6qDWGtt9zZwmVLqKbw7zIeva7G95DPVqtY+X92SWX53hN8Kb4/r4+88ba8SIBh/MbakxbbTWlcCN3q7GD/QWnvJZ6plrbVXt2RWoBzh5wODmkwnAQUm1eIPpL06T9rOM9JenunR9gqUwPf1cX18jbRX50nbeUbayzM92l5+F/jucX2+Bs5QSuUrpW7SWjuA24GPgO3AG1rrrWbW6SukvTpP2s4z0l6eMaO9ZCwdIYQIEn53hC+EEKJzJPCFECJISOALIUSQkMAXQoggIYEvhBBBQgJfCCGChAS+CBhKqUuUUlopNcI9naGM5y43LJ+ulJrs5ZoWNh3tUCn1kFLqPG/WIEQDCXwRSBYB6zDuTgTIAOY2WT4d8GrgAwsxhrkFQGt9v9Z6tZdrEAKQG69EgFBKRQI7gRkYt6KnA7lAL+AQ8BrwE4yx6oswRh/cATwNJLs3c6fW+kul1IPueanu73/SWi91D2X7ntZ6jHufdwGRWusHlVI3A7cAoe79XofxC+c9oMz9dRlwn3sbK5VSs4DHMQYxXA/8l9a6VimVB7wIzAdswBVa6x3d3GQiCMkRvggUC4F/aa13AaXAGOB+4HWtdYbW+ncY4f5H9/QXwJ/d0xMxwvi5JtsbAVyAMT75A0opWzv7f1trPVFrPQ7jlvibtNZfYfzy+Zl7n3saVlZK2TEegHGV1nosRuj/V5PtFWutxwNPAXd1oj2EOI0EvggUi4AV7p9XuKfbcx7whFJqI0YwRyulotzL3tda12qti4GjQGI72xqjlPpCKfU9sBgY3c76ZwD73L+gwDiin9Zk+dvu7zlASgfeixDtCpTx8EUQU0rFATMxQlcDVozx/R9o56UWYJLWuvqU7QHUNpnlxPi/4qD5QZK9yc/LgYVa601KqSUY5wvaLLud5Q37b9i3EF0mR/giEFwOvKS1Hqy1TtFaDwL2YfS/RzVZr/yU6Y8xRiYEjKt62tnPESDB/bSmMJo/DjIKOOzu+lncxj4b7ABSlFJD3dPXAf9uZ/9CdIkEvggEi4BVp8x7C+gHjFJKbVRKXYXx6MFL3NNTgTuALKXUZqXUNtp5ApPWuh54CPgG42Rs0xOp97nnf3LK/BXAz5RSG5RSaU22VYPxhKw33d1ALoxzDEL0GLlKRwghgoQc4QshRJCQwBdCiCAhgS+EEEFCAl8IIYKEBL4QQgQJCXwhhAgSEvhCCBEkJPCFECJI/H87GOwHzy8BFAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.loglog(ETeV,np.exp(-1. * tauCalc), lw = 2, label = 'Own calculation')\n",
    "plt.loglog(ETeV,np.exp(-1. * tauLIV), lw = 2, label = 'Own calculation with LIV')\n",
    "plt.loglog(ETeV,np.exp(-1. * tau.opt_depth(z0,ETeV)), lw = 2, ls = '--', \n",
    "           color = 'red', label = 'Model')\n",
    "\n",
    "plt.legend(loc = 0)\n",
    "plt.gca().set_ylim(1e-2,2)\n",
    "plt.gca().set_xlabel('Energy (TeV)')\n",
    "plt.gca().set_xlabel('Attenuation')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:basic]",
   "language": "python",
   "name": "conda-env-basic-py"
  },
  "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}