Monte_Carlo_Gamma_Ray_Simulation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Monte Carlo simulation of low energy gamma rays passing through a medium\n",
    " -*- coding: utf-8 -*-\n",
    "\"\"\"\n",
    "Monte Carlo simulation of low energy gamma rays passing through a medium\n",
    "Created on Mon Jun 17 22:11:48 2019"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loading approximation data for various cross sections...\n",
      "Loaded approximation data.\n"
     ]
    }
   ],
   "source": [
    "import monfgp\n",
    "import monfgc\n",
    "import monfsgc\n",
    "import monfsgp\n",
    "import math\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import time as time\n",
    "import ipywidgets as ipw\n",
    "\n",
    "print(\"Loading approximation data for various cross sections...\")\n",
    "MesCr = np.loadtxt('sef.txt')\n",
    "falphaZ=MesCr[:,0]\n",
    "fb0=MesCr[:,1]\n",
    "fb1=MesCr[:,2]\n",
    "fb2=MesCr[:,3]\n",
    "ffaZ=MesCr[:,4]\n",
    "print(\"Loaded approximation data.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Medium properties:\n"
     ]
    }
   ],
   "source": [
    "print(\"Medium properties:\")\n",
    "mat = input(\"Enter material of the medium ('NaI' or 'Other'): \")\n",
    "rho = float(input(\"Enter density of the medium: \"))\n",
    "neq = int(input(\"Enter no of elements per molecule or compound: \"))\n",
    "\n",
    "nam = [0]*neq\n",
    "ZA = [0]*neq\n",
    "A=[0]*neq\n",
    "alphaZ=[0]*neq\n",
    "b0 = [0]*neq\n",
    "b1 = [0]*neq\n",
    "b2 = [0]*neq\n",
    "faZ = [0]*neq\n",
    "\n",
    "for n in range(neq):\n",
    "    nam[n] = int(input(\"\\nEnter number of atoms per molecule of element {0}:\".format(n+1)))\n",
    "    ZA[n] = int(input(\"Enter ZA of element {0}:\".format(n+1)))\n",
    "    A[n] = int(input(\"Enter A of element {0}:\".format(n+1)))\n",
    "    if mat == 'NaI':\n",
    "        flag = 0\n",
    "    else:\n",
    "        print(\"Element %d \\n %d atoms; Z = %d A = %d \\n\", n+1, nam(n), ZA(n), A(n)) \n",
    "        \n",
    "    alphaZ[n]=ZA[n]/137.036\n",
    "    for j in range(14):\n",
    "        if (alphaZ[n] >= falphaZ[j] and alphaZ[n] < falphaZ[j+1]):\n",
    "            b0[n]=((fb0[j+1]-fb0[j])/(falphaZ[j+1]-falphaZ[j]))*(alphaZ[n]-falphaZ[j])+fb0[j]\n",
    "            b1[n]=((fb1[j+1]-fb1[j])/(falphaZ[j+1]-falphaZ[j]))*(alphaZ[n]-falphaZ[j])+fb1[j]\n",
    "            b2[n]=((fb2[j+1]-fb2[j])/(falphaZ[j+1]-falphaZ[j]))*(alphaZ[n]-falphaZ[j])+fb2[j]\n",
    "            faZ[n]=((ffaZ[j+1]-ffaZ[j])/(falphaZ[j+1]-falphaZ[j]))*(alphaZ[n]-falphaZ[j])+ffaZ[j]\n",
    "            break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NaI (Accurate cross section are used for NaI)\n",
      "\n",
      "\n",
      "Medium properties:\n",
      "\n",
      "Medium Geometry:\n",
      "\n",
      "Enter geometry of the medium (p for plane,c for cylinder):p\n",
      "Enter thickness of the plane(cm):40\n",
      "parallel plane layer of  40.0  cm width.\n",
      " \n"
     ]
    }
   ],
   "source": [
    "if mat == 'NaI':\n",
    "    print(\"NaI (Accurate cross section are used for NaI)\")\n",
    "else:    \n",
    "    print(\"Material = %s \\n density = %3.1f g/cm**3 \\n molecule consisting of %d elements\\n\",mat,rho,neq)\n",
    "\n",
    "print(\"\\n\")\n",
    "print(\"Medium properties:\\n\")\n",
    "print(\"Medium Geometry:\")\n",
    "\n",
    "geo_flag = 1\n",
    "while geo_flag:\n",
    "    geo = input(\"\\nEnter geometry of the medium (p for plane,c for cylinder):\")\n",
    "    if geo=='p' or geo=='c':\n",
    "                geo_flag = 0\n",
    "    else:\n",
    "                print(\"Please enter valid input 'p' or 'c' :\")\n",
    "\n",
    "geom = [0]*3\n",
    "geom[0] = geo\n",
    "\n",
    "if geo == 'p':\n",
    "    geom[1]= float(input(\"Enter thickness of the plane(cm):\")) # thickness of the plane\n",
    "    geom[2]= 0.0                                               # this value is not used\n",
    "elif geo == 'c':\n",
    "    geom[1]= float(input(\"Enter height of the cylinder(cm):\")) # height of the cylinder\n",
    "    geom[2]= float(input(\"Enter radius of the cylinder(cm):\")) # radius of the cylinder\n",
    "\n",
    "if geom[0]=='p':\n",
    "    H=geom[1]\n",
    "    print(\"parallel plane layer of \",geom[1],\" cm width.\")\n",
    " \n",
    "if geom[0]=='c':\n",
    "    H=geom[1]\n",
    "    R=geom[2]\n",
    "    print(\"cylinder of \",geom[1],\" cm height and \",geom[2],\" cm radius\")\n",
    "\n",
    "print(\" \")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Beam Geometry:\n",
      "\n",
      "Enter beam geometry (p for parallel,i for isotropic):p\n",
      "Enter radius of the beam(cm):0.1\n",
      "Enter diaphram radius(cm):80\n",
      "Parallel beam of  0.1  cm radius \n",
      "Input diafragm of  80.0  cm \n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(\"Beam Geometry:\")\n",
    "\n",
    "geoh_flag = 1\n",
    "geoh=[0]*4 \n",
    "\n",
    "while geoh_flag:\n",
    "    geoh[0] = input(\"\\nEnter beam geometry (p for parallel,i for isotropic):\")\n",
    "    if geoh[0]=='p' or geoh[0]=='i':\n",
    "        geoh_flag = 0\n",
    "    else:\n",
    "        print(\"Please enter valid input 'p' or 'i' :\")\n",
    "\n",
    "     \n",
    "if geoh[0] == 'p':\n",
    "    geoh[1]= input(\"Enter radius of the beam(cm):\")\n",
    "    geoh[2]= 0.0                                           \n",
    "    geoh[3]= input(\"Enter diaphram radius(cm):\") \n",
    "elif geoh[0] == 'i':\n",
    "    geoh[1]= input(\"Enter radius of the source(cm):\")\n",
    "    geoh[2]= input(\"Enter the source - medicum distance(cm):\")\n",
    "    geoh[3]= input(\"Enter diaphram radius(cm):\")\n",
    "\n",
    "haz=geoh[0]\n",
    "rf=float(geoh[1])\n",
    "d=float(geoh[2])\n",
    "Rdia=float(geoh[3])\n",
    "    \n",
    "Rmin=Rdia\n",
    "\n",
    "if geom[0]=='c':\n",
    "   if R <= Rdia:\n",
    "        Rmin=R\n",
    "   \n",
    "if haz=='p':\n",
    "    print(\"Parallel beam of \",rf,\" cm radius \\nInput diaphragm of \",Rdia,\" cm \\n\")\n",
    "else:\n",
    "    print(\"Isotropic source of \",rf,\" cm radius at a distance of \",d,\" cm \\ninput diaphragm of \",Rdia,\" cm\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Radiation properties:\n",
      "Enter photon energy (MeV): 1\n",
      "Enter number of photons: 20\n",
      "\n",
      "\n",
      "20  gamma rays of  1.0  MeV energy \n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(\"Radiation properties:\")\n",
    "\n",
    "E0 = float(input(\"Enter photon energy (MeV): \"))\n",
    "nfot = int(input(\"Enter number of photons: \"))\n",
    "print(\"\\n\")\n",
    "print(nfot,\" gamma rays of \",E0,\" MeV energy \\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sg for without plot, t for tracks plotted, and p for interaction points plotted\n",
      "Enter your choice (sg/t/p): t\n"
     ]
    }
   ],
   "source": [
    "# Read graphic options\n",
    "op1='sg'\n",
    "op2='xx'\n",
    "\n",
    "print(\"sg for without plot, t for tracks plotted, and p for interaction points plotted\")\n",
    "opt_flag = 1\n",
    "while opt_flag:\n",
    "    opt = input(\"Enter your choice (sg/t/p): \")\n",
    "    if opt == 'sg' or opt == 't' or opt == 'p':\n",
    "        opt_flag = 0\n",
    "    else:\n",
    "        print(\"Please enter correct choice:\")\n",
    "\n",
    "if opt!='sg':\n",
    "    op1='g'\n",
    "    op2=opt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Plotting tracks of photons\n",
      "monfgp running\n",
      "Simulation is over\n",
      "\n",
      "Average values:\n",
      "\n",
      "\n",
      "Transmission coefficient:  0.2\n",
      "Backscattering coefficient:  0.05\n",
      "absorption coefficient:  0.75\n",
      "\n",
      "Fraction of photons leaving the medium without interacting:  0.1\n",
      "Fraction of photons which suffer more than one interaction:  0.0\n",
      " \n",
      "FILES AND PLOTS\n",
      "  \n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqsAAAFUCAYAAAADN3WgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xd0FGXfxvFrdjebSic06UU6CaGJICBVeBGQIiKCYFcEfVQQFDuCKAqPYMGCgNixoD5IUVqo0pEmvSMlEAipW+b9Y81CTIKUhJ3A93NODuzM7Mxvc2By5Z67GKZpCgAAALAiW6ALAAAAALJDWAUAAIBlEVYBAABgWYRVAAAAWBZhFQAAAJZFWAUAAIBlOf5lP/NaAQAA4EowstpIyyoAAAAsi7AKAAAAyyKsAgAAwLIIqwAAALAswioAAAAsi7AKAAAAyyKsAgAAwLIIqwAAALAswioAAAAsi7AKAAAAyyKsAgiY8uXLKzQ0VBERESpevLj69++vM2fOSJJatGihjz76KNv3JiYmKiIiQh06dMjyvMWLF1diYqJ/20cffaQWLVr4XxuGofDwcEVERPi/Xn/99Wyvt2rVKnXs2FGFChVSwYIFVaNGDT377LM6efKkJGny5Mmy2+0ZzhcREaFDhw75zzF58mTVrl1bYWFhKlGihB5++GHFx8f797/44osyDENvv/12hmuPGzdOhmHoxRdflCQtWLBANpst07WWLVummjVr+l/b7XaFhIT4X48cOTJDnfnz51dUVJR+/vlnSVJycrKqVKmiqVOnZrj+Sy+9pCZNmsjr9Wb7/QGA3EJYBRBQP/30k86cOaM1a9Zo5cqVGjFixAW9b/r06QoODtacOXN0+PDhTPvdbrf++9//nvcc69ev15kzZ/xfQ4YMyfK4pUuXqkWLFmrSpIm2bt2q+Ph4zZo1Sw6HQ+vXr/cf17hx4wznO3PmjEqVKiVJevPNN/X000/rjTfe0KlTp7R8+XLt3btXbdq0UVpamv8c119/vaZMmZLh+lOnTtX111+fYVupUqUyXatx48batGmT//VNN92kCRMm+F8/88wzGeqMj4/XI488ojvuuEPx8fEKDQ3Vxx9/rCeeeEJHjhyRJG3ZskVvvfWWPv74Y9ls/MgAcOVx5wFgCdddd53at2+vjRs3XtDxU6ZM0UMPPaQ6deros88+y7R/8ODBGjNmTIaWy0s1ZMgQ9e/fX8OGDVPx4sUlSWXLltVLL72UobU2O6dPn9YLL7yg8ePH65ZbblFQUJDKly+vr7/+Wnv37tW0adP8xzZo0EBJSUnatGmTJGnTpk1KTk5WgwYNLvtz/JPNZlOfPn2UmJio7du3S5KaNWumnj176tFHH5Vpmrr//vs1bNgwVatWLcevDwAXgrAKwBL279+vmTNnqm7duv967L59+7RgwQL17t1bvXv3zvTYWpLq16+vFi1aaMyYMZdVV2JiopYtW6Zu3bpd8jmWLl2qlJQUde3aNcP2iIgItW/fXnPnzs2wvU+fPv7PNGXKFPXt2/eSr30+Ho9Hn3zyiYKCglSuXDn/9tGjR2vlypXq1q2bUlJSNHjw4Fy5PgBcCMIqgIDq0qWLChYsqKZNm6p58+b+R9XnM3XqVNWpU0c1atRQr169tGnTJq1duzbTcS+//LLGjx+vY8eOZXmemJgYFSxY0P81e/bsTMecPHlSXq9XJUqU8G8bMmSIChYsqPDw8AzdFpYvX57hfJUqVZIkHT9+XEWLFpXD4ch0/pIlS+r48eMZtt1111364osv5HK59OWXX+quu+7K9L5Dhw5luFbBggUz9NE9n/Q6Q0JC9NRTT2natGkqVqyYf39ERITeeecdff/99/r4449lt9sv6LwAkBsIqwAC6ocfflB8fLz27t2rd999V6Ghof/6nqlTp6p3796SfH03mzdvnqmfpyTVqlVLHTt21GuvvZbledasWaP4+Hj/V7t27TIdU6hQIdlstgz9Yl9//XXFx8frtttuk9vt9m+/4YYbMpxv586dkqSiRYvq+PHjGY5Nd/jwYRUtWjTDtrJly6py5cp65plnVKVKFZUpUybT+0qVKpXhWvHx8QoPD8/yc/5Tep0nT55Up06dFBsbm+mYmjVrZvgTAAKFsAogT1m6dKm2b9+uUaNGqUSJEipRooRWrFihL774Issw+NJLL+nDDz/UwYMHL+l64eHhatSokb777rtLrrlx48YKDg7OdI7ExET98ssvatWqVab39O3bV2+++WaudQGQfC2o7777rj799NMsW6YBwAoIqwAsy+12KyUlxf/lcrk0ZcoUtWnTRps3b9a6deu0bt06bdy4UUlJSfrll18ynaNy5crq2bNnpumgLsbrr7+uSZMm6bXXXtPRo0clSQcOHNDu3bsv6P0FChTQCy+8oIEDB2rWrFlyuVzas2ePevToodKlS6tPnz6Z3tOzZ0/NmTNHt99++yXXfSGKFCmi++67Ty+//HKuXgcALhVhFYBlPfzwwwoNDfV/9erVS19//bUGDhzob1UtUaKEKlSooD59+mTZFUCSnn/++Sz7c0ZFRWWYp/Txxx/P8v1NmzbVvHnztGjRIl1//fUqWLCgbrnlFrVo0UIDBw70H7ds2bJMc5+uXLlSkq+f68iRI/XUU08pf/78atSokcqUKaPffvtNwcHBma4ZGhqq1q1bZ9st4tChQ5mu9e233/7r9zQrjz/+uGbOnKkNGzZc0vsBIDcZpmmeb/95dwIAAAA5xMhqIy2rAAAAsCzCKgAAACyLsAoAAADLIqwCAADAsgirAAAAsCzCKgAAACyLsAoAAADLIqwCAADAsgirAAAAsCzCKgAAACyLsAoAAADLIqwCAADAsgirAAAAsCzCKgAAACyLsAoAAADLIqwCAADAsgirAAAAsCzCKgAAACyLsAoAAADLIqwCAADAsgirAAAAsCzCKgAAACyLsAoAAADLIqwCAADAsgirAAAAsCzCKgAAACyLsAoAAADLIqwCAADAsgirAAAAsCzCKgAAACyLsAoAAADLIqwCAADAsgirAAAAsCzCKgAAACyLsAoAAADLIqwCAADAsgirAAAAsCzCKgAAACyLsAoAAADLIqwCAABLME0z0CXAggirAAAgoEzTlNvtVmpqqrxeb6DLgcU4Al0AAAC4NpmmKY/HI7fbLdM0/V/AuQirAADgivpnSDUMQzabjVZVZImwCgAAroj0x/0ejydDSAXOh7AKAAByVXpIdbvdkkRIxUUhrAIAgFyRVUg1DCPAVSGvIawCAIAcRUhFTiKsAgCAHEFIRW4grAIAgMvi9Xr9o/slQipyFmEVAABcMpfLRUhFriKsIk9p3bq1vvnmm0CXAQD4m91ul9PpJKQi1xBWkafExcXJ4eCfLQBYhc1mI6giVzHJGQAAACyLsAoAAADLIqwCAADAsgirAAAAsCzCKgAAACyLsAoAAADLIqwCAADAsgirAAAAsCzCKgAAACyLsAoAAADLIqwCAADAsgirAAAAsCzCKgAAACyLsAoAV9jTTwfr6aeDA10GAOQJjkAXAADXmj/+oJ0AAC4Ud0wACKAZMxy6554QJSUFuhIAsCbCKgAE0PLldv3vfw6Fhga6EgCwJsIqAATQ3r2GypXzyjACXQkAWBNhFQACaO9em8qVMwNdBgBYFmEVAALENNPDqjfQpQCAZRFWAQTEqlU2JScHuorAOnlSOn3aUPnyhFUAyA5hFcAVd+aM1LVrmKKiwjVxYpBSUgJdUWDs2eO7BdMNAACyR1gFcMVFREiff56sSpW8Gjw4RHXrhuujj4KUmhroyq6svXt9t2BaVgEge4RVAAHRtKlHM2cm66efklSmjFdPPBGimJhwTZ4cJJcr0NVdGelhlT6rAJA9wiqAgDEMqXlzj2bPTtb33yepeHFTgwaFqF69cE2b5pDbHegKc9fevYYKFTKVP3+gKwEA6yKsAgg4w5BatfLot9+S9M03SSpUyNQjj4Sqfv1wffGFQx5PoCvMHXv22OgCAAD/grAKwDIMQ2rXzqMFC5L05ZdJiogw9eCDoWrYMEzffHP1hVamrQKAf0dYBWA5hiF16ODRokVJ+uyzZDmd0r33huqGG8L03XcOea+CfGea0r59BjMBAMC/IKwCsCybTbr1VreWLEnSlCnJMgypX79Q3XhjmH78MW+H1qQkQ2lphsLC8vCHAIArgLAKwPJsNum229xatixJH3+crLQ0Q3fdFapmzcI0c6ZdZh5snDx92pAkORwBLgQALI6wCiDPsNulHj3c+v33RE2cmKyEBEN33BGmFi3CNGtW3gqtNpuv2P/7v6t8ygMAuEyEVQB5jsMh9erl1urViXr33WSdOGHo9tvD1LJlmH79NW+E1pQUQ4ZhqlKlPFAsAAQQYRVAnuVwSHfd5daaNYl6++0UHT1qqGvXMLVtG6YFC6wdWlNTpVKlTAUHB7oSALA2wiqAPC8oSOrXz6U1axL11lsp2r/fUKdOYerQIVSLF9sDXV6WUlIMpq0CgAtAWAVw1QgOlu67z6V16xL1xhsp2rnTpg4dwnTrraGaPt1aoTUlRUxbBQAXgLAK4KoTEiI9+KBL69cnatSoFP3xh0333BOm+vXD9Pvvgb/t7dxpU1qaaFkFgAsQ+Ls2AOSS0FBpwACXVq5M1C23uHTsmKHWrcPVvXuo1qwJ3O0vIcGQRDcAALgQhFUAV73ISOnrr1O0aVOiXnwxVStX2tWiRbjuuCNE69df+dtg+mIG5cvTDQAA/o1hnn+4LHdS5J6EBKljR+n48Qt+y85du1S+fPncqwnXhNOecI0/2UdjT/RTvLeAbouYo+eLTlDtkO1X5Pr5/lyrZDNEeyq1UOmgI1fkmkBuMSTfGskXIl8+uT7/XCpdOsvdXq9XTqdTdru1+pjjisnyHxJrpyBwHA6pXDlp8+asA2uNGlL16hlugrsOHlTZatWuYJG4GkVIGqYFejjtd72981a9vfNWfb+nrbpft1jDq36lGvkP5Or1PdtskmmqRK3C8hqFcvVaQG4zJBm2v59QmKaMXbtkbNgg4x+NYWaBAvLefLMUFnbli0SeRssqAs80pa1bpYULpUWLfH8eOuTbV7So1KyZ/6te//5aEBsb2Hpx1TlxQnrnHafee8+pxESpWze3hg1LVZUquXMLLFw4Qh6PdOrUmVw5P3DFmKbsf/6poOXLZYuNlW3xYhmHD/t2RUbK27SpzKZN5b3pJpk1a/qWoTsPWlaveVm2rBJWYT2mKe3a5Quu6eF1925JUoLNptA2beRp0kTuJk3kjY72TbIJ5IC4OEPjxwfp/fedSkmRbr/draefTs3xVaby54+QYRBWkQd5PLJt3Cj74sWyL10q+9KlssXFSZLM666Tt2lTX0C96SaZVateePeAvxFWr3mEVeRh+/dLsbGaPmiQbitSRPZt2yRJZliYPI0aydO4scwiRQJcJK4WRxNCNXZBA01cEqU0j129623W0DYrVKHIqcs+d1KaQ0WGDZQkJb859rLPB+Q22+7dcvz4o8ySJWXbskXG6dOSJG/58vI0aSJv06ZS8+ZShQoXHU7/ibB6zSOsIu+rW7euFi1aJOPoUd9v9UuW+H6z37gxU/8o4HJMVR+V1CHN1P/pPT0sj+zqr0/0rF5VOe27rHPb5ZZkyiOeCiDvMA1Drv795bnxRnluvFHm34OkbDabHI6cGQJDWL3mEVaR96WH1UwSEmSkpFz5gnBVOpNoqEG7MjJNacH3B2Wa0pvvFdQnX+WXaUp9eyToqYdPqnQpzyWdv+D1FSRJ8dt252TZQO5ITJSRlCTv9df7Bsb+A2EVOYiwirwv27AK5LANG2xq2zZM1at7NXNmkkJDpQMHDL35plNTpwbJNKV77nHpySfTVLLkxd0qCxWKkCSdPEmfVeR9hFXkoCzDKosCAEAW6tTx6qOPUrR6tV2PPBIi05RKlzY1dmyqypb1yu2WJk0KUlRUuIYODdaRI5fXVw8AkDXCKgBko2NHt158MVXffhuk1193SpIWLbJr5067JEPvvZei7t3dmjgxSHXqhGv48GAdP05oBYCcRFgFgPP4z3/SdMcdLr36arC++sqhQYNCVL68V3a7qc2bbXr33RStXJmozp3dmjAhSLVrh+uxx4K1fDm3VwDICfRZRZ5Cn1UEQmqq1LFjmFavtsntNvTTT0kaPdqphARDsbFJ/uO2bbPptdecmj7doSJFTO3enZjl+eiziqsJfVaRg+izCgCXIjhY+uijZJmmVLGiR82be9SypUfr19t17NjZe+v113s1aVKKbrnFHcBqAeDqQlgFgAtQrpyp2NhE/fqrryW1ZUtfIJ0/P3ML0M03exQXZ9OhQ/RfBYDLRVgFgAtUq5apokV9f4+K8qpwYa/mzcv8+DMmxjf/6tq1PMoEgMtFWAWAS2C3Sy1aeDRvnl3/7Ppfu7ZvANbq1dxiAeBycScFgEvUqpVbf/1l05YtGW+lYWFSjRperVlDyyoAXC7CKgBcoptv9j3unzcvcyiNifFo7drMra4AgItDWAWAS1S6tKmqVT367bes+q16dfKkod27GWQFAJeDsAoAl6FlS4+WLLErJSXjdgZZAUDOIKwCwEXavNlQ3brhats2VLNmOZSSYmjZsoyhtEYNr4KDTfqtAsBlypklJ3BVWrVqlWJjY3Xo0CGFhoaqVq1aat26tQoXLhzo0oCACg83tXOnoZMnbYqJ8WrvXkPz5jn8fVglKSjINyvAmjW0CQDA5eAuikwmT56smJgYjRo1SsnJyapataqKFSumxYsXq02bNrr77ru1b9++QJcJBEy5clLfvi6lpBj65JNkNW3qyXaQ1bp1dnk8mc9hmtIPPzi0eDEtrwBwPrSsIpPExEQtWbJEoaGhWe5ft26dtm/frrJly17hygDr6N/fpalTnfrqqyC1bOnRiy8G68gRQ8WLnx3+HxPj0QcfOLVtm03Vq3v9233h1VDfvqHq2tWlpk2zSLMAAEm0rCILAwYMyDaoSlJ0dLRatWp1BSsCrCcmxquoKI8++SRIN9+c9dKr9er5AmrWXQF8y7eOH5+SxT4AQDpaVpGt3bt3a/z48dqzZ4/cbrd/+48//hjAqgBrMAzpnntceuyxEKWmSkWK+JZeveOOs/9XqlTxKiLCN8iqd293pnNERXkzbQMAZERYRba6dOmie++9V7feeqtsNhrhgX/q3t2lZ58N1uTJTt1889mlV42/p1a12aToaE+GGQG85FMAuCiEVWQrJCREgwYNCnQZgGXlyyfdfrtLn38epBEjUjR9epD++MMmp9NUtWq+vqsxMV69/36Q0tIkp1N6772gAFcNAHkLYRXZeuyxx/TSSy+pbdu2Cg4O9m+PiYkJYFWAtfTv79KkSU6dOOF7+tCxY6gSEw0dPXpGdrtvkFVamlMbN9r0888OjRkTLIk1WAHgQhFWka0//vhDn376qebNm+fvBmAYhubNmxfgygDriIryql49j7791qFq1Tw6ftyQy2Xo228duv12t38lq0ceCdbmzQ7165emyZNpXQWAC0VYRba+//577dq1S06nM9ClAJZ2zz1pGjAgVLfd5tKPP/puq6NHO9Wtm1tbtthkGKY2b3booYdSNXo0YRUALgajZpCtqKgoxcfHB7oMwPK6dnUrf35TcXGGPB7f6Krt2+1q0iRMPXuGKSTEN1vA6NFp/sFXAIALQ8sqsnXkyBFVq1ZNDRo0yNBnlamrgIzCw6WePV2aMiVIQUGmoqK8+vNPQ9u22TRmTIoOHTL03/86lZwshYWdfd+WLYaqV6f/KgCcD2EV2XrppZcCXQKQZ9xzj0sffuhUpUoeJSdLzz2XpiFDQlSlilelS0sej6ENG+y64Yazq1U1aRKuDRsSVbp04AKr2y39979B6t7dpXLlAlYGAGSLsIpslS1bViVLllRISIgkKTk5WUeOHAlwVYA11azpVaNGHu3YYSguzqYOHZI1dqxXY8Y49fHHvlWqVq+2ZQirbrehNm3CNH9+kkqUuHKB1TSltWttmj49SF9+6dDx4zb9+addH3zAaloArIc+q8hWjx49MiwGYLfb1aNHjwBWBEjJydLu3dbs+Nm/f5ri4nz/Z/77X6cGDUpTbKxDe/bYVKqUV2vXZlyONTjY1JEjhjp1CtXx47n/mbZvNzRypFMxMeFq0SJcH3wQpOLFfSH58cdTc/36AHApCKvIltvtzjATgNPpVFpaWgArwrVu1SqbatYMV8OG4dqxw3qB9bbb3AoP94W/+HhD/fq5VKSIV2+84VTduhlXspKkW291KyxM2rXLppiYcK1alfOf6eBBQ2+/HaRmzcJUr16ERo92qnRpryZMSNH27WdUtapX113npe8sAMsirCJbkZGRGQZTzZgxQ0WLFg1gRbhW7d1rqH//ELVsGS6325DNJrVuHably+3//uYrKDRUcjpNSaaGDUtVeLg0YIBLc+c6VLKkqR07bDp3go2ePV06fdpQ27ZuxccbOnQoZ27JcXHSpElBat8+VDVqhGv48BDZ7dKoUSnaujVRP/2UrL59XSpYUIqNtatJEw+zFACwLMIqsvX+++9r5MiRKlu2rMqWLavRo0frgw8+CHRZuIacPCkNHx6sevXCNXOmQ0OGpGrTpjNavDhRBQv6Vov66itrdb1/5JE0SYY+/ND3VOL++9NUoICpjRt9t9t1684G7JYtPSpSxKs//rApKMhUhw6erE55Qc6ckb7+2qEePUJVpUqEHn88RMeOGXrmmTStXXtG8+cnacAAl0qWPNuCun27TceO2XTTTZd+XQDIbda6y8NSKlWqpOXLl+vMmTMyTVP58uULdEm4RqSlSR99FKTRo4MVHy/17u3W8OGpKlXKF7Ty5TP122+JuuuuUN1/f6h27EjVM88Edg5T05R27jTUvbtbr74qffONQ6NGpapAAemBB9L0xhu+6d/O7QqwcaNNLpe0d69d11/vkeMi78hpadJvv9n19ddB+uUXh5KSDJUu7dWAAS716OFS7dre835PYmN9tTRt6r7ozwsAVwotq8hk2rRp8nq9/tcREREZgurOnTu1ePHiQJSGq5xpSjNmONSwYbiGDg1RdLRHixcn6d13U/xBNV2RItKMGcnq3dul0aODdc89IUoJ4GD2t95yqmHDcOXPL4WHmzp+3KYFC3xh8JFHXAoLMxUebmrNmrO33eBgyemU7HZTlSt7szt1Bh6PL2QOGhSsypUj1LNnmObPt6tXL5dmzUrSxo2JeuWVVNWpc/6gKvnOU7KkVxUr0l8VgHXRsopM4uLiVLduXdWrV0/16tVTZGSkUlJStGPHDi1cuFBFixbVa6+9FugycRU6dUoaODBEJUt69e23SWrd+vx9KZ1O6d13U1Slilcvvhisffts+uKLZBUrduXD1y23uPXSS8GaPt2hRo08WrjQrk8+Cfr7Ub+pe+5x6Z13grRy5dmW1Ro1vFqzJlFly+ZTw4bZh1WvV1q/3qZvvgnSt986dPiwTeHhpjp2dKt7d5datvQo6CJXcDVNackSu5o1o78qAGszTPO8N3V+3b5GeTwezZs3T0uWLNHhw4cVGhqq6tWrq3379ipbtmzA6qpbt64WLVoUsOsj923ebNP113sv+pH4jBkOPfBAiCIjTX39dbJq1Liwlsqc1Lx5mExTatfOrdGjnXI4pM2bE1WihKnDhw3VrOkbIJZ+az19+oxWrLCpTZtwffVVktq3z9x3dOVKm9q2DZPHYygoyFTbtm517+5W+/buDKthXax16ww1axaht95K0X33uS79RLjm2Ww2OS72P2w2vF6vnE6n7HZrDZ7EFZPlr860rCJLdrtdbdq0UZs2bQJdCq4xlxoyO3d2q3TpJN1xR6jatAnTlCnJat36yg4cuvNOlwYPDlHv3r5BVm639OmnQRo8OE0lS5pq186t//0vYxPoli2+H8rVq2f9uXfutMnjMfSf/6Tq8cfTVKhQztT68ce+AWA//OBQr14uhYfnzHkBIKfRZxXAVaNePa/mz09S+fJe9egRqg8/vMhn45epe3eXgoJMbd3qC6AVK3o0ZUqQPH9n5uefT9U/H1ht3ux7pF+2bNYPsubNc6hQIVPPP59zQVWS3n47Vc88k6LFi+3q2DHsiixKAACXgrAK4KpSurSpWbOS1KaNR08+GaIhQ4L9YTG3FSkidejg1o8/OlS2rFdFi5rat8+m335Lbz01VaBAxlC6datNVat6Zcvibuz1SnPn2tW6tVs5/VTUMKShQ1367LMUbdpkU+vWYdq1i8AKwHoIqwCuOvnySV98kawBA9L0/vtO9ewZqoSEK3PtO+906dgxm667zqvDh22KjPTqk0/OtvA2a5YxOW/ebMu268Pq1TbFxdnUrl3uTS31f//n1k8/JSk+Xpoy5cq2RAPAhaDPKrKVmpqqb7/9Vnv27JHbffaH5fPPPx/AqoAL41uxKVWVK3v11FPBats2TF9/nawyZXJ33Gjr1h5FRnp1+rSh/ftt6tLFpR9+cGjhQruaN/eoeXOPfvopSJKpffsMHT1qU/XqWS9jPHu2QzabqVatcnce1EaNvIqNTco0PRgAWAEtq8hW586dNWPGDDkcDoWHh/u/gLzk3ntd+vbbZB04YNPNN4dp1arcve0FBUk9e7r155++65Qu7ZVkaPdu3yP2mJizLav79vmOyW5w1dy5DjVo4FWRIrlasiSpTBkzx7saAEBOoGUV2Tpw4IBmzZoV6DKAy9aypUdz5ybp9ttD1aFDmCZOTNFtt+Vea+Wdd7o0YYJThmHq9GlfSC1a1Levdm2v0gdZbd7sC6tZdQP46y9Da9fa/x6UBQDXLlpWka0bb7xRf/zxR6DLAHJEtWpezZuXpKgor+6+O1RvvOHU+aeZvnS1ankVHe1RSIi0fbvvNnvqlG9fcPDZ47ZssalAAVMlS2YuZO5cXzNn27YshQrg2kbLKjKpXbu2DMOQ2+3WJ598oooVKyo4OFimacowDG3YsCHQJQKXpGhRUz/9lKQBA0L0yivB2rHDprffTskQIHPKnXe6NGRIiDZs8IXO9BbWc23ZYlP16lmvIDVnjkOlSnn/bokFgGsXYRWZ/Pzzz4EBiVG2AAAe3klEQVQuAcg1ISHSRx/5lmgdOTJYe/ca+uyzFBUpkrPNrD16uDVsmKnERF8SzTqs2nXbbZlXj0pL882v2q2bi6VQAVzz6AaATMqVK6dy5cpp+PDh/r+fuw3I63xzjKZp0qRkrV5tV8uWYdq2LWdvh0WKmGra1DeYKijIzDKsnjxpqFq1zC2ny5bZlZBgqF27K7sCFwBYEWEV2dq0aVOG1x6PR6tXrw5QNUDO697drZ9/TlJCgtS6dZgWLMjZ4fAPPOCbkspmk06fzvqYrAZXzZnjkNNpqnlz+qsCAGEVmYwaNUr58uXThg0blD9/fuXLl0/58uVTsWLF1Llz50CXB+SoRo18S7SWLOlV166hmjw55ybGv+UWj4KCTLlcWXcDkLKetmr2bLuaNPEoIiLHSgGAPIuwikyGDRumhIQEDR48WKdPn1ZCQoISEhIUFxenUaNGBbo8IMeVK2dqzpwkNW/u0aBBIbr77hAdP3755w0KkmrW9MrrlY4dyxxWixb1KjIyY1/Z3bsNbdtmz9VVqwAgL2GAFbI1cuRIfffdd1q8eLEMw9BNN92kLl26BLosIFcUKCB9802yHnwwRN98E6Tvvw9ShQpe1a3rUUyMR9HRXkVFeVSgwMWdt1s3l9atC9GOHZnbBrLrAiAxZRUApCOsIlsDBgzQjh071KtXL0nS+++/r7lz5+qdd94JcGVA7nA4pPHjU1ShgldOp7Rhg02rVtn13XdnuwZUruybQzUmxqO6db2qU8ejfPkk01SmkftHjkjff+977/HjmVtWsxpcNXu2Q5UqeVW5MkufAoBEWMV5LFy4UBs3bpTx90/gu+++W7Vr1w5wVUDuCguThg9Py7Dt+HFDa9fatHatXevW2bRsmV3Tp/tCqGGYCgmR+vVzafTojKtN5c8vrVljk2TK7Ta0ebPtnNZUI1PLamKiFBtr1733Zp7OCgCuVYRVZKtq1arat2+fypUrJ0nav3+/6tSpE+CqgCuvaFFTbdp41KbN2amkjhzxBdg1a+waO9apJUsyzyQQGuqbwiouzpBk6vPPgzRixNlA+8/BVYsW2ZWaatAFAADOwQArZCsuLk7Vq1dXixYt1KJFC9WoUUPHjh1Tp06d1KlTp0CXBwRU8eKmbrnFo2eeSdNjj6Xpjz9s2rcv86P+unW9kgxJhj7/3CH3OTm0WrWM86jOmeNQeLipJk2YXxUA0tGyimy9/PLLgS4ByBP69nXpjTec+vTTID37bMYuBB07uvTrr75b7fHjNv32W3oLrKlChc4eZ5q+/qo33+zOleVfASCvIqwiW82bN9fevXu1fft2tW7dWsnJyXK73cqXL1+gSwMspWxZU61aeTRtWpCGDk2T/ZweATfeePZRf2iorytAVrZssenAAZuGDEnLcj8AXKvoBoBsffjhh+revbsefPBBSdKBAweYugrIRr9+Lh08aNPcuRn7rlap4pXD4RvZny+fV//7X9ZtBLNnM2UVAGSFsIpsvfPOO1qyZIny588vSapSpYqOHj0a4KoAa2rf3q1ixbyaMiVjy6ndLlWo4GtdPXXKUFpa1itZzZ5tV506HpUqxZRVAHAuwiqyFRwcLKfT6X/tdrv901gByCgoSOrd26VZsxz666+M/0+io31hNTXVpurVMw+eOnlSWrHCTqsqAGSBsIpsNW/eXCNHjlRycrLmzp2rHj166NZbbw10WYBl9e3rksdjaNq0jK2r0dFnA2rNmh75Zgc4a948hzwegyVWASALhFVk67XXXlNkZKRq166tiRMnqkOHDhoxYkSgywIsq1IlU82auTV1apC850yh2qjR2bDqW8kq46P+2bMdKlTIVP36mVe0AoBrHbMBIFs2m01dunRRly5dFBkZGehygDzh7rtduvfeUC1YYFfLlr6QGhWVHkJN/fHH2QFYCxfa1bSpR3Pn2tW6tTvDLAIAAB9aVpGJaZp68cUXVbRoUVWrVk1Vq1ZVZGQk864CF+DWW90qVMjMMNAqOFhyOEyFh6evZuXTuXOoJk50KC7ORhcAAMgGYRWZjBs3TkuWLNHKlSsVFxenEydOaMWKFVqyZInGjh0b6PIASwsJkXr1cunnnx1/P/L3KVLEVLlyps7tr+r1Gtq0ySabzVSrVoRVAMgKYRWZTJ06VV988YUqVKjg31axYkVNmzZNU6dODWBlQN7Qr59LLpdvedV0+fObKlEi87RU69c71KCBV0WKXMkKASDvIKwiE5fLpaJFi2baHhkZKZfLFYCKgLwlIcE3A8DkyUEyTd9Sqm63oUOHDP8CAZIUEmJqwwY7XQAA4DwYYIVMzp1b9WL2AZCSkqRWrcL9r8uVi1Bamq/PamqqTZUre7V1q12SqdBQUykpTFkFAOdDWEUm69ev969adS7TNJWSkhKAioC8w+GQpk9P0s6dNj33XLBCQ03FxxuqUcPU5s02NWzo1tatvodaLpehokW9qlWLKasAIDuEVWTi8WReYQfAhXE6pbZtPZI82rbNps8+C9ITT6TprbeCJUlhYYZ8g6xMnTkjnTlj0223heqHH5IDWTYAWBZ9VgEgl/Tv71JKiqHISFNRUR4Zhql9+85dvcqQYZjq35++4ACQHcIqAOSSOnW8qlvXo2nTgvT++ykyDGn+fIfOXcHqwQfT1LkzfVYBIDuEVQDIRf36ubRpk11JSVL79m4lJ6e3rPq6Ajz/fFogywMAyyOsAkAu6t7dpfBw34pWvXplfNxfsKCpiIgAFQYAeQQDrAAgF+XLJ3Xr5tKXXwbpq6+CZBimzL97AVSvzmBGAPg3tKwCQC5KSpKOHjWUlmaocGFTjnOaCGJimLIKAP4NYRUAcsmOHYZatw7T7NkORUZ6lS+fqRo1PPL1V5W+/DJIkybxgAsAzoewCgC54IcfHGrePFyHDhmaPj1Zgwenads2u4YPPzugKi7OppdfDg5glQBgfYRVAMhBaWnS008Hq2/fUFWr5tXixUlq08ajnj1dCgkxNWtW+tRVpr79NkmTJ7MYAACcD2EVAHLQBx8E6b33nHr44TT98kuSSpf2jaYqVEjq3Nmtb74J8h+bP7+pFi3otwoA50NYBYAc9MADLn33XZJGj06V05lxX79+Lp0+fXae1R49Qq94fQCQ19CzH3mbxyPbxo2yr1ghIykp0NUAckrqIEkbM273eqXfY5tIuvnvLaaeqj9PznHLrmh9wMUydu+WY8YMpUycKE/LllJQ0L+/CchBhmma59t/3p3AldYgOloL3nxTjqVLZV+8WPbly2WcPh3osoDzOq4i6qNPNUvt5VCa3HJKMmXycAt5jBkWJk+jRvLceKM8TZvKU6+ebGFhcjhypu3L6/XK6XTKbrfnyPmQ5xhZbiSswtJSUqQVK6RFi6SFC5U8f75Cvb4+fp6qVeVp0sR307zxRpmFCwe4WCCz5b8Hqd+DBRR3wqbRrySo2Y2pimlaVJJ0aOcxRYRzm4X1GUePyr5mjezLlsm+dKlsGzfKME2ZTqe89evLbNZM3qZNZd5wgy5nWTbC6jWPsIo84MwZadkyaeFCX0BdscI3vNowpDp19PnBg7pt7FhfOI2MDHS1wHm9916QnnkmWGXLmpo6NVlRUb5ftPLn9/0wP3jwjPLlC2SFwCU6cUL2FSvkWLJE9iVLZFu3TobHI9PhkBkT4wuuN90kb+PGUsGCF3xawuo1j7AKC4qPlxYv9recavVqyfP3EpQ1a0rNmknNm0tNmkiFCqlx48aaM2dOYGsGLtCnX4Rq7rxgjX/zlArk53aKq5PNZpMjOVm2FStkxMbKtnixjPXrfS2vhiGzTh15b7pJZtOm8jZpIp2noYGwes0jrMJijh+XKlTwtaYCV6H0G2iWd1/gGmQ6nXKtXi2zSpUs9xNWr3lZ3i6ZDQCBU6iQ9PbbvtB6gcaOG6eHH344F4sCAFwMm2HIZruwwYJmvnwyy5TJ5YpwtaFlFXlK3bp1tWjRokCXAQD4m81mYzYA5JQsW1aZNwUAAACWRVgFAACAZRFWAQAAYFmEVQAAAFgWYRUAAACWRVgFAACAZRFWAQAAYFmEVQAAAFgWYRUAAACWRVgFAACAZRFWAQAAYFmEVQAAAFgWYRUAAACWRVgFAACAZRFWAQAAYFmEVQAAAFgWYRUAAACWRVgFAACAZRFWAQAAYFmEVQAAAFgWYRUAAACWRVgFAACAZRFWAQAAYFmEVQAAAFgWYRUAAACWRVgFAACAZRFWAQAAYFmEVQAAAFgWYRUAAACWRVgFAACAZRFWAQAAYFmEVQAAAFgWYRUAAACWRVgFAACAZRFWAQAAYFmEVQAAAFgWYRUAAACWRVgFAACAZTkCXQDwb0zT1KpVqzRixAglJibK7XYHuiQAwN/sdrtsNpsMw5BhGIEuB1chwiosyzRNLVmyRK+++qokafDgwXr11Vfl9XoDXBkA4FzpYZXQitxgmKZ5vv3n3QlcDtM0s7yhmaap+fPn69VXX1VERISGDBmi6Ohobn4AYFHpWSL9vn5ucL0YXq9XTqdTdrs9N8qE9WX5D4aWVVxRXq9X7733nvr06aP8+fNn2jd79my99tprioyM1KhRo1SzZk1CKgBYXPp92jAMmaYpr9frD66XElqBc9Gyiivu+eeflyRVqVJFffr0kdvt1s8//6zXX39dFSpU0JNPPqlq1aoFuEoAwOUwTfOiW1ppWb3m0bKKwPB4PBluPM8//7z++OMP9ezZU0lJSfr4449Vs2ZNTZw4UZUqVQpgpQCAnJIeUtNbWr1er2w2Gy2tuGiEVeSa9BtTelBdu3atatWqJcMwtHHjRp08eVLjx4/XF198QUgFgKsUoRWXi7CKHHfujUiSNm7cqIcfflgJCQkqUqSIjhw5opYtW2rGjBkaMmSIzpw5E+CKAQC5jdCKS0VYRY5LD6nvvvuuDh8+rIIFC6p+/fr69ddftWfPHg0cOFC9e/dWeHi4rr/+en355ZeKiory37gAAFev84VWICuEVVy2f/ZJnTNnjhYsWKB9+/bpxIkTmjVrlnr16qX//e9/+uGHH7R582bt3r1bNWrU0EMPPaTHH39cLpdLdrudOVQB4BpjmqZ/sRfDMOR0OgNcEayG2QBwyc6dJ9U0Tc2YMUNdunTRgAEDtGbNGo0ePVr33nuvTp8+LbfbrfLly8vtduvgwYOKiIhQoUKF5HK5FBcXp2LFisnh4HcnALhWpeeRYsWKadasWQGuBgHCbAC4fKdPn9bIkSP12muv+YPqxIkTNX78eNWvX1/t2rXTXXfdpePHj0uStm/frn379um2227T1KlTVbNmTb3xxhtKSkrSs88+q+3bt2vq1KkaNWpUID8WAACwKFpWcdF69Oihl156STVq1NDBgwf12GOP6fXXX1fFihX9xwwaNEjly5dX//79VahQIQ0cOFAHDx7Ud999p9TUVAUHB/uP3bNnj8qXLx+ATwIAACwky5ZVejPjvFJSUiSdfTxz8OBBVa5c2d+qmpSUpB07dqhIkSIZju/YsaM2bdqkjRs3SpKefvppDRgwQJIUFBQkSf7+qQRVAACQHcIqsuTxePTss8/q448/lnR2Kb3rrrtOK1eu9D/mT0hIUOvWrbV69WpJUkhIiNLS0tS2bVuZpqkDBw7I7XardOnSatWqlaSzswUw8hMAAPwb0gKylJycrAULFqhbt27+bR6PR5LUrFkzvfLKK5KkChUqqEyZMnr11Ve1detWffHFF+rUqZO2bNmicePGqVevXgycAgAAl4wUAUm+x/ndu3dXo0aNdP/996tUqVIqU6aMDhw4oBIlSkg62xL67LPPKjIyUuvXr1dUVJQeeeQRud1ujRgxQn/99ZdeeOEFVa9e3X/uc2cNAAAAuBgMsIJ/QuYtW7Zo0qRJOnTokMqVK6dt27Zp+vTpGY51u91yOBx65ZVX9Oeff2r06NG67rrrJEmJiYkKDw8PxEcAAAB5X5YtW4TVa1R2rZ0ej0crV67UsGHDFBcXp++++06VK1f2H5/+p8vl0ogRI3Tq1Cnde++9ql279r+eGwAA4DwIq/C1oqYvdZeV9KC5atUqff/994qKitLtt9+e5TFer1eff/65Tpw4oQEDBmRYxQoAAOAiEVZx1pIlS7RlyxY1bdpU1apVy/KYt956SydOnND999+vcuXKZdtimpycrNDQ0NwuGQAAXN2YZ/ValD6CX/K1iJ4+fVqPP/64hg0bptTUVPXr188/7dS5x0lS06ZNdfDgQS1dulSSsm2NJagCAIDcQli9yqU/mt+wYYMMw1BiYqIKFSqkWbNmqWzZsvrzzz+1fft2/wT90tlQ2rBhQ/Xu3VsdO3YMSO0AAAB0A7jK/fbbbxozZozi4+P1448/as2aNXrjjTd0+PBhNWjQQEOHDlW1atUYFAUAAAItyyDCPKtXifTpp86VkpKiyZMnq1+/furZs6ckqXr16vJ4PBo5cqQ6d+4sSZo3b54kqWXLlle2aAAAgH9BWL1KpAfVX375RYUKFdINN9ygTZs2acWKFfr0008lSampqSpbtqzatWunr7/+Whs2bND+/fu1dOlSPf3004EsHwAAIEv0Wb1KzJ07V7169dKIESM0aNAgzZ49W1WqVFGVKlUUGxsrSQoODpYkPfbYY3riiSfkcrlUqVIlrV69Wn369Alk+UCusNvtio6OVs2aNRUVFaW33nrL3z971apVGjRoULbv3bNnjz7//PMrVWoG+/fvV4UKFXTixAlJ0smTJ1WhQgXt3bs307HJyclq3rx5hsGUlyItLU3NmjWT2+2+rPMAQE4jrOYxWf1AWrp0qR588EHddNNNWrJkie68807NnTtXW7duVdeuXfXss88qLS1Na9euVc+ePbVv3z7Vq1dPL7/8sp5++mkFBwdnGGAFXC1CQ0O1bt06bdq0SXPnztXMmTP10ksvSZLq16+vt99+O9v3BjKslilTRg8//LCGDh0qSRo6dKgeeOABlStXLtOxkyZNUteuXS97nmOn06lWrVrpq6++uqzzAEBOI6zmEekD4dJ/IE2YMEG//PKLTpw4ocaNG6tWrVr6888/JUndunVTYmKitm7dqnvuuUeVKlXSXXfdpfvuu0+NGzdW1apV/edL//Of/V2Bq02xYsX0wQcfaMKECTJNUwsWLPDPdLFw4UJFR0crOjpadevWVUJCgoYOHarY2FhFR0dr7Nix2rNnj2666SbFxMQoJibGP6XbggUL1KJFC3Xv3l3VqlVT7969/f+vVq5cqRtvvFFRUVFq2LChEhIS5PF4NHjwYDVo0EB16tTRxIkTs6z3P//5j5YvX65x48Zp8eLFevLJJ7M87rPPPvP3P5ek119/XbVr11ZUVJQ/7LZo0UL/+c9/1KxZM1WvXl0rV65U165dVaVKFQ0fPtz/3i5duuizzz67/G82AOQk0zTP94UASk5ONpctW+Z/7fV6zZ9//tmsXbu2+cADD5hjxowx27Zta5qmac6cOdO89dZbzcOHD5umaZrvvfee+cADD5hbt241TdP0bweuJeHh4Zm2FSxY0Pzrr7/M+fPnm//3f/9nmqZpduzY0Vy8eLFpmqaZkJBgulyuDPtN0zQTExPN5ORk0zRNc9u2bWa9evVM0zTN+fPnm/nz5zf3799vejwe84YbbjBjY2PN1NRUs0KFCubvv/9umqZpnjp1ynS5XObEiRPNV155xTRN00xJSTHr1atn7tq1K8v6Z82aZUoy58yZk+X+1NRUs3jx4v7XM2fONBs3bmwmJiaapmmacXFxpmmaZvPmzc0hQ4aYpmma48aNM0uWLGkeOnTITElJMa+77jrz+PHjpmmaptvtNosWLfqv31cAyCVZ5lEGWFmQ+fc0UqmpqZo1a5Y2bdqksLAwFS9eXGXKlNGMGTPkdDr1/PPPa+7cuZo3b57at2+vn376SVOnTtWQIUPUrl07HThwQEFBQZKkEiVKSMp61gDgWmJmMV1fkyZN9MQTT6h3797q2rWrSpcunekYl8ulRx99VOvWrZPdbte2bdv8+xo2bOh/T3R0tPbs2aMCBQqoZMmSatCggSQpf/78kqQ5c+Zow4YNmj59uiTp1KlT2r59uypUqJDpmr/88otKliypjRs3qk2bNpn2Hz9+XAULFvS//vXXX9W/f3+FhYVJkgoXLuzf16lTJ0lS7dq1VbNmTZUsWVKSVLFiRe3fv19FihSR3W6X0+lUQkKC8uXLd75vIwBcMaQWC/F6vfJ4PP75ToODg7Vu3To99dRTmjdvnqKjo1WnTh1t375dHTt2VKdOnTR69Gg99dRTkqRWrVppxowZio+PV4UKFTRixAhVrFgxwzUIqriW7dq1S3a7XcWKFcuwfejQofroo4+UnJysG264QVu3bs303rFjx6p48eJav369Vq1apbS0NP++9MGLkq+rjtvtznbuYtM0NX78eK1bt07r1q3T7t271bZt20zHrVu3TnPnztXy5cs1duxYHT58ONMxoaGhSklJyXDu7OZLTq/RZrNlqNdms2UYVJWamqqQkJAszwEAgUBysYBdu3ZJ8v3QsNvtOnXqlIYPH66EhAQ1b95czZo1U48ePVS4cGF5vV5t27ZNt99+uzp37qyoqCitW7dOsbGx6ty5s9555x0VLFjQP2Aqq1Yk4Fp07NgxPfTQQ3r00UczBbqdO3eqdu3aevrpp1W/fn1t3bpV+fLlU0JCgv+YU6dOqWTJkrLZbPr000//dfR9tWrVdOjQIa1cuVKSlJCQILfbrXbt2um9996Ty+WSJG3btk2JiYkZ3muaph5++GGNGzdOZcuW1eDBg/2/lJ6rUKFC8ng8/sDatm1bTZo0SUlJSZLkn03gQsXFxSkyMtL/RAYArICwGkBer1eLFi3SN9984w+V7777rm688Ub99ddf8nq9GjRokO644w599NFHcrlcstls+uuvv3Ts2DGNGTNGs2fP1kMPPaSQkBA5HA5FR0dLOtuCyqpUuJYlJyf7p65q3bq12rZtqxdeeCHTcePGjVOtWrUUFRWl0NBQtW/fXnXq1JHD4VBUVJTGjh2rRx55RFOmTNENN9ygbdu2KTw8/LzXdjqd+uqrrzRw4EBFRUWpTZs2SklJ0X333acaNWooJiZGtWrV0oMPPphpuqgPP/xQZcuW9T/6f+SRR7R161YtXLgw03Xatm2rxYsXS5JuueUWderUSfXr11d0dLTGjBlzUd+v+fPnq0OHDhf1HgDIbSy3GgB79uxR+fLl/a9TU1O1adMmxcTEaNCgQbr//vtVu3ZtSb5+cnv37tVzzz2nvn37qn379oqNjdXBgwc1adIkDR8+XM2aNQvQJwEQaGvXrtVbb73lX/zjcnTt2lWjRo1S1apVc6AyALhoWbawEVavoE8++UTjx49X6dKlVaBAAbVu3Vp33323PvvsM73//vuKjY3V/fffr127dum6667TqVOnFBERoRdffFGrV6/WyJEjVaBAAfXp00cPPPBAhnOfr68agKvbpEmTdPfdd1/WXKtpaWn68ssv1bdv3xysDAAuCmE1UH7//Xe1a9dO3bt3V79+/VS7dm1Nnz5dgwcP1i+//KKKFSvqySef1K233qpu3bpp5cqVKl68uBITEzV8+HANHDhQN998sz799FNVrVpVDRs29J+b0f0AAOAqQVgNlPj4eBUuXFjJyckKDg6Wy+VSUFCQhg8frpUrV2r27NmaNm2avv/+e3399dey2+2Kj4/Xzz//rAkTJmjs2LFq3LhxoD8GAABAbiKsBtJ9992n5ORkffbZZ0pLS5PT6ZTL5VJkZKSWLl2qkiVL6rHHHlOHDh3Upk0b9e7dWwULFtRzzz2nmjVrBrp8AACA3EZYDaSjR4+qcuXK2rFjh4oVK6aUlBSFhISoe/fuatGihR599FG9//77mjNnjr777jsdO3ZMkZGRkuiPCgAArgmE1UC766675HQ6NWnSJEm+WQDatWun0aNHq1GjRjp69KhsNpuKFi0q6ewcqQRVAABwDSCsBtqBAwdUu3ZtnTx5UmfOnNFzzz2nEydOaMKECRmWNqQlFQAAXIOyDD8MI7+CSpcurZYtW6pw4cJq166dJOmtt97KtAY3QRUAAMCHltUr7MCBA5owYYKGDRumAgUKSKIlFQAAQHQDsB7mSAUAAPCjG4CVmKZJUAUAAPgXpKUA4bE/AADAvyOsAgAAwLIIqwAAALAswioAAAAsi7AKAAAAyyKsAgAAwLIIqwAAALAswioAAAAsi7AKAAAAyyKsAgAAwLIIqwAAALAswioAAAAsi7AKAAAAyyKsAgAAwLIIqwAAALAswioAAAAsi7AKAAAAyyKsAgAAwLIIqwAAALAswioAAAAsi7AKAAAAyyKsAgAAwLIIqwAAALAswioAAAAsi7AKAAAAyyKsAgAAwLIIqwAAALAswioAAAAsi7AKAAAAyyKsAgAAwLIIqwCA/2/Hjk0AgGEgiBHw/jNngUA6c4XUfuPyMECWWAUAIEusAgCQJVYBAMgSqwAAZM1nPytXAADAg88qAABZYhUAgCyxCgBAllgFACBLrAIAkCVWAQDIun1K4jGTnDRXAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Run the actual simulation to generate results \n",
    "if op2=='t':\n",
    "    print(\"\")\n",
    "    print(\"Plotting tracks of photons\")\n",
    "    \n",
    "# Data Input\n",
    "if op2 != 't':\n",
    "    print(\"\")\n",
    "    print(\"Computing ....\")\n",
    "    print(\"\")\n",
    "    \n",
    "tic = time.time()\n",
    "#\n",
    "NS=0;NT=0;NR=0;NA=0;NM=0;Emin=0.01;  #fix var name syntax error\n",
    "\n",
    "with open(\"data.txt\") as finp:\n",
    "    intern = finp.readlines()[10].strip().split()\n",
    "for i in range(4):\n",
    "    intern[i]=int(intern[i],10)\n",
    "\n",
    "Mz=intern[0]\n",
    "ME=intern[1]\n",
    "P=intern[2]\n",
    "B=intern[3]\n",
    "\n",
    "A1 = np.array(A)\n",
    "pmol=A1.dot(np.transpose(nam))\n",
    "N=(rho/pmol)*6.022e23\n",
    "\n",
    "if geo=='p':\n",
    "    capaz=np.zeros([1,ME]) \n",
    "    Edep=capaz\n",
    "    Et=capaz\n",
    "    Er=capaz\n",
    "else:\n",
    "    capa=np.zeros([Mz,P])\n",
    "    Edep=np.zeros([1,ME])\n",
    "    Et=Edep\n",
    "    Er=Edep\n",
    "\n",
    "# \n",
    "numtht=np.zeros([1,B]) \n",
    "numthr=numtht\n",
    "time.sleep(2)\n",
    "\n",
    "if op1=='g':\n",
    "    if geo=='p': # \n",
    "        NT,NS,NR,NA,nfot = monfgp.monfgp(H,nfot,d,opt,op2, Rmin, rf, geoh, E0, ZA, N, nam, mat ,pmol, faZ, b0, b1, b2, NS, Emin,NT,NR,NM,NA)  # PLANE PARALLEL LAYER -- GRÁFICAL MODE /P\n",
    "    else:\n",
    "        NT,NS,NR,NA,nfot = monfgc.monfgc(H,nfot,d,opt,op2, Rmin, rf, geoh, E0, ZA, N, nam, mat ,pmol, faZ, b0, b1, b2, NS, Emin,NT,NR,NM,NA,R)  # CYLINDER -- GRÁFICAL MODE /C\n",
    "else:\n",
    "    if geo=='p':\n",
    "        NT,NS,NR,NA,nfot = monfsgp.monfsgp(H,nfot,d,opt,op2, Rmin, rf, geoh, E0, ZA, N, nam, mat ,pmol, faZ, b0, b1, b2, NS, Emin,NT,NR,NM,NA,ME,Edep,Et,B,numtht, Er, numthr)  # PLANE PARALLEL LAYER -- ONLY CALCULATIONS /P # Make function for this\n",
    "    else:\n",
    "        NT,NS,NR,NA,nfot = monfsgc.monfsgc(H,nfot,d,opt,op2, Rmin, rf, geoh, E0, ZA, N, nam, mat ,pmol, faZ, b0, b1, b2, NS, Emin,NT,NR,NM,NA,R)  # CYLINDER -- ONLY CALCULATIONS /C  # Make function for this\n",
    "\n",
    "# RESULTS DISPLAY\n",
    "\n",
    "toc=time.time()\n",
    "T=toc-tic\n",
    "\n",
    "print(\"Simulation is over\")\n",
    "print(\"\\nAverage values:\")\n",
    "print(\"\\n\")\n",
    "if geo=='p':\n",
    "    cftr=(NT+NS)/nfot\n",
    "    cfref=NR/nfot\n",
    "    cfabs=NA/nfot\n",
    "    print(\"Transmission coefficient: \",cftr)\n",
    "    print(\"Backscattering coefficient: \",cfref)\n",
    "    print(\"absorption coefficient: \",cfabs)\n",
    "    print(\"\")\n",
    "    print(\"Fraction of photons leaving the medium without interacting: \",NS/nfot)\n",
    "    print(\"Fraction of photons which suffer more than one interaction: \",NM/nfot)\n",
    "else:\n",
    "    cfabs=NA/nfot\n",
    "    print(\"absorption coefficient: \",cfabs)\n",
    "    print(\"\")\n",
    "    print(\"Fraction of photons which suffer more than one interaction: \",NM/nfot)\n",
    "\n",
    "if op1=='sg':\n",
    "    print(\"Computing time / input photon(mseg) = \",T*1000/nfot)\n",
    "    print(\"\")\n",
    "\n",
    "# Recording of files and plots\n",
    "\n",
    "print(\" \")\n",
    "print(\"FILES AND PLOTS\")\n",
    "print(\"  \")\n",
    "#\n",
    "dz=H/Mz \n",
    "dE=(E0-Emin)/ME \n",
    "dA=90/B \n",
    "gr=np.pi/180\n",
    "\n",
    "Mint=math.floor(Mz/5)\n",
    "m=[1,Mint,2*Mint,3*Mint,4*Mint,5*Mint] \n",
    "color=['k' 'b' 'r' 'g' 'c' 'y']\n",
    "\n",
    "if op1=='sg': \n",
    "     # Energy transfered to the medium\n",
    "     # File: 'Eabs.out' (Transferred energy per unit depth in KeV cm^(-1).\n",
    "     if geo=='c':\n",
    "         with open('Dosis.txt','w') as fdosis:\n",
    "             print(\"File Dosis.txt Transferred energy vs radius (cm) and depth (cm)\")\n",
    "             xxd1=np.linspace(0.5*dz,H-0.5*dz,Mz)\n",
    "             #\n",
    "             dR=R/P\n",
    "             yyd1 = np.empty((Mz, P))\n",
    "             xxd2=np.linspace(0.5*dR,R-0.5*dR,P)\n",
    "             for k in range(P):\n",
    "                 for r in range(Mz):\n",
    "                     yyd1[r][k]=(capa[r][k]*1000/(2*np.pi*xxd2[k]*dR*dz))/nfot\n",
    "             fdosis.write(\"          \")\n",
    "             fdosis.write(\"\\n\")\n",
    "             for k in range(Mz):\n",
    "                 fdosis.write('%.4f ',xxd1(k))\n",
    "                 for r in range(P): \n",
    "                    fdosis.write(\"%.4f  \" % xxd2(r))\n",
    "                    fdosis.write('%.6f  ' % yyd1(k,r))\n",
    "                 fdosis.write(\"\\n\")\n",
    "             \n",
    "     # Graph: 'transferred energy vs radius'\n",
    "             plt.figure()\n",
    "             ax5 = plt.axes()\n",
    "             ax5.set_xlim(0,R)\n",
    "             print(\"Plot Dosis.emf transferred energy vs radius at several depths\")\n",
    "             print(\"colour depth (cm)\")\n",
    "             for k in range(len(m)):\n",
    "                 xx = np.repeat(xxd2-0.5*dR, 2)[:-1]\n",
    "                 yy = np.repeat(yyd1[m[k],:], 2)[1:]\n",
    "                 plt.plot(xx,yy,color(k))\n",
    "                 #\t    plot([xxd2(P)-0.5*dR,R],[yyd1(m(k),P),yyd1(m(k),P)],color(k))\n",
    "             print(\"%s        %5.3f \\n\",color(k),xxd1(m(k)))\n",
    "             ax5.set_title(\"transferred energy vs radius at several depths\")\n",
    "             ax5.set_xlabel(\"Radial distance (cm)\")\n",
    "             ax5.set_ylabel(\"de/dV (KeV/(phton·cm^3))\")\n",
    "             \n",
    "     # transferred energy spectrum     \n",
    "     # File: 'eabs.out'\n",
    "         with open(\"Eabs.txt\",\"w\") as feabs:\n",
    "             print(\" \")\n",
    "             print(\"File Eabs.txt Spectrum of transferred energy (MeV)\")\n",
    "             xxea=np.linspace(0.5*dE+Emin,E0-0.5*dE,ME)\n",
    "             yyea=(Edep/dE)/nfot\n",
    "             feabs.write('#.4f \\t #.8f \\n',[xxea,yyea])\n",
    "        \n",
    "     # Graph: 'eabs'   \n",
    "         plt.figure()\n",
    "         ax0 = plt.axes()\n",
    "         ax0.set_xlim(Emin,E0)\n",
    "         xx = np.repeat(xxea-0.5*dE, 2)[:-1]\n",
    "         yy = np.repeat(yyea, 2)[1:]        \n",
    "         plt.plot(xx,yy,\"b\")\n",
    "         plt.plot([xxea(ME)-0.5*dE,E0],[yyea(ME),yyea(ME)],\"b\")\n",
    "         ax0.set_title(\"Spectrum of transferred energy\")\n",
    "         ax0.set_xlabel(\"Energy (MeV)\")\n",
    "         ax0.set_ylabel(\"dn/dE (MeV^ (-1))\")\n",
    "         print(\"Plot Eabs.emf Spectrum of the transferred energy\") \n",
    "         plt.savefig('eabs.png')\n",
    "        \n",
    "     if geo=='p':\n",
    "             \n",
    "     # Energy distribution \n",
    "     # File: 'etrans.out'\n",
    "         with open(\"Etrans.txt\",\"w\") as fetrans:\n",
    "             xxet=np.linspace(0.5*dE+Emin,E0-0.5*dE,ME)\n",
    "             yyet=(Et/dE)/nfot\n",
    "             for i in range(len(xxet)):\n",
    "                 xxet1 = xxet[i]\n",
    "                 yyet1 = yyet[0][i]\n",
    "                 fetrans.write(\"%f\\t%f\\n\" % (xxet1,yyet1))\n",
    "             print(\"File Etrans.txt Energy spectrum (MeV) of transmitted photons\")\n",
    "        \n",
    "     # Plot: 'etrans'\n",
    "         plt.figure()\n",
    "         ax1 = plt.axes()\n",
    "         ax1.set_xlim(Emin,E0)\n",
    "         xx = np.repeat(xxet-0.5*dE, 2)[:-1]\n",
    "         yy = np.repeat(yyet, 2)[1:]\n",
    "         plt.plot(xx,yy,\"b\")\n",
    "         plt.plot([xxet[ME-1]-0.5*dE,E0],[yyet[0][ME-1],yyet[0][ME-1]],\"b\")   \n",
    "         ax1.set_title(\"Energy distribution of trasmitted photons\")\n",
    "         ax1.set_xlabel(\"Energy (MeV)\")\n",
    "         ax1.set_ylabel(\"dn/dE (MeV^ (-1))\")\n",
    "         print(\"Plot Etrans Energy spectrum (MeV) of transmitted photons\") \n",
    "         plt.savefig('etrans')\n",
    "     \n",
    "     # Energy distribution of backscattered photons \n",
    "     # File: 'erefle.out'\n",
    "         with open(\"Eretro.txt\",\"w\") as ferefle:\n",
    "             xxer=np.linspace(0.5*dE+Emin,E0-0.5*dE,ME)\n",
    "             yyer=(Er/dE)/nfot\n",
    "             for i in range(len(xxer)):\n",
    "                 xxer1 = xxer[i]\n",
    "                 yyer1 = yyer[0][i]\n",
    "                 ferefle.write(\"%f\\t%f\\n\" % (xxer1,yyer1))\n",
    "             ferefle.write(\" \")\n",
    "             print(\"File Eretro.out Energy distribution (MeV) of backscattered photons\")  \n",
    "        \n",
    "     # Graph: 'Eretro'  \n",
    "         plt.figure()       \n",
    "         ax2 = plt.axes()\n",
    "         ax2.set_xlim(Emin,E0)\n",
    "         xx = np.repeat(xxer-0.5*dE, 2)[:-1]\n",
    "         yy = np.repeat(yyer, 2)[1:]\n",
    "         plt.plot(xx,yy,\"b\")\n",
    "         plt.plot([xxer[ME-1]-0.5*dE,E0],[yyer[0][ME-1],yyer[0][ME-1]],\"b\")\n",
    "         ax2.set_title(\"Energy distribution of backscattered photons\")\n",
    "         ax2.set_xlabel(\"Energía (MeV)\")\n",
    "         ax2.set_ylabel(\"dn/dE (MeV^ (-1))\")\n",
    "         print(\"Grafica Eretro Energy distribution of backscattered photons\")\n",
    "         plt.savefig('Eretro.png')\n",
    "     \n",
    "     # Angular distribution of transmited photons \n",
    "     # File: 'atrans.out'\n",
    "         with open(\"Atrans.out\",\"w\") as fatrans:\n",
    "             xxat=np.linspace(0.5*dA,90-0.5*dA,B)\n",
    "             yyat=(numtht/(dA*gr*2*np.pi*np.sin(xxat*gr)))/nfot           \n",
    "             for i in range(len(xxat)):\n",
    "                 xxat1 = xxat[i]\n",
    "                 yyat1 = yyat[0][i]\n",
    "                 fatrans.write(\"%f\\t%f\\n\" % (xxat1,yyat1))\n",
    "             fatrans.write(\" \")\n",
    "             print(\"File Atrans.txt Angle distribution of transmited photons\")   \n",
    "         \n",
    "     # Graph: 'atrans'\n",
    "         plt.figure()\n",
    "         ax3 = plt.axes()\n",
    "         ax3.set_xlim(0,90)\n",
    "         xx = np.repeat(xxat-0.5*dA, 2)[:-1]\n",
    "         yy = np.repeat(yyat, 2)[1:]       \n",
    "         plt.plot(xx,yy,\"b\")\n",
    "         plt.plot([xxat[B-1]-0.5*dA,90],[yyat[0][B-1],yyat[0][B-1]],\"b\")\n",
    "         ax3.set_title(\"Angle distribution of transmited photons\")\n",
    "         ax3.set_xlabel(\"theta (º)\")\n",
    "         ax3.set_ylabel(\"dn/dw (str^ (-1))\")\n",
    "         print(\"Plot Atrans  Angle distribution of transmitted photons\")  \n",
    "         plt.savefig('atrans')\n",
    "     \n",
    "     #  Angular distribution of photons \n",
    "     #  File: 'aretro.out'\n",
    "         with open(\"Aretro.out\",\"w\") as farefle:\n",
    "             xxar=np.linspace(0.5*dA+90,180-0.5*dA,B)\t\n",
    "             yyar=(numthr/(dA*gr*2*np.pi*np.sin(xxar*gr)))/nfot\n",
    "             for i in range(len(xxar)):\n",
    "                 xxar1 = xxar[i]\n",
    "                 yyar1 = yyar[0][i]\n",
    "                 farefle.write(\"%f\\t%f\\n\" % (xxar1,yyar1))\n",
    "             farefle.write(\" \")\n",
    "             print(\"File Aretro.out Angle distribution of backscattered photons\")   \n",
    "        \n",
    "     # Graph: 'aretro'                  \n",
    "         plt.figure()\n",
    "         ax4 = plt.axes()\n",
    "         ax4.set_xlim(90,180)\n",
    "         yyar=(numthr/(dA*gr*2*np.pi*np.sin(xxar*gr)))/nfot         \n",
    "         xx = np.repeat(xxar-0.5*dA, 2)[:-1]\n",
    "         yy = np.repeat(yyar, 2)[1:]          \n",
    "         plt.plot(xx,yy,\"b\")\n",
    "         plt.plot([xxar[B-1]-0.5*dA,180],[yyar[0][B-1],yyar[0][B-1]],\"b\")\n",
    "         ax4.set_title(\"Angle distribution of backscattered photons\")\n",
    "         ax4.set_xlabel(\"theta (º)\")\n",
    "         ax4.set_ylabel(\"dn/dw (str^ (-1))\")\n",
    "         print(\"Plot Aretro Angle distribution of backscattered photons\")\n",
    "         plt.savefig('aretro.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}