Characteristic quantities.ipynb 7.72 KB
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### This Code aims to draw characteristic quantities versus EGMF and redshift\n",
    "\n",
    "Because for redshift upper than 1 there is no TeV photons, we will use the GeV band\n",
    "\n",
    "Characteristic quantities are: mean delay and mean time delay\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[  0.00000000e+00   1.00000000e-08   1.00000000e-09   1.00000000e-10\n",
      "    1.00000000e-11   1.00000000e-12   1.00000000e-13   1.00000000e-14\n",
      "    1.00000000e-15   1.00000000e-16   1.00000000e-17   1.00000000e-18]\n",
      " [  3.08000000e-02   2.31669167e+00   2.46272189e+00   2.31761730e+00\n",
      "    2.35409581e+00   2.19476260e+00   2.19649884e+00   2.30721120e+00\n",
      "    2.42420603e+00   2.75069472e+00   2.78753184e+00   4.18989200e-01]\n",
      " [  1.00000000e-01   5.03021487e+00   5.17248664e+00   5.17318998e+00\n",
      "    5.02288445e+00   5.06413069e+00   5.17880568e+00   5.20734079e+00\n",
      "    5.04978855e+00   4.95592646e+00   2.65673680e+00   3.95433192e-01]\n",
      " [  5.00000000e-01   5.42277169e+00   5.47854199e+00   5.52659330e+00\n",
      "    5.67044186e+00   5.53964475e+00   5.61695768e+00   5.55298921e+00\n",
      "    5.37630821e+00   4.77957703e+00   2.00147238e+00   2.29919715e-01]\n",
      " [  1.00000000e+00   4.29019100e+00   4.29747752e+00   4.43425212e+00\n",
      "    4.34672581e+00   4.36502152e+00   4.39432230e+00   4.30122813e+00\n",
      "    4.07557647e+00   3.40072648e+00   9.92819407e-01   9.84245819e-02]\n",
      " [  2.00000000e+00   4.36694310e+00   4.32006031e+00   4.29913362e+00\n",
      "    4.38270454e+00   4.28456317e+00   4.35261766e+00   4.23559394e+00\n",
      "    3.81561053e+00   2.54938196e+00   4.67573547e-01   5.24824684e-02]]\n",
      "========================================\n",
      "[[  0.00000000e+00   1.00000000e-08   1.00000000e-09   1.00000000e-10\n",
      "    1.00000000e-11   1.00000000e-12   1.00000000e-13   1.00000000e-14\n",
      "    1.00000000e-15   1.00000000e-16   1.00000000e-17   1.00000000e-18]\n",
      " [  3.08000000e-02   1.60470842e+15   1.57648280e+15   1.58144197e+15\n",
      "    1.63046762e+15   1.56922509e+15   1.60240329e+15   1.54999595e+15\n",
      "    1.55111818e+15   1.36769539e+15   3.97493300e+14   5.75147193e+12]\n",
      " [  1.00000000e-01   5.41158326e+15   5.68576605e+15   5.70251767e+15\n",
      "    5.52322164e+15   5.44810792e+15   5.64840914e+15   5.62853346e+15\n",
      "    5.25726175e+15   4.48577135e+15   9.91017717e+14   1.36177755e+13]\n",
      " [  5.00000000e-01   2.41804767e+16   2.37623266e+16   2.44916015e+16\n",
      "    2.49931585e+16   2.48178986e+16   2.50055740e+16   2.43188573e+16\n",
      "    2.33421908e+16   1.80598946e+16   2.16234523e+15   2.41590545e+13]\n",
      " [  1.00000000e+00   3.42065092e+16   3.48746289e+16   3.55725057e+16\n",
      "    3.56122244e+16   3.43933196e+16   3.52897430e+16   3.31823713e+16\n",
      "    3.21234051e+16   2.09827453e+16   1.19849654e+15   9.39172158e+12]\n",
      " [  2.00000000e+00   6.02478400e+16   5.73622982e+16   6.04580297e+16\n",
      "    5.94871078e+16   5.89450726e+16   6.02765357e+16   5.42390170e+16\n",
      "    4.79504034e+16   2.43448368e+16   3.96613703e+14   5.74073019e+12]]\n"
     ]
    }
   ],
   "source": [
    "from matplotlib.pyplot import figure, show\n",
    "from numpy import zeros, size, nditer, average\n",
    "from modules.read import ReadResults\n",
    "from modules.constants import degre\n",
    "\n",
    "Redshifts=[\"0.0308\",\"0.1\",\"0.5\",\"1\",\"2\"]\n",
    "EGMFs=[\"08\",\"09\",\"10\",\"11\",\"12\",\"13\",\"14\",\"15\",\"16\",\"17\",\"18\"]\n",
    "powerlaw_index=2\n",
    "\n",
    "theta_mean=zeros((size(Redshifts)+1,size(EGMFs)+1))\n",
    "dt_mean=theta_mean.copy()\n",
    "\n",
    "it=nditer(theta_mean, flags=['multi_index'], op_flags=['readwrite'])\n",
    "while not it.finished:\n",
    "    i=it.multi_index[0]\n",
    "    j=it.multi_index[1]\n",
    "    if i==0:\n",
    "        if j==0:\n",
    "            theta_mean[i,j]=0\n",
    "            dt_mean[i,j]=0\n",
    "        else:\n",
    "            theta_mean[i,j]=10**(-float(EGMFs[j-1]))\n",
    "            dt_mean[i,j]=10**(-float(EGMFs[j-1]))\n",
    "    else:\n",
    "        if j==0:\n",
    "            theta_mean[i,j]=float(Redshifts[i-1])\n",
    "            dt_mean[i,j]=float(Redshifts[i-1])\n",
    "        else:\n",
    "            fileId=\"z=\"+Redshifts[i-1]+\"-EGMF\"+EGMFs[j-1]+\"-lambda_B=1Mpc-Dominguez\"\n",
    "            weightini,time_delay,theta,Esource = ReadResults(\"Simulations/\"+fileId,cols=[1,3,6,7])\n",
    "            weight_source = (Esource/min(Esource))**(1-powerlaw_index)\n",
    "            weight = weightini* weight_source\n",
    "            cond= (Esource<1e3) & (Esource>1e0) # GeV band\n",
    "            theta_mean[i,j]=average(theta[cond],weights=weight[cond])*degre\n",
    "            dt_mean[i,j]=average(time_delay[cond],weights=weight[cond])\n",
    "    \n",
    "    it.iternext()\n",
    "\n",
    "print theta_mean\n",
    "print \"========================================\"\n",
    "print dt_mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#%matplotlib inline\n",
    "from matplotlib.pyplot import figure, show\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from numpy import log10\n",
    "\n",
    "fig = figure()\n",
    "ax = fig.add_subplot(111,projection='3d')   \n",
    "\n",
    "B=theta_mean[0,1:]\n",
    "i=1\n",
    "for z in Redshifts:\n",
    "    ax.plot(theta_mean[i,1:],log10(dt_mean[i,1:]),log10(B),\"--*\",label=\"z=\"+z)\n",
    "    i+=1\n",
    "    \n",
    "ax.legend(loc=\"best\")\n",
    "ax.set_xlabel(\"$\\\\theta$ [deg]\")\n",
    "ax.set_ylabel(\"Time delay [s]\")\n",
    "ax.set_zlabel(\" log B [Gauss]\")\n",
    "show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#%matplotlib inline\n",
    "from matplotlib.pyplot import figure, show\n",
    "from numpy import log10\n",
    "\n",
    "fig = figure()   \n",
    "ax1 = fig.add_subplot(121)   \n",
    "ax2 = fig.add_subplot(122)\n",
    "\n",
    "B=theta_mean[0,1:]\n",
    "i=1\n",
    "for z in Redshifts:\n",
    "    ax1.plot(B,theta_mean[i,1:],'--*',label=\"z=\"+z)\n",
    "    ax2.plot(B,dt_mean[i,1:],'--*',label=\"z=\"+z)\n",
    "    i+=1\n",
    "\n",
    "ax1.set_xscale('log')   \n",
    "ax1.set_yscale('log')\n",
    "ax1.grid(b=True,which='major')\n",
    "ax1.legend(loc=\"best\")\n",
    "ax1.set_xlabel(\"B [Gauss]\")\n",
    "ax1.set_ylabel(\"average arrival angle [degre]\")\n",
    "ax2.set_xscale('log')   \n",
    "ax2.set_yscale('log')\n",
    "ax2.grid(b=True,which='major')\n",
    "ax2.legend(loc=\"best\")\n",
    "ax2.set_xlabel(\"B [Gauss]\")\n",
    "ax2.set_ylabel(\"average time delay [s]\")\n",
    "show()\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
}