{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Beware this code take time to run!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.01        0.01217357  0.01481958  0.01804072  0.021962    0.0267356\n",
      "  0.03254677  0.03962104  0.04823295  0.05871672  0.07147921  0.08701572\n",
      "  0.1059292   0.12895366  0.15698264  0.19110393  0.23264171  0.28320803\n",
      "  0.34476529  0.41970246  0.51092774  0.62198149  0.75717355  0.92175056\n",
      "  1.12209954  1.36599578  1.66290459  2.02434862  2.46435506  3.        ] [ 19199.20665593  17508.27031736  15966.26022101  14426.43951218\n",
      "  12796.96868216  10940.54707206   8932.04599858   6836.51600451\n",
      "   4771.76094894   2954.40799888   1607.70442167   1091.53593533\n",
      "    835.45280584    675.81811682    562.01738481    471.70846909\n",
      "    399.57803019    338.47728559    284.08836902    234.08272762\n",
      "    182.49932448    161.88596902    139.6835118     117.23818033\n",
      "    102.09606623     94.83681866     87.28526624     82.58799388\n",
      "     92.24970053    102.09606623]\n"
     ]
    }
   ],
   "source": [
    "from modules.analytic import distance, lambda_gg\n",
    "from numpy import logspace, log10, shape, zeros, nditer\n",
    "\n",
    "z_tab=logspace(-2,log10(3),num=30)\n",
    "E_save=zeros(shape(z_tab)[0])\n",
    "it=nditer(E_save, op_flags=['writeonly'])\n",
    "\n",
    "for z in z_tab:\n",
    "    diff_save=1e5\n",
    "    for E in logspace(1,5,num=1000):\n",
    "        diff = abs(lambda_gg(E,z)[0]-distance(z)[1])     \n",
    "        if diff < diff_save:\n",
    "            it[0] = E\n",
    "            diff_save = diff\n",
    "    it.iternext()\n",
    "    \n",
    "print z_tab, E_save    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[  1.91992047e+04  -1.00002070e+00]\n"
     ]
    }
   ],
   "source": [
    "from numpy.random import normal\n",
    "from scipy.optimize import curve_fit\n",
    "\n",
    "def Ecut(z,E0,alpha):\n",
    "    return E0*(z/1e-2)**alpha\n",
    "\n",
    "E_approx=Ecut(z_tab,E_save[0],-1) + 0.2 * normal(size=len(z_tab))\n",
    "\n",
    "param, pcov = curve_fit(Ecut, z_tab, E_approx)\n",
    "\n",
    "print param"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEWCAYAAAB1xKBvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmcFdWZ//HPwy6iEvdE2gFFwG3E6OASlxZHZYxrDEaI\nUYyimFETVNSJUZpR4xK3KKMGHUSTAC4xETHueB2MMS4IGlkElfwaTVCjKIsgy/P749yGFunuuktV\n3er7fb9e9yVVt27VgWPfp895zmLujoiISEvapF0AERHJBgUMERGJRAFDREQiUcAQEZFIFDBERCQS\nBQwREYlEAUNERCJRwBARkUgqKmCYWa2ZTTWz283s4LTLIyIi61RUwADWAIuBjsCClMsiIiKNxB4w\nzGysmS00szfWOz/AzGab2Vwzuzh/eqq7HwlcAoyKu2wiIhJdEi2Mu4EBjU+YWVtgdP78LsAgM9vZ\n1y1stYjQyhARkQrRLu4HuPtUM+u+3ul+wDx3nw9gZhOBY82sD3AE0BW4Ne6yiYhIdLEHjCZsB9Q3\nOl4A7OPu1wC/b+nDZqYldkVEiuDuVuxn00p6l/yF7+6pv0aOHJn6vQr5XJRrm7um0Peaur6c/26V\nUHeVUn/Fvl/I+Uqou3KXoxLqrqVrinlvQ+dLlVbAeA+oaXRcQ4Gjourq6sjlcuUsU8Fqa2tTv1ch\nn4tybXPXFPpeU9fPnz+/xXLErZx1V8r9yll/xb5fyPlKqDvQz16U9xqfz+Vy1NXVtViGllg5ok6L\nDwk5jEfcfff8cTtgDnAo8D7wEjDI3WdFvJ8nUW6Jx5AhQxg3blzaxZAiqO6yzczwSu6SMrMJwAtA\nLzOrN7PT3H0VcA7wBDATuC9qsJDsGzJkSNpFkCKp7qpbIi2McjMzHzlyJLW1tWXvWhARaW1yuRy5\nXI5Ro0aV1MLIbMDIYrklyOVyCvQZpbrLtorvkhIRkdYhswGjEkZJSXH0G2p2qe6yKVOjpMpNXVIi\nIoVTl5RkjlqG2aW6q24KGCIiEklmA4ZyGNmlfvDsUt1lk3IYGSy3iEialMOQzFHLMLtUd9VNAUNE\nRCJRl5SISJWo2i4pJb1FRKJR0juD5ZZA6xFll+ou26q2hSEiIslSC0NEpEqohSEiIonIbMBQ0ju7\nVG/ZpbrLJiW9M1huCZQ4zS7VXbaV2iWlgCEiUiWUwxARkUQoYEji1A+eXaq76qaAISIikSiHISJS\nJao2h6FhtSIi0WhYbQbLLYGGZmaX6i7bqraFISIiyVILQ0SkSqiFISIiiVDAkMRpsEJ2qe6qmwKG\niIhEohyGiEiVUA5DREQSoYAhiVM/eHap7qpbZgPGcYeN4O67c6xZk3ZJREQqW9XP9F7c4WuM7zyU\nX/iF9D5gKw44AA44APbeGzp1SruEIiKVp2pzGF3mzeDMwUuYY725od3FLHn3Q4YPhy22CIHjkktg\n8mT4+OO0Syoi0jpkNmBQUwP/8z+0eX0GvbdbwpUnvs7LL8PChTBqFHTuDL/8JXTvDrvtBsOGwW9+\nA/PnQwYbVa2K+sGzS3VX3dqlXYCS5QNHgy5d4NBDwwtg1SqYMQOefx4mTYIRI6BtW9Z2YR10UAgo\nbbIbOkVEEpHZHEakci9aBCtXwlZbrT3lDu++GwLI1Knw3HPwySdwyCHQv3947bQTWNG9fCIilanU\nHEbrDhgPPwynnQZDh8KFF34pcDRWXw/PPgtTpsAzz4Sg0hA8+veH7bcv819ARCQFVZv0juTYY0N/\n1JIl0KcPXHwxfPjhVy6rqYFTToFx4+D//T/I5UJ31eOPh1FXPXvCT34SAsrKlYn/LVod9YNnl+qu\nurXugAHrchzTp68LHB980OTlZiFAnHkmTJwYkui/+11onFxyCWyzDQweHN5btCjBv4eISMpad5fU\nhvzzn2HsbZHefz8M1500Cf7v/6BfPzjmGDj6aOjRo+jbiojETjmMFC1dCk89FYLH5MnwL/8S0iWD\nBsEmm6RdOhGRL1MOo1yGD4eLLmq2u2p9G28Mxx0HY8fC3/8OV1wR8h7bbx+6tF55JcbyZpj6wbNL\ndVfdKi5gmNnGZvaymX070QdfcAEsWxZyHAUGDghzOwYMgIcegpkzw4TBgQPhm9+EO+6Azz6Lp9gi\nIkmpuC4pMxsFLAZmufujTVwTX5fUggVwzTUwfjycfTZcdVXRt1qzBp5+GsaMCcN1TzghtDz69Stj\neUVEIqr4LikzG2tmC83sjfXODzCz2WY218wuzp87DJgJfHXsa1K6dYPRo+H112HXXUu6VZs2cPjh\n8OCDMGtWGH114olw2GHw4otlKq+ISEKS6JK6GxjQ+ISZtQVG58/vAgwys52Bg4F9gcHAULMU51t3\n6xbGz5bJttuGYblz54auqoEDw8iq114r2yMyQ/3g2aW6q26xBwx3nwp8st7pfsA8d5/v7iuBicCx\n7v4zdx8OjAfGVMRQqA25++6CcxwN2rcP3VJz54aWxpFHhuAxc2aZyygiUmZpLT64HVDf6HgBsE/D\ngbvf09INhgwZQvfu3QHo2rUrffv2pba2Flj3W1Asx6tXk3v4YTjvPGrPPhsuvJBc/tu+0Pudd14t\np58O55+fY//94eijaxk5EhYsiLH8FXDccK5SyqPj6Me1tbUVVR4dN3+cy+UYN24cwNrvy1IkkvQ2\ns+7AI+6+e/74BGCAuw/NH58M7OPu50a8X/qNj/r6kByfMAHOOCOsVbX11kXf7rPP4Oab4ZZb4Pjj\n4fLLwyR1EZFyqfikdxPeAxp/HdYQWhmR1dXVrY2kqWhYcmTGjDCDb8yYkm636aYhSLz1Fmy5ZRiO\ne9ddrXPvjlTrTUqiusumXJa2aN1AC6MdMAc4FHgfeAkY5O6zIt4v/RZGzP7617Ag4rbbhsDxjW+k\nXaLyadwdJdmiusu2im9hmNkE4AWgl5nVm9lp7r4KOAd4gjCM9r6owSJT3IveI3a33eAvfwlzNvr2\nDdNCWkuM1BdOdqnuqlvFTdyLwsx85MiRa5NwFWv2bPjWt9blOJrYj6Mlr74aWhs77wy33170bUSk\nSuVyOXK5HKNGjdLigxWtITk+cWJJgWP5chg5Eu69NwSN446LoawJUbdGdqnusq3iu6Sq3ob243j+\n+YJv06kTXHttmDU+YkRocXyy/uwWEZEYNdnCMLPFET7/D3ffqbxFallmuqQ2pL4+DIPaaKOib7F0\nadg88OGH4fe/D7sCiog0JfYuKTOb7u59m/1whGvikKkuqRj94Q9h1vj48fDv/552aUSk0sXZJfWd\nCJ+Pco1Ecf/9Te453pTjjgtdVIMHh49nhcbyZ5fqrro1FzAuNLMDmvuwu79T5vJUr/32CzmO3r0L\nChwHHRR2/Rs+PCTDRUTi0lzAeAv4hZn9zcyuM7M9kypUFKnP9C63xjPHGweOZcta/Ogee8DUqXDD\nDTBqVOXP18hc3knWUt1lU2IzvfOztE8Cvgd0JqwkO8Hd3yr56UWqihxGfX2Y4n355WE7vwgWLgy7\n/u2/f1iTKuLHRKRKlJrDKGgeRr6VcTewu7un9nVUFQGjSJ9+GnIbW28d5mx07Jh2ib5KY/mzS3WX\nbbHPwzCzdmZ2jJmNBx4HZqNkd7qmTm0yx7HZZvDYY7ByJRx1FCyOMjhaRCSCJgOGmR1uZmMJK8sO\nBSYDO7r7Se7+cFIFbEqry2EU4qmnmk2Od+oEDzwAPXpA//4FDbxKhH5DzS7VXTbFnsMwsynABOB3\n7l7cCnoxUZcUX96PY+jQDS454g6XXgqTJ8MLL0CXLimVVUQqQmxdUu7e393vdPePzexAMzst/8Ct\nzKxHsQ+UMll/VNVpp33lEjO46ir4t38LS4msWZNCOTegaluGrYDqrrpFyWHUARcD/5U/1QH4TYxl\nkkI0BI5Jkzb4thncdlsYQTVqVMJlE5FWJcqw2hnAnsCr7r5n/tzr7v6vCZSvqTKpSyqqFSugY0cW\nLgx7a1x/PQwcmHahRCQNSaxWu8Ld13ZmmNnGxT5MErZ8OfTqBRddxDb2Ab//PfzoR/Daa2kXTESy\nKErAeMDMfgV0NbMzgWeAu+ItVsuqepRUVJ06haXUly6FPn345n0Xc9fVH3LccaGLKi2qt+xS3WVT\nont6m9nhwOH5wyfc/amSn1wCdUkVodFGTg8e+Etu/uhkpkyBDh2SL4omf2WX6i7bEp3pXSkUMEpQ\nX8+aFSs5YcQObLEF3HlnSIyLSOsXWw7DzM4ws4saHb9nZovNbImZnV3sAyVlNTW06bkDv/41vPQS\njB6ddoFEJCuay2EMA8Y2Ov7A3TcBtgIGxVoqiV2XLmHHvquugqefBt55p+D9OIqlfvDsUt1Vt+YC\nhrn7R42OHwBw98+B4vcXlYrRowdMnAjf/z68u7Dzuj3HEwocIpItzS0NMs/de27gfBtgnrvvEHfh\nmqIcRnndcUdYDv3FF2HTT9clxznjDLjoIthii7SLKCJlEOc8jKfM7Mr1HmbAFcCTxT6wXDSstnyG\nDYN99oFLLuHLS44sXarlbkVagSQWH+xCmG/xb8CM/Ok9gFeAM9w9tW8StTDKb9Ei2G23sJbhgQfG\n+ywNzcwu1V22ldrCaNfUG+6+BDjJzHYEdgUcmOXu84p9mFSurl3DiKnTTw+Ni42ay1LNnw+dO4dd\nmkSkajTXwvi6u/+92Q9HuCYOamHEZ+BA6NkTrr66mYvGjAn9V2ecEZZVV+AQyYQ4cxiPRvh8lGsk\nQ269FcaObWG9qTPPXJfj6NMnJMY/+CCxMopIOpoLGHvkJ+o1+QK2Saqgkoxtt4XrrgtdU6tWNXPh\n+snxfv3giy8iPUODFbJLdVfdmttAqa27b9LCa7skCyvJOOWUsHnfDTdEuLghcMycmc7CVCKSGK0l\nJRs0fz7svXfY2rVXrxJu5K7FqkQqRBL7YUgV6t4dLrssbBde0tauRx2lmeMirYQChjTpnHNCWmLM\nmBJucscdX1lyRP3g2aW6q25R9vS+0cx2TaIwhdBM7/i1bQt33RVaGgsWFHmThhzH9OnrAseECWUt\np4g0L7ENlMxsKDAEaE9YvXaCu39a8pNLoBxGsv77v+Hll2HSpDKkI+rrYd48OOSQspRNRKJLbAMl\nM+tDCByDgeeBO9392WIfXAoFjGR98QXstRdceimcdFLapRGRYiWS9DaztkAfYGfgQ8LaUueb2X3F\nPliyo0MH+N//hZ/8BD76qOXrW7LBrsRVq8L0ciXHK5q6gatblBzGTcAc4EjgKnffy92vdfejgb5x\nF1AqQ79+Yd+MYcNg9eoYHvD55yFRov04RCpWlBzGacD97r50A+91dfdFcRWumTKpSyoFy5aFUbLd\nusHdd4ekeNnVr7cfx4UXhlmEIlKy2HMYZrYXYaXaxj4F/ubuzS0eERsFjPQsWwbHHhu+w++9F9o1\nud5xiRoCx8EHw4knxvQQkeqSRMB4EdgLeD1/anfgTWAz4Gx3f6LYhxdLASNdn38Oxx8Pm24Kv/0t\ntG9f2Oe1p0J2qe6yLYmk9/tA33zuYi9C3uId4DDgumIfLNm10Ubwhz+ENQdPOinymoPls2yZchwi\nKYgSMHq7+5sNB+4+E+jj7m/z1a4qqRKdOsFDD8HKlaHHqJCgUfJvqH/+s5LjKVHrorpFCRhvmtnt\nZnawmdWa2W3ATDPrCKyMuXxSwTp2hAcfhDZt4IQTYMWKhB586KFfnjmuwCGSiCgB41TgbeAnwI8J\n3VGnEoJF//iKJlnQoQPcd19ocRx/PCxf3vJnyjKWf0NLjrz9dun3lWZpHkZ1azZgmFk74I/ufr27\nH59/Xe/uy9x9jbsvLmdhzKxPvjXzgJkNK+e9JT7t24floTbbLIyg+vzzBB/eEDjefBN22CHBB4tU\nnyijpJ4BTkhyvoWZtQHucfcfNPG+RklVoFWrYMgQ+Mc/QlK8S5e0SyQijSUxSmop8IaZjTWzW/Ov\nWwoo4FgzW2hmb6x3foCZzTazuWZ2caPzRwOTgT9GfYZUhnbt4J57oGdP6NsXnn8+7RIBV16pHIdI\nmUQJGA8BlwHPAa8Ar+ZfUd0NDGh8Ir821ej8+V2AQWa2M4C7P+LuRwLfL+AZUiHatg1bYFx/PQwc\nCCNGfDWvkWg/+KmnKjleRsphVLcW5+m6+zgz6wxs7+6zC32Au081s+7rne4HzHP3+QBmNhE41sy2\nBr4DdAQebe6+Q4YMoXv3cNuuXbvSt2/ftUP+Gv6n1nF6x127wuuv13L22dC7d46f/hTOOiu8P336\n9OTKU1NDbuBAOPhganM56NOH3BFHwBlnUNu/f8X8e+lYx3Ec53I5xo0bB7D2+7IUUXIYxwC/ADq6\ne3cz2xMY5e7HRH5ICBiPuPvu+ePvAke4+9D88cnAPu5+bsT7KYeREe5hWagf/xjOPjsskd6hQ4oF\nqq+HP/4RzjorxUKIpCOJHEYdsA/wCYC7vwaUOhxF3/ZVwgwGDQqjX195BfbdF954o+XPxaamRsFC\npEhRAsbKDYyQWlPic98Dahod1wAFbQKqLVqz5RvfgMmT4T//Ew44IMe118a0THopHnxQOY4W6Gcu\nm3Jl2qI16kzv7wPtzGwnM7sVeKHE574C7GRm3c2sA/A9YFIhN6irq1vbZyfZYAannw6/+hU88QTs\ntx9MnZp2qfLcw7AuJcelFaqtrU0sYJwL7AqsACYAnxFmfUdiZhMIAaaXmdWb2Wn5ZdHPAZ4AZgL3\nufusQgsv2XTSSbU8/TSccw6cfDIcc0yYd5cqM7j55tB3tnixAkcT9EtadYu8p3clMTMfOXIktbW1\n+h8445Yvh9tuC1tfHHUUjBoV0gypa9iPo0MHuOmmtEsjUpJcLkcul2PUqFGx74fRG7gQ6M66Ybju\n7qmtI6VRUtmW28CeCosWwXXXhTkcZ5wBl1wCm2+eTvm+xD20PgTQfhhZl8QoqQeAacDPgBGNXiJl\n07Ur/PznYQTVokXQu3cIIImuS7UhTQWLzz5LthwiFSDqKKnb3f0v7v5K/lXITO9YaJRUdjX3G+p2\n28GYMSEZ/uKL0KsX3H57BQSOxt57D3r0qMoch1oX2VSuUVJRuqTqgA8JS4Ss3fHA3T8u+elFUpdU\n9XjxxdDyeOmldZP/unZNu1Ssy3FMnBj60C68MGx0LlLBkuiSGkLIYbzAunWkUm9hSHYV0jLcd1+Y\nNAmefhpmzYIddwy/2P/97/GVL5L19+Po3TuMFW7l1Kqvbi0GDHfv7u491n8lUTiRBrvtBvfeC6++\nGrqndt01TNieNy/lgjUEjhkzoF+/lAsjEq8WA4aZbWxml5nZnfnjnczsqPiL1jzlMLKrlH7w7t3h\nlltgzhzYZpsw+e9734Np08pWvOLU1MDXvpZyIeKnHEY2JZnDuJ/QBXWKu+9qZhsDL7j7HiU/vUjK\nYUiDxYvhzjvhxhtDgvz88+HII8M+4xXh6afhySdDjmPrrdMujVS5JHIYO7r7tcAXAO6+tNiHiUB5\n+8E32SQEiXffDbnnkSNh553DfI5ly8r2mOL16QNLl4ZCtYJRVWrVV7coAWOFmW3UcGBmO9JotJRI\nJWjfHgYPDivi3nknPP546L762c9STpB36/bl5HjDkiOaxyEZFHV588eBbmY2HpgCXNzsJxKgHEZ2\nxdkPbgYHHRT2FP/Tn+CTT2CXXcJe4zNmxPbYljUeVeUeIlwGKYeRTYnlMADMbEtg3/zhi+7+UclP\nLoFyGFKIjz8OK+SOHh3yHOeeGxY8bNfifpMirUsSOQzc/SN3n5x/pRosJPuSbhluvjn813+FPMdZ\nZ4UEeY8eYULgBx8kWpTmvfpqxec41KqvbpUylkQkdh06wEknhW0vJk2Cd94J8+1OOSXMJE/dk09q\nWXWpaE12SZlZD3d/N+HyRKIuKSmXjz+GsWPDEutbbhn26DjxROjUKaUCackRiVGpXVLNBYxX3X0v\nM5uS5lLmG6KAIeW2ejU89ljIc0ybFr6rhw2D7bdPqUANgWPGjNAkEimDOHMYbc3sUsJOeeeb2QWN\nXucX+8By0Sip7KrEemvbNmzg9Pjj4ft52TLYc0/4zndgypQwsClRDaOqKuzfqhLrTloW+yip/MZJ\nxwM/Bu5Y/313H1Xy04ukFka2ZWUTniVL4De/Ca2ONWtCd9UPfhAmC6Zu1apUhnllpe5kw+Lskvqx\nu//SzC539/8uuoQxUMCQJLnDc8+FwDFlStiH/Ec/CvnpVKxeDbvvDkcfrRyHFCTOLqkf5v97fLE3\nF2kNzKC2Fh58MKQUNt00HB9xRFjRPPHfXdq2DQ9uWFZdo6okIc0FjJlmNhfobWZvrPd6PakCSuuT\n5X7wmhq48kr4299g0CAYMSL8sj92LCxfnnBBGpZVbwgcd94Z+2OzXHdSumZnepvZtsCTwNHAl5ox\n7j4/1pI1Q11S2daa+sHd4ZlnwmTAadPCjoBnn53CwrT19WF9ql13jfUxranuqlFsOYz1HtIB6JU/\nnOPuK4t9YDkoYEglmjkTbr4ZHngAvvtdGD48rGMlUiliXxrEzGqBucBt+ddcMzu42AeWi4bVSqXZ\nZRcYMwbeeiv0GPXvD//xHylPo/joo7AuinIcVS3JDZSmAYPcfU7+uBcw0d2/WfLTi6QWRrZVS7fG\n8uXw61/D1VeHpdYvvzwkyxP14YdQV1e2mePVUnetVRKLD7ZrCBYA7v4WoHU+RVrQqRMMHRq2kz3l\nlPDngw8OOY/Eft/ZaqsN78fxkdYQlcJFaWHcDawGfkNIfH8faOPuP2z2gzFSC0OyaNWq8Iv+lVfC\nFluEFsfhh4dhu4lpWHJk6FDo2zfBB0sliD3pbWYdgXOAb+VPTQVuc/fUdt1TwJAsW706JMavuAI2\n3jgEjm9/O+HAIVUpiS6pYe5+g7t/J/+6CRhW7ANFqn2wQtu2YZn1N96Aiy6CSy+FAw6AN99MuWD/\n+EeLyfFqr7tqFyVgDNnAudPKXA6RqtOmTRh++9prYbmR2lq47LKEJwA29swz2o9DmtXcWlKDgMHA\ngYRuqAabAKvd/dD4i7dh6pKS1uj998P2sX/9a9hSNpXBSNqPo1WLc/HBfwF6ANcAF7NupvdiYIa7\nryr2oaVSwJDW7OGHw8q4RxwB110XtphNXEPgeOghmDsXunRJoRBSbrHlMNz9b+6ec/d93f25/J9z\n7v5qmsGigSbuZZfqrXnHHhvyGRttFFb6mDgxxf045sz5UrBQ3WVTkhP3Fjc67AC0B5a4+6YlP71I\namFkmyZ/Rffii2EEbLducPvtYQJgmtbWnbuGdWVQImtJNXpYG+AYYF93v6TYh5ZKAUOqycqVcMMN\noXvqkEPgmGPCMNwtt0yxUKecAl//unIcGZPEsNq13H2Nu/8BGFDsA0WkMO3bwyWXhDWqjjkGJk2C\nHXeEAw+EX/wi9Bol7qqrYPFi7cdRZaIsPnhCo9dAM7sG+DyBskkrpX7w4my5JZx6Kvzud7BwIfz0\np/DOO3DooeF7e8QImDo1TAyMy9q6q6mB224L+3E0BI4rrojvwVIRorQwjgaOyr8OJ4ySOjbOQolI\n8zp1Civh3n57GNA0YUKYNX7eeSHfMWJEQhMBGwcOLTXS6hWUw6gUymGING3OHLjnHrj33pBmGDIk\n7A6YyvBcqShJ7Idxj5l1bXT8NTMbW+wDRSRevXvDz38etpG96ir4059ghx1g4EB49NGwCGJi3OGW\nW5TjaCWidEnt4e6LGg7c/RMgtb0wJPuUw0hG27ZhNdzx42H+fDjssBBAampCl9U77xR+z4LrbsUK\nmD1byfFWIkrAMDPbvNHB5kDb+IokIuXWtSuceSa88AI8+2yYQtGvH5xwQmiBxNbD26nTl5PjffqE\nIV8KHJkUZeLeKcClwP2E5UEGAle5+73xF6/JMimHIVKiJUtCruOmm8L+HOefHwJIuzi3R6uvD1sQ\n/uu/wjAtep20RCbumdmuQH/AgSnuPrPYB5aDAoZI+axeDZMnw403hq6r884L6w5utlnaJZNyS2Ti\nnru/6e63uvvouIOFmR1rZmPMbKKZHRbnsyQdymFUlrZtw/pVzz0X5nhMmwY9esDw4fD221++Nta6\nW7VKXVUVrqCZ3klw94fd/UzCJk3fS7s8ItVk773ht78NKYcOHWC//WD//WH0aPjgg5gfPm2a9uOo\ncIkEDDMba2YLzeyN9c4PMLPZZjbXzC5e72M/A0YnUT5JlhYerHw1NXDttfDee2FTp7/8BXr1gmuu\nqeXee+Gzz2J4aL9+MH16SK4ocFSkRCbumdmBwBLgXnffPX+uLTAH+HfgPeBlYBAwm7AHx5Pu/kwT\n91MOQyRhy5bBI4+EYbq5XNivY/DgMOO8Y8cyP6zxRk5TpsAee5T5AdUp0dVqS2Fm3YFHGgWM/YCR\n7j4gf9yw+u1S4FRCAJnu7r/awL0UMDJMy5tnV0PdffxxyHeMHx+6r/bZJ6wM0vDq2TPkRkq2YEGY\nrl6Wm60zf36IQ1OmhD8//3xZb1+xSg0YcQ6ga8l2QH2j4wXAPu5+LnBrSx8eMmQI3fObA3Tt2pW+\nffuu/RJqSMzpuDKPp0+fXlHl0XFxx0OH1jJ0KDz0UI7Zs2H16lruvx+GD8/xySfQt28tffvCRhvl\n6NkThgyppXPnAp/XrVtZyvvRR7BiRS1TpsCjj+b44gsYMKCW/v1D+XK59P894zjO5XKMGzcOYO33\nZSnSbGGcAAxw96H545NZFzBaupdaGCIV7NNP4fXXQ0qi4TVrFuy0U0hVNLx23bXIeR+33x6aBk3s\nx/Hpp6H18PTT4b8LF4Y90vv3D6+dd67O/Z+y3MJ4D6hpdFxDaGWISMZttlnYr+PAA9edW7EiBJGX\nXgpdQDfeGFIVe+755SDSvXuEL/Ojjgo5jt69YehQ1px/ITPe34rHH4fHHw8DrvbfPyyNcvrpIQVS\n5l6tqpRmC6MdIel9KPA+8BIwyN1nRbiXjxw5ktra2rXNMMmOnHIYmVXuuvv0U3jllRBEGl6LFoUN\nonr2/OqrWzdokx/b+dFH8KeJ9Ww8+hr2emsCD35tKLMGX8nh327PQQdB585lK2bm5XI5crkco0aN\nqvykt5lBPVGXAAAHOUlEQVRNAA4GtgA+AC5397vN7D+AmwlrU/2vu18d8X7qksowBYzsSqLuFi8O\nEwbnzfvq65//DJMKO3UK1xxyCAwYAEfuXs/2f74PLrigOvuaIsrMKKlyUsAQqU7LloVVdj/7LEwy\n7NAh7RJlS5ZzGCWpq6tTl5RIlencGXbbrcAPPfZYiC4bSI5Xi4YuqVKphSGJU5dUdmWy7i66CO66\nC4YObXJUVbVIZPFBEZHMuu66MLtwyRJt5FSizLYwNEpKRArWsORIw2YgVSJTo6TKTV1SIlIS96oc\nTaUuKcmcciTfJB2tpu6qMFiUgwKGiIhEktmAUVdX13p+26kyyjtll+oum3K5HHV1dSXfRzkMEZEq\noRyGZI5ahtmluqtuChgiIhJJZgOGchjZpX7w7FLdZZNyGBkst4hImpTDkMxRyzC7VHfVTQFDREQi\nUZeUiEiVqNouKSW9RUSiUdI7g+WWIJN7Kgigusu6qm1hiIhIstTCEBGpEmphiIhIIhQwJHEarJBd\nqrvqpoAhIiKRZDZgaFhtdmmUTXap7rJJw2ozWG4RkTQp6S2Zo5ZhdqnuqpsChoiIRKIuKRGRKqEu\nKRERSYQChiRO/eDZpbqrbgoYIiISiXIYIiJVQjkMERFJRGYDhmZ6Z5fqLbtUd9mkmd4ZLLcE2oQn\nu1R32VZql5QChohIlVAOQ0REEqGAIYlTP3h2qe6qmwKGiIhEohyGiEiVUA5DREQSoYAhiVM/eHap\n7qqbAoaIiESiHIaISJVQDkNERBJRUQHDzHqY2V1m9kDaZZH4qB88u1R31a2iAoa7v+vuZ6RdDonX\n9OnT0y6CFEl1V91iDxhmNtbMFprZG+udH2Bms81srpldHHc5pHIsWrQo7SJIkVR31S2JFsbdwIDG\nJ8ysLTA6f34XYJCZ7ZxAWcqqnM3zYu9VyOeiXNvcNYW+V8ndF+UuWyXUX7HvF3q+Euhnr+X34qi/\n2AOGu08FPlnvdD9gnrvPd/eVwETgWDPb3MzuAPpmodWh/2mbf6+p6+fPn99iOeKmgFHc+UqoO9DP\nXpT34ggYiQyrNbPuwCPuvnv++LvAEe4+NH98MrCPu58b8X4aUysiUoRShtW2K2dBClDSF34pf2ER\nESlOWqOk3gNqGh3XAAtSKouIiESQVsB4BdjJzLqbWQfge8CklMoiIiIRJDGsdgLwAtDLzOrN7DR3\nXwWcAzwBzATuc/dZcZdFRESKl8m1pEREJHkVNdNbREQqV6sJGGZ2rJmNMbOJZnZY2uWRwmgdsWwy\ns43N7J78z97gtMsj0RXzM9fquqTMrCtwvdakyiYze8DdB6ZdDonGzH4AfOzuj5rZRHc/Ke0ySWEK\n+ZmruBZGGdae+hlh2RFJgdYOy74C63A7oD7/59WJFlS+Iu6fv4oLGBSw9pSZ/cDMbjKzb1hwLfCY\nu2tJzfQUVX8plFOaVsj6bwtYN6eqEr9Pqk2sa/dVXAUXsvaUu//a3Ye7+/vAucChwHfN7KxkSy0N\niq2/rK0j1poVUofAQ8AJZnYbmkuVurjX7ktraZBCNW72QvitZp/GF7j7LcAtSRZKIotSfx8Dw5Is\nlBRkg3Xo7suAH6ZTJImoqbor+Geu4loYTWhdmfnqo/rLPtVhdpWt7rISMLT2VLap/rJPdZhdZau7\nrAQMrT2Vbaq/7FMdZlfZ6q7iAobWnso21V/2qQ6zK+66a3UT90REJB4V18IQEZHKpIAhIiKRKGCI\niEgkChgiIhKJAoaIiESigCEiIpEoYIiISCQKGCJFMrM6M7ugmPfN7E+N/vwLM/urmV1nZqea2dfj\nKK9IqbKyWq1IYszMALzlWa1Fv+/u32p0OBT4mru7mT0L/BX4e5SyiiRJLQwRIL/Ozhwzuwd4A7jM\nzF4ysxlmVtfoukvz100Fejc6f56ZvZm/fnyjW+9iZs+a2dtmdm6j65fk/zsJ6AJMM7MTgb2B35rZ\nNDPrFOtfWqRAamGIrNMT+AGwGfBdd+9nZm2Ah83sQGAZYeG2PYD2wDTCwm4AFwPd3X2lmW2aP2dA\nH6AW2BSYY2a3uftq8q0Pdz/GzBa7+54AZnY2cIG7T4v/rytSGLUwRNb5m7u/BBwBHG5mrwGvEloS\nOwEHAA+5+3J3X8yXV/x8HRhvZt9n3d7WDkx295Xu/k/gA2CbCOWw8vx1RMpLAUNknaWN/ny1u++Z\nf/Vy97H5842/zK3R8beB/wG+Cbyc30cZ4ItG168mWqteK4JKRVLAEPmqJ4AfmtnGAGa2nZltBfwf\ncJyZdTKzTYCjAM8nybd39xxwCaFLqwvFtRQWE7qvRCqOchgi6zTkFZ4ys52BP+cHTC0GTnb318zs\nPmAGoXvppfzn2gK/NrPNCEHil+7+qZk5TbcWvIk/jwPuMLNlwP7uvrw8fzWR0mk/DBERiURdUiIi\nEokChoiIRKKAISIikShgiIhIJAoYIiISiQKGiIhEooAhIiKR/H8KlBdMSMD06wAAAABJRU5ErkJg\ngg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f6ae6708150>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from matplotlib.pyplot import figure, show\n",
    "%matplotlib inline\n",
    "fig = figure()   \n",
    "ax1 = fig.add_subplot(111)\n",
    "ax1.plot(z_tab,E_save,drawstyle='-')\n",
    "ax1.plot(z_tab,Ecut(z_tab,param[0],param[1]),\"r--\")\n",
    "ax1.set_xscale('log')   \n",
    "ax1.set_yscale('log')\n",
    "ax1.grid(b=True,which='major')\n",
    "ax1.set_xlabel(\"redshift\")\n",
    "ax1.set_ylabel(\"cutoff energy [GeV]\")\n",
    "show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$E_{cut}(z) \\simeq 2.10^4 \\left(\\frac{z}{10^{-2}}\\right)^{-1} \\textrm{ GeV}$$ \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}