{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import maxwell\n",
    "\n",
    "me = 9.11e-28 # g\n",
    "R = 8.3144598e7 # erg/K/mol \n",
    "k = 1.38064852e-16 # erg/K\n",
    "c = 3.e10 # cm/s\n",
    "pi = 3.1415\n",
    "sigma=6.65e-25 # Thomson cross section in cgs [cm^2]\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def lorentz_transform(p,v):\n",
    "    #p is the four-element momentum vector p = h*nu/c * [1,kx, ky, kz], where kx^2+ky^2+kz^2=1\n",
    "    #The units are assumed to be in the cgs system\n",
    "    return(np.dot(lorentz_matrix(v),p))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def lorentz_matrix(v):\n",
    "    # v is assumed to be a 3-element vector with the x, y, and z components of the velocity in cgs units (cm/s)\n",
    "    beta=np.sqrt(np.sum(v**2))/c\n",
    "    beta_vec=v/c\n",
    "    gamma=1./np.sqrt(1. - beta**2)\n",
    "    matrix=np.zeros((4,4))\n",
    "    # general matrix form of Lorentz transform\n",
    "    # Instead of writing each element separately, one can write fewer lines of code with the following vector expressions\n",
    "    matrix[0,0]=gamma\n",
    "    matrix[1:,0]=-gamma*beta_vec\n",
    "    matrix[0,1:]=-gamma*beta_vec\n",
    "    #  The outer product of two coordinate vectors u and v is a matrix w, such that w[i,j] = u[i]*v[j]\n",
    "    matrix[1:,1:]=(gamma-1)*np.outer(beta_vec,beta_vec)\n",
    "    for i in range(1,4):\n",
    "        matrix[i,i]+=1\n",
    "    return(matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.07832773 -0.10783277 -0.21566555 -0.32349832]\n",
      " [-0.10783277  1.00078328  0.00156655  0.00234983]\n",
      " [-0.21566555  0.00156655  1.00313311  0.00469966]\n",
      " [-0.32349832  0.00234983  0.00469966  1.0070495 ]]\n"
     ]
    }
   ],
   "source": [
    "c = 3.e10 #cm/s \n",
    "print(lorentz_matrix(np.asarray([0.1*c,0.2*c,0.3*c])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def random_direction():\n",
    "    #Returns a random unit vector (k=[kx,ky,kz], where kx^2+ky^2+kz^2=1), drawn from an isotropic distribution\n",
    "    phi=2.*np.pi*np.random.rand()\n",
    "    # note: random orientation=even distribution in cos(theta)\n",
    "    cos_theta=2.*np.random.rand()-1\n",
    "    sin_theta=np.sqrt(1-cos_theta**2)\n",
    "    return(np.array([sin_theta*np.cos(phi),sin_theta*np.sin(phi),cos_theta]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def scatter(p_in,v):\n",
    "    v_vec=v*random_direction()\n",
    "    # first, transform to electron frame\n",
    "    p_0_prime=lorentz_transform(p_in,v_vec)\n",
    "    # neglect recoil, so energy unchanged\n",
    "    h_nu_prime=p_0_prime[0]\n",
    "    # assume isotropic scattering, draw random new direction\n",
    "    k_1_prime=random_direction()\n",
    "    # new four momentum of photon in electron frame\n",
    "    p_1_prime=h_nu_prime*np.array([1.,k_1_prime[0],k_1_prime[1],k_1_prime[2]])\n",
    "    # transform back to observer frame\n",
    "    p_out=lorentz_transform(p_1_prime,-v_vec)\n",
    "    \n",
    "    return(p_out)\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def move(n,sigma,H,x_vec,p_in,v):\n",
    "    # n: number density of electrons [cm^{-3}]\n",
    "    # sigma: the Thomson electron scattering cross section [cm^2]\n",
    "    # H: the height of the corona [cm]\n",
    "    # x_vec: a three element vector with the x,y,z coordinates of the photon [cm]\n",
    "    # p_in: the four momentum vector of the photon p = h*nu/c * [1,kx, ky, kz], where kx^2+ky^2+kz^2=1\n",
    "    # v: the magnitude of the velocity of the electron [cm/s]\n",
    "    \n",
    "    # draw random probability\n",
    "    p=np.random.rand()\n",
    "    # and find corresponding scattering depth\n",
    "    tau=-np.log(p)\n",
    "    # figure out how far and where the photon went\n",
    "    direction=p_in[1:]/np.sqrt(np.sum(p_in[1:]**2))\n",
    "    dist=tau/n/sigma\n",
    "    x_vec+=dist*direction\n",
    "    # determine if it escaped or got absorbed\n",
    "    if (x_vec[2]>H):\n",
    "        p_out=p_in\n",
    "    elif (x_vec[2]<0):\n",
    "        p_out=0*p_in\n",
    "    # if it neither, execute inverse Compton scatter\n",
    "    else:\n",
    "        p_out = scatter(p_in,v)\n",
    "    return(p_out)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# this is the version for mono-energetic electrons and photons\n",
    "def monte_carlo(n_phot,H,tau,v,h_nu_zero):\n",
    "    # run Monte Carlo\n",
    "    h_nus=np.zeros(n_phot) # final energies\n",
    "    out_z=np.zeros(n_phot) # final z-position\n",
    "    n=tau/H/sigma # electron density, given tau and H\n",
    "    for i in range(0,n_phot):\n",
    "        # draw random initial up-going photon direction\n",
    "        k_zero=random_direction()\n",
    "        k_zero[2]=np.abs(k_zero[2])\n",
    "        # initial photon four-momentum\n",
    "        p_in=h_nu_zero*np.array([1,k_zero[0],k_zero[1],k_zero[2]])\n",
    "        h_nus[i]=h_nu_zero\n",
    "        # initial position of photons\n",
    "        x_vec=np.zeros(3)\n",
    "        # keep scattering until absorbed or escaped\n",
    "        while (x_vec[2] < H and x_vec[2]>=0):\n",
    "            p_in=move(n,sigma,H,x_vec,p_in,v)\n",
    "            h_nus[i]=p_in[0]\n",
    "            out_z[i]=x_vec[2]\n",
    "    return(h_nus[out_z>0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# this is the version for mono-energetic electrons and photons\n",
    "def monte_carlo_maxwellian(n_phot,H,tau,Te,h_nu_zero):\n",
    "    # run Monte Carlo\n",
    "    h_nus=np.zeros(n_phot) # final energies\n",
    "    out_z=np.zeros(n_phot) # final z-position\n",
    "    sigma=6.65e-25 # Thomson cross section in cgs [cm^2]\n",
    "    n=tau/H/sigma # electron density, given tau and H\n",
    "    for i in range(0,n_phot):\n",
    "        # draw random initial up-going photon direction\n",
    "        k_zero=random_direction()\n",
    "        k_zero[2]=np.abs(k_zero[2])\n",
    "        # initial photon four-momentum\n",
    "        p_in=h_nu_zero*np.array([1,k_zero[0],k_zero[1],k_zero[2]])\n",
    "        h_nus[i]=h_nu_zero\n",
    "        # initial position of photons\n",
    "        x_vec=np.zeros(3)\n",
    "        # keep scattering until absorbed or escaped\n",
    "        while (x_vec[2] < H and x_vec[2]>=0):\n",
    "            #generate a velocity magnitude from a maxwellian distribution\n",
    "            sample = maxwell.rvs()\n",
    "            #We need to transform to physical units by doing the transformation sample**2 = me*v**2/(kT)\n",
    "            v = np.sqrt(sample**2 * k * Te / me)\n",
    "            p_in=move(n,sigma,H,x_vec,p_in,v)\n",
    "            h_nus[i]=p_in[0]\n",
    "            out_z[i]=x_vec[2]\n",
    "    return(h_nus[out_z>0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def plot_mc(tau,H,v,h_nu_zero):\n",
    "    h_nus=monte_carlo(n_phot,H,tau,v,h_nu_zero)\n",
    "    h_nus/=h_nu_zero\n",
    "    plt.hist(h_nus,bins=np.logspace(-1,2,num=100),log=True,label=r'$\\tau=${:4.1f}'.format(tau))\n",
    "    plt.xscale('log')\n",
    "    plt.xlim(0.5,100)\n",
    "    plt.xlabel(r'$h\\nu/h\\nu_{0}$')\n",
    "    plt.ylabel(r'$N(h\\nu)$')\n",
    "    plt.legend()\n",
    "    plt.show()\n",
    "    print('Fraction of transmitted photons: {0:5.3f}'.format(float(h_nus.size)/float(n_phot)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def plot_mc_maxwellian(tau,H,Te,h_nu_zero):\n",
    "    h_nus=monte_carlo_maxwellian(n_phot,H,tau,Te,h_nu_zero)\n",
    "    h_nus/=h_nu_zero\n",
    "    plt.hist(h_nus,bins=np.logspace(-1,2,num=100),log=True,label=r'$\\tau=${:4.1f}'.format(tau))\n",
    "    plt.xscale('log')\n",
    "    plt.xlim(0.5,100)\n",
    "    plt.xlabel(r'$h\\nu/h\\nu_{0}$')\n",
    "    plt.ylabel(r'$N(h\\nu)$')\n",
    "    plt.legend()\n",
    "    plt.show()\n",
    "    print('Fraction of transmitted photons: {0:5.3f}'.format(float(h_nus.size)/float(n_phot)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAERCAYAAAC6kZqPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFHtJREFUeJzt3X+MXeWd3/H3ZwHHDSWsAs62eLBsdlgWs4UQTUikKD+2\nZYuBdViRtGsvWjWLg8UqRKn6zzrSqlVURQI1laoIstSbIKdNYpdl04CDgbbaIKhEAiY0ibFLsQgp\nY7qLgYakq7CA+faPuV6G8Qxz75m599xjv1+SpTnPPefc70wO95PnOc89T6oKSZKa+KW2C5AkdZch\nIklqzBCRJDVmiEiSGjNEJEmNGSKSpMYMEUlSY4aIJKkxQ0SS1NjJbRewXM4888xau3Zt22VIUqc8\n+uijz1fVqqbHHzchsnbtWvbu3dt2GZLUKUl+spTjx3I4K8lHkjyY5NYkH2m7HknS/EYWIkluS/Jc\nkn1z2jckeSLJwSTbes0F/D9gJTA9qholSYMZZU9kB7BhdkOSk4BbgMuB9cDmJOuBB6vqcuCPgM+N\nsEZJ0gBGdk+kqh5IsnZO8yXAwap6CiDJLuCqqtrfe/3/Am9b6JxJtgJbAdasWbPcJUvqgFdffZXp\n6WlefvnltksZaytXrmRiYoJTTjllWc/b9o311cAzs7angfcluRq4DPhl4OaFDq6q7cB2gKmpKRdG\nkU5A09PTnHbaaaxdu5YkbZczlqqKF154genpadatW7es5247ROZVVd8EvtnPvkk2AhsnJyeHW5Sk\nsfTyyy8bIItIwhlnnMHhw4eX/dxtz846BJw9a3ui19a3qtpdVVtPP/30ZS1MUncYIIsb1t+o7RB5\nBDg3ybokK4BNwF2DnCDJxiTbX3rppaEUKEla2MiGs5LsBD4CnJlkGvhXVfWVJDcA9wEnAbdV1eOD\nnLeqdgO7p6amrlvumodp7ba737T99I1XtlSJJDU3ytlZmxdo3wPsaXpe74lIUnvG8sb6ILraE5E0\nHHN7+UvV5ijBvffey2c+8xmOHDnCJz/5SbZt2zbvftdeey3f/va3ede73sW+ffvm3WdY2r4nIkma\nx5EjR/jUpz7FPffcw/79+9m5cyf79++fd99PfOIT3HvvvSOucEbnQ8Qb65La9rOf/YyLL76YCy64\ngLe//e28+93v5v3vfz+vv/5643M+/PDDTE5Ocs4557BixQo2bdrEnXfeOe++H/rQh3jnO9/Z+L2W\nwuEsSVqid7zjHTz22GM8/PDDfP7znz/mw/6DH/wgP//5z4857gtf+AKXXnrpvOc8dOgQZ5/9xjcg\nJiYm+N73vre8hS+DzoeIJI2Lffv2ccEFFxzT/uCDD7ZQzWh0PkScnSVpXOzfv5/3vOc9x7Q36Yms\nXr2aZ55546lQ09PTrF69evmKXSadDxGHsySNi2effZYrrrjimPYmPZH3vve9PPnkk/z4xz9m9erV\n7Nq1i2984xvLUeay6nyISNJsbU7Jveyyy9iyZQs7duzgwx/+8JLOdfLJJ3PzzTdz2WWXceTIEa69\n9tq/HSq74oor+PKXv8xZZ50FwObNm7n//vt5/vnnmZiY4HOf+xxbtmxZ8u/Tj1QdHw+/nZqaqi4t\nj+s31qXlceDAAc4///y2y+iE+f5WSR6tqqmm53SKrySpsc6HiE/xlaT2dD5EJOl4GZYfpmH9jQwR\nSZ22cuVKXnjhBYPkLRxd2XDlypXLfm5nZ0nqtImJCaanp4eyat/x5Oga68ut8yHilw2lE9spp5yy\n7OuGq3+dH87yxroktafzISJJao8hIklqzBCRJDVmiEiSGjNEJEmNdT5EfHaWJLWn8yHiFF9Jak/n\nQ0SS1B5DRJLUmCEiSWrMEJEkNWaISJIaM0QkSY0ZIpKkxsY2RJKcmmRvkt9uuxZJ0vxGFiJJbkvy\nXJJ9c9o3JHkiycEk22a99EfA7aOqT5I0uFH2RHYAG2Y3JDkJuAW4HFgPbE6yPslvAfuB50ZYnyRp\nQCNbHreqHkiydk7zJcDBqnoKIMku4Crg7wKnMhMsv0iyp6pen3vOJFuBrQBr1qwZXvGSpHm1vcb6\nauCZWdvTwPuq6gaAJJ8Anp8vQACqajuwHWBqaqqGW6okaa62Q+QtVdWOxfZJshHYODk5OfyCJElv\n0vbsrEPA2bO2J3ptffMpvpLUnrZD5BHg3CTrkqwANgF3DXIC1xORpPaMcorvTuAh4Lwk00m2VNVr\nwA3AfcAB4PaqenyQ89oTkaT2jHJ21uYF2vcAe0ZVhyRp+bQ9nLVkDmdJUns6HyIOZ0lSezofIvZE\nJKk9nQ8ReyKS1J7Oh4gkqT2GiCSpsc6HiPdEJKk9nQ8R74lIUns6HyKSpPYYIpKkxsb6UfD96MKj\n4Nduu7vtEiRpKDrfE/GeiCS1p/MhIklqjyEiSWrMEJEkNdb5EPHLhpLUns6HiDfWJak9nQ8RSVJ7\nDBFJUmOd/7Lh8WK+LyQ+feOVLVQiSf2zJyJJaswQkSQ11vkQcYqvJLWn8/dEqmo3sHtqauq6tms5\nygcuSjpRdL4nIklqjyEiSWrMEJEkNWaISJIaM0QkSY0ZIpKkxgwRSVJjYxkiSc5PcmuSO5L8Ydv1\nSJLmN7IQSXJbkueS7JvTviHJE0kOJtkGUFUHqup64J8CHxhVjZKkwYyyJ7ID2DC7IclJwC3A5cB6\nYHOS9b3XPgrcDewZYY2SpAGMLESq6gHgxTnNlwAHq+qpqnoF2AVc1dv/rqq6HLhmoXMm2Zpkb5K9\nhw8fHlbpkqQFtP3srNXAM7O2p4H3JfkIcDXwNt6iJ1JV24HtAFNTUzW8MiVJ82k7ROZVVfcD9/ez\nb5KNwMbJyclhliRJmkfbs7MOAWfP2p7otfWtqnZX1dbTTz99WQuTJC2u7RB5BDg3ybokK4BNwF2D\nnMD1RCSpPaOc4rsTeAg4L8l0ki1V9RpwA3AfcAC4vaoeH+S89kQkqT0juydSVZsXaN/DEqbxek9E\nktrT9nDWktkTkaT2dD5EJEntGcspvoM4noez5q7V/vSNV7ZUiSTNr/M9EYezJKk9A4dIklN7z7yS\nJJ3gFg2RJL+U5PeS3J3kOeB/Av8nyf4k/yZJq+NIfk9EktrTzz2R7wD/DfgssK+qXgdI8k7gN4Gb\nkvznqvra8MpcWFXtBnZPTU1d18b7w7H3LiTpRNFPiFxaVa/ObayqF4E/B/48ySnLXpkkaewtOpx1\nNECSnJbkA0l+daF9JEknlkGm+H4H+AHw60kOA9f2eiOtOp6n+ErSuOt7dlZVTVXVlqr6APA1YE+S\ns4ZXWt91OcVXkloy0BTf3kyt1cBPgHuZ6Z1Ikk5QfQ9nJflL4O8AfwU82/u3e0h1SZI6YJB7Ir9W\nVT8bWiWSpM4ZJESmklwD/BTYB/yQme+N/M1QKuuTN9YlqT2D3BO5jZnhq+8C5wD/EhhoAalh8Ma6\nJLVnkJ7IT6rqW72f/2wYxUiSuqWfZ2f9hyT/HPhukn8xgpokSR3Rz3DWDiDArwC/n+QnSe5K8q+T\n/JOhVidJGmuLDmdV1V8Af3F0O8nJwPnARcAlOLQlSSesRUMkSaqqjm5X1WvAj3r/vjbfPpKkE0M/\nw1nfSfLpJGtmNyZZkeQfJvkq8M+GU97iXE9EktrTT4hsAI4AO5McXYzqx8CTwGbg31XVjiHW+Jac\n4itJ7elniu+vVNWXgC/11g05E/hFVf10uKVJksZdPyFyT5J3MbMs7g+ZuRfywyQ/8jEoknRi62d2\n1vokbwPWA/8AuBC4Crgwyd9U1boh1yhJGlN9fWO993ysx5IcBH7BzJDWrzHTM5EknaD6meJ7HnAl\n8NvAKuC/Al8HtlbVK8MtT5I0zvrpiRwAHgNuAu5s+6m9kqTx0U+I/CHwG8CngJuTvMAbXzb80ayH\nMkqSTjD93Fj/97O3k0zwxg32jwFDCZEkv8PMMNo7gK9U1X8ZxvtIkpobaI11gKqarqp7quqmqvr9\nQY5NcluS55Lsm9O+IckTSQ4m2dZ7n29V1XXA9cDvDlqnJGn4Bg6RJdrBzDfg/1aSk4BbgMuZmUa8\nOcn6Wbv8ce91SdKYGWmIVNUDwItzmi8BDlbVU73ZXruAqzLjJuCeqvr+fOdLsjXJ3iR7Dx8+PNzi\nJUnHGHVPZD6rgWdmbU/32j4NXAp8PMn18x1YVduraqqqplatWjX8SiVJbzLI8rgjVVVfBL642H5J\nNgIbJycnh1+UJOlNxiFEDgFnz9qe6LX1pap2A7unpqauW+7Cxs3abXe/afvpG69sqRJJmjEOw1mP\nAOcmWZdkBbAJuKvfg11PRJLaM9IQSbITeAg4L8l0ki29lRJvAO5j5tvxt1fV4/2e0/VEJKk9Ix3O\nqqrNC7TvAfaMshZJ0tKNw3DWkjicJUnt6XyIOJwlSe0Zh9lZSzLqKb5zZ0hJ0onMnogkqbHOh4gk\nqT2GiCSpsc6HiLOzJKk9nQ8R74lIUns6HyKSpPZ0PkQczpKk9nQ+RBzOkqT2dD5EJEntMUQkSY0Z\nIpKkxgwRSVJjnQ8RZ2dJUns6HyLOzpKk9nQ+RCRJ7en8eiInsvnWNnn6xitbqETSicqeiCSpMUNE\nktSYISJJaqzzIeIUX0lqT+dDxCm+ktQeZ2ctYr4ZUJKkGZ3viUiS2mOISJIaM0QkSY0ZIpKkxgwR\nSVJjhogkqbGxDJEk5yT5SpI72q5FkrSwkYVIktuSPJdk35z2DUmeSHIwyTaAqnqqqraMqjZJUjOj\n7InsADbMbkhyEnALcDmwHticZP0Ia5IkLcHIQqSqHgBenNN8CXCw1/N4BdgFXNXvOZNsTbI3yd7D\nhw8vY7WSpH60fU9kNfDMrO1pYHWSM5LcClyc5LMLHVxV26tqqqqmVq1aNexaJUlzjOWzs6rqBeD6\nfvZNshHYODk5OdyiJEnHaLsncgg4e9b2RK+tbz7FV5La03aIPAKcm2RdkhXAJuCuQU7geiKS1J5R\nTvHdCTwEnJdkOsmWqnoNuAG4DzgA3F5Vjw9yXnsiktSekd0TqarNC7TvAfY0Pa/3RCSpPW0PZy2Z\nPRFJak/nQ0SS1J7Oh4g31iWpPZ0PEYezJKk9nQ8RSVJ7xvIb64Nwdtabrd1295u2n77xypYqkXQi\n6HxPxOEsSWpP50NEktQeQ0SS1FjnQ8QpvpLUns6HiPdEJKk9nQ8RSVJ7DBFJUmOGiCSpMb9sOMvc\nL+pJkt5a53si3liXpPZ0PkQkSe0xRCRJjRkikqTGDBFJUmOGiCSpsc6HiM/OkqT2dD5EnOIrSe3p\nfIhIktpjiEiSGjNEJEmNGSKSpMYMEUlSY4aIJKkxQ0SS1NhYrieS5FTgS8ArwP1V9fWWS5IkzWNk\nPZEktyV5Lsm+Oe0bkjyR5GCSbb3mq4E7quo64KOjqlGSNJhRDmftADbMbkhyEnALcDmwHticZD0w\nATzT2+3ICGuUJA1gZMNZVfVAkrVzmi8BDlbVUwBJdgFXAdPMBMn/4C2CLslWYCvAmjVrBq7pRFgO\nd77f8ekbr2yhEknHo7ZvrK/mjR4HzITHauCbwMeS/Amwe6GDq2p7VU1V1dSqVauGW6kk6RhjeWO9\nqv4a+IN+9k2yEdg4OTk53KIkScdouydyCDh71vZEr61vPsVXktrTdog8ApybZF2SFcAm4K5BTuB6\nIpLUnlFO8d0JPAScl2Q6yZaqeg24AbgPOADcXlWPD3JeeyKS1J5Rzs7avED7HmDPqOqQJC2ftoez\nlszhLElqT+dDxOEsSWpP50PEnogktSdV1XYNyyLJS8CTAxxyOrBY8ix1n7d67Uzg+UXOPU76+VuM\ny3ss5TyDHNvvvovt1/R1r6Hhvk/T8wx6XNufRedV1WmLnHthVXVc/AO2L/f+S91nkdf2tv03G+bf\nt833WMp5Bjm2330X26/p615D43kdnWifRZ0fzpplwcejLGH/pe4zaE3jbBS/y3K9x1LOM8ix/e67\n2H5Lfb0rRvV7tH0dnVCfRcfNcFbXJNlbVVNt16Hu8hrScljqdXQ89US6ZnvbBajzvIa0HJZ0HdkT\nkSQ1Zk9EktSYISJJaswQkSQ1ZohIkhozRMZAklOTfDXJnya5pu161E1JzknylSR3tF2LuivJ7/Q+\ni/5Tkn+82P6GyJAkuS3Jc0n2zWnfkOSJJAeTbOs1Xw3cUVXXAR8debEaW4NcR1X1VFVtaadSjbMB\nr6Nv9T6Lrgd+d7FzGyLDswPYMLshyUnALcDlwHpgc5L1zCwL/ExvtyMjrFHjbwf9X0fSQnYw+HX0\nx73X35IhMiRV9QDw4pzmS4CDvf/H+AqwC7gKmGYmSMD/TTTLgNeRNK9BrqPMuAm4p6q+v9i5/cAa\nrdW80eOAmfBYDXwT+FiSP+H4eU6Shmfe6yjJGUluBS5O8tl2SlOHLPR59GngUuDjSa5f7CQjWx5X\nC6uqvwb+oO061G1V9QIz49hSY1X1ReCL/e5vT2S0DgFnz9qe6LVJg/A60nJYluvIEBmtR4Bzk6xL\nsgLYBNzVck3qHq8jLYdluY4MkSFJshN4CDgvyXSSLVX1GnADcB9wALi9qh5vs06NN68jLYdhXkc+\nxVeS1Jg9EUlSY4aIJKkxQ0SS1JghIklqzBCRJDVmiEiSGjNEJEmNGSKSpMYMEalPSf5Rkv84wP63\nJvlAk2OlrjBEpP5dBDw2wP7vB77b8FipEwwRqX8XAX8vyQNJ/neSS5OcnuSvju6Q5NFe2/nA/6qq\nI4MeO9pfSVoaQ0Tq30XA4ar6EPAZ4Jqqegl4e5Kja/P8ALiQmSVH7214rNQZhojUhySnAGcA/7bX\ndArw097Pfwn8/d7Pv97bvoxeiAx6bJJTk3w1yZ8muWY4v5G0PAwRqT/nAz+oqtd72xcC+3o/Pwuc\nleTjwPPMLOzzy1X17KDHVtWTwNXAHVV1HfDRYf5S0lK5PK7Un4uYGW466kLgzt7PzwJXMDOEdQXw\nm8B3Gh4LMyvM/aj38xGkMWZPROrPRcAPZ23/Bm/uTfwe8LGqep7574f0eyzANDNBAv43qjHnolTS\nMkvyfeB9VfVqw+NPBW4GXgb+e1V9fTnrk5aTISJJasyusiSpMUNEktSYISJJaswQkSQ1ZohIkhoz\nRCRJjRkikqTGDBFJUmP/H/62G3TfYf+HAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1a17a7f810>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fraction of transmitted photons: 0.859\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAERCAYAAAC6kZqPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFHFJREFUeJzt3X+MXeWd3/H3ZwHHDSWsAs62zGDZ1CyL2UKIJiRSlB/b\nssXAOqxIdtcOWjWLg8UqRKn6zzrSqlVURQI1laoIstSbIKdNYpeyaYDEQFttEKxEAiZsEmOXYhFS\nxnQXA42TrsIC5ts/5rIMwwxz75m599zjeb8kS3Oee84z35kc5pPnPM85J1WFJElN/FLbBUiSussQ\nkSQ1ZohIkhozRCRJjRkikqTGDBFJUmOGiCSpMUNEktSYISJJauzEtgtYLqeffnqtW7eu7TIkqVMe\nfvjhZ6tqTdPjj5sQWbduHfv27Wu7DEnqlCQ/WcrxY3k5K8mHktyf5OYkH2q7HknS/EYWIkluSfJM\nkv1z2jcleSzJoSQ7es0F/D9gNTA9qholSYMZ5UhkF7BpdkOSE4CbgEuBjcDWJBuB+6vqUuCPgM+O\nsEZJ0gBGNidSVfclWTen+SLgUFU9AZBkD3BFVR3off5/gbcs1GeS7cB2gLVr1y53yZI64KWXXmJ6\nepoXXnih7VLG2urVq5mcnOSkk05a1n7bnlifAJ6atT0NvCfJlcAlwC8DNy50cFXtBHYCTE1N+WIU\naQWanp7mlFNOYd26dSRpu5yxVFU899xzTE9Ps379+mXtu+0QmVdVfQP4Rj/7JtkMbN6wYcNwi5I0\nll544QUDZBFJOO200zhy5Miy99326qzDwJmztid7bX2rqjuravupp566rIVJ6g4DZHHD+h21HSIP\nAWcnWZ9kFbAFuGOQDpJsTrLz6NGjQylQkrSwkV3OSrIb+BBwepJp4F9X1ZeTXAfcA5wA3FJVjw7S\nb1XdCdw5NTV1zXLXPEzrdnz7ddtPXn95S5VIUnOjXJ21dYH2vcDepv06JyJJ7RnLifVBdHUkImk4\n5o7yl6rNqwR33303n/70pzl27Bif+MQn2LFjx7z7XX311XzrW9/iHe94B/v37593n2Fpe05EkjSP\nY8eO8clPfpK77rqLAwcOsHv3bg4cODDvvh//+Me5++67R1zhjM6HiBPrktr2s5/9jAsvvJDzzjuP\nt771rbzzne/kve99L6+88krjPh988EE2bNjAWWedxapVq9iyZQu33377vPt+4AMf4O1vf3vj77UU\nXs6SpCV629vexiOPPMKDDz7I5z73uTf8sX//+9/Pz3/+8zcc9/nPf56LL7543j4PHz7MmWe+dgfE\n5OQk3/ve95a38GXQ+RCRpHGxf/9+zjvvvDe033///S1UMxqdDxFXZ0kaFwcOHOBd73rXG9qbjEQm\nJiZ46qnXngo1PT3NxMTE8hW7TDofIl7OkjQunn76aS677LI3tDcZibz73e/m8ccf58c//jETExPs\n2bOHr3/968tR5rLqfIhI0mxtLsm95JJL2LZtG7t27eKDH/zgkvo68cQTufHGG7nkkks4duwYV199\n9d9dKrvsssv40pe+xBlnnAHA1q1buffee3n22WeZnJzks5/9LNu2bVvyz9OPVB0fD7+dmpqqLr0e\n1zvWpeVx8OBBzj333LbL6IT5fldJHq6qqaZ9usRXktRY50PEp/hKUns6HyKSdLxclh+mYf2ODBFJ\nnbZ69Wqee+45g+RNvPpmw9WrVy97367OktRpk5OTTE9PD+WtfceTV9+xvtw6HyLebCitbCeddNKy\nvzdc/ev85Swn1iWpPZ0PEUlSewwRSVJjhogkqTFDRJLUmCEiSWqs8yHis7MkqT2dDxGX+EpSezof\nIpKk9hgikqTGDBFJUmOGiCSpMUNEktSYISJJaswQkSQ1NrYhkuTkJPuS/FbbtUiS5jeyEElyS5Jn\nkuyf074pyWNJDiXZMeujPwJuHVV9kqTBjXIksgvYNLshyQnATcClwEZga5KNSX4TOAA8M8L6JEkD\nGtnrcavqviTr5jRfBByqqicAkuwBrgD+PnAyM8HyiyR7q+qVuX0m2Q5sB1i7du3wipckzavtd6xP\nAE/N2p4G3lNV1wEk+Tjw7HwBAlBVO4GdAFNTUzXcUiVJc7UdIm+qqnYttk+SzcDmDRs2DL8gSdLr\ntL066zBw5qztyV5b33yKryS1p+0QeQg4O8n6JKuALcAdg3Tg+0QkqT2jXOK7G3gAOCfJdJJtVfUy\ncB1wD3AQuLWqHh2kX0ciktSeUa7O2rpA+15g76jqkCQtn7YvZy2Zl7MkqT2dDxEvZ0lSezofIo5E\nJKk9nQ8RRyKS1J7Oh4gkqT2GiCSpsc6HiHMiktSezoeIcyKS1J7Oh4gkqT1j/RTf48W6Hd9uuwRJ\nGorOj0ScE5Gk9nQ+RJwTkaT2dD5EJEntMUQkSY0ZIpKkxjofIk6sS1J7Oh8iTqxLUns6HyKSpPYY\nIpKkxgwRSVJjhogkqTFDRJLUWOdDxCW+ktSezoeIS3wlqT0+Cn5MzPe4+Cevv7yFSiSpf50fiUiS\n2mOISJIa83LWEPgmQ0krhSMRSVJjhogkqTFDRJLU2FiGSJJzk9yc5LYkf9h2PZKk+Y0sRJLckuSZ\nJPvntG9K8liSQ0l2AFTVwaq6Fvhd4H2jqlGSNJhRjkR2AZtmNyQ5AbgJuBTYCGxNsrH32YeBbwN7\nR1ijJGkAIwuRqroPeH5O80XAoap6oqpeBPYAV/T2v6OqLgWuWqjPJNuT7Euy78iRI8MqXZK0gLbv\nE5kAnpq1PQ28J8mHgCuBt/AmI5Gq2gnsBJiamqrhlSlJmk/bITKvqroXuLeffZNsBjZv2LBhmCVJ\nkubR9uqsw8CZs7Yne2198ym+ktSetkPkIeDsJOuTrAK2AHcM0oHvE5Gk9oxyie9u4AHgnCTTSbZV\n1cvAdcA9wEHg1qp6dJB+HYlIUntGNidSVVsXaN/LEpbxOiciSe0Zy4n1QVTVncCdU1NT17Rdy3Kb\n+zRgX1Ilady0PSciSeqwzoeIE+uS1J7Oh4gT65LUnoFDJMnJvWdeSZJWuEVDJMkvJflYkm8neQb4\nn8D/SXIgyb9N0uqyKC9nSVJ7+lmd9R3gfwCfAfZX1SsASd4O/AZwQ5L/WlVfHV6ZCxuH1Vm+U13S\nStVPiFxcVS/Nbayq54E/A/4syUnLXpkkaewtejnr1QBJckqS9yX5RwvtI0laWQa52fA7wA+AX0ty\nBLi6NxpplXesS1J7+l6dVVVTVbWtqt4HfBXYm+SM4ZXWd10u8ZWklgy0xLe3UmsC+AlwNzOjE0nS\nCtX35awkfwX8PeCvgad7/+4cUl2SpA4YZE7kV6vqZ0OrRJLUOYOEyFSSq4CfAvuBHzJz38jfDqWy\nPjmxLkntGWRO5BZmLl99FzgL+FfAQC+QGgYn1iWpPYOMRH5SVd/sff1fhlGMJKlb+nl21n9M8i+A\n7yb5lyOoSZLUEf1cztoFBPgV4PeT/CTJHUn+TZLfGWp1kqSxtujlrKr6c+DPX91OciJwLnABcBFe\n2pKkFWvREEmSqqpXt6vqZeBHvX9fnW8fSdLK0M/lrO8k+VSStbMbk6xK8k+SfAX458Mpb3G+T0SS\n2tNPiGwCjgG7k7z6MqofA48DW4F/X1W7hljjm3KJryS1p58lvr9SVV8Evth7b8jpwC+q6qfDLU1z\nzX351ZPXX95SJZI0o58QuSvJO5h5Le4PmZkL+WGSH/kYFEla2fpZnbUxyVuAjcA/Bs4HrgDOT/K3\nVbV+yDVKksZUX3es956P9UiSQ8AvmLmk9avMjEwkSStUP0t8zwEuB34LWAP8d+BrwPaqenG45UmS\nxlk/I5GDwCPADcDtbT+1t21zJ7claSXrJ0T+EPh14JPAjUme47WbDX8066GMkqQVpp+J9f8wezvJ\nJK9NsH8EGEqIJPltZi6jvQ34clX9t2F8H0lScwO9Yx2gqqar6q6quqGqfn+QY5PckuSZJPvntG9K\n8liSQ0l29L7PN6vqGuBa4PcGrVOSNHwDh8gS7WLmDvi/k+QE4CbgUmaWEW9NsnHWLn/c+1ySNGZG\nGiJVdR/w/Jzmi4BDVfVEb7XXHuCKzLgBuKuqvj9ff0m2J9mXZN+RI0eGW7wk6Q1GPRKZzwTw1Kzt\n6V7bp4CLgY8muXa+A6tqZ1VNVdXUmjVrhl+pJOl1Bnk97khV1ReALyy2X5LNwOYNGzYMvyhJ0uuM\nw0jkMHDmrO3JXltffIqvJLVnHELkIeDsJOuTrAK2AHf0e7DvE5Gk9ow0RJLsBh4AzkkynWRb702J\n1wH3MHN3/K1V9Wi/fToSkaT2jHROpKq2LtC+F9g7ylokSUs3thPr/VrJE+vzPcfLF1VJGqVxmBNZ\nEi9nSVJ7Oh8iTqxLUns6HyKORCSpPZ0PEUlSewwRSVJjnQ8R50QkqT2dDxHnRCSpPZ0PEUlSezof\nIl7OkqT2dP6O9aq6E7hzamrqmmH0P99d4ZKkGZ0fiUiS2mOISJIaM0QkSY0ZIpKkxjofIq7OkqT2\ndD5EvNlQktrT+RCRJLXHEJEkNWaISJIaM0QkSY0ZIpKkxjofIi7xlaT2dD5EXOIrSe3pfIhIktrT\n+UfB6/XmPrr+yesvb6kSSSuBIxFJUmOGiCSpMUNEktSYISJJaswQkSQ1NpYhkuSsJF9OclvbtUiS\nFjayEElyS5Jnkuyf074pyWNJDiXZAVBVT1TVtlHVJklqZpQjkV3AptkNSU4AbgIuBTYCW5NsHGFN\nkqQlGFmIVNV9wPNzmi8CDvVGHi8Ce4Ar+u0zyfYk+5LsO3LkyDJWK0nqR9tzIhPAU7O2p4GJJKcl\nuRm4MMlnFjq4qnZW1VRVTa1Zs2bYtUqS5hjLx55U1XPAtf3sm2QzsHnDhg3DLUqS9AZtj0QOA2fO\n2p7stfXNp/hKUnvaHok8BJydZD0z4bEF+NggHSznSGTuwwslSW9ulEt8dwMPAOckmU6yrapeBq4D\n7gEOArdW1aOD9OtIRJLaM7KRSFVtXaB9L7C3ab/OiUhSe9qeE1kyRyKS1J7Oh4gkqT2dD5Ekm5Ps\nPHr0aNulSNKK0/kQ8XKWJLWn8yEiSWpP2/eJLJmrs97cfPe+PHn95S1UIul41PmRiJezJKk9nQ8R\nSVJ7DBFJUmOdDxGX+EpSezofIs6JSFJ7Oh8ikqT2GCKSpMYMEUlSYyv6ZkNfQiVJS9P5kYgT65LU\nns6HiCSpPYaIJKkxQ0SS1JghIklqzBCRJDW2opf4rlT9LG32nSOS+tH5kYhLfCWpPZ0PEUlSewwR\nSVJjhogkqTFDRJLUmCEiSWrMEJEkNWaISJIaG8ubDZOcDHwReBG4t6q+1nJJkqR5jGwkkuSWJM8k\n2T+nfVOSx5IcSrKj13wlcFtVXQN8eFQ1SpIGM8rLWbuATbMbkpwA3ARcCmwEtibZCEwCT/V2OzbC\nGiVJAxjZ5ayqui/JujnNFwGHquoJgCR7gCuAaWaC5C95k6BLsh3YDrB27drXfearb5fX3N+nz9aS\nBO1PrE/w2ogDZsJjAvgG8JEkfwLcudDBVbWzqqaqamrNmjXDrVSS9AZjObFeVX8D/EE/+/oUX0lq\nT9sjkcPAmbO2J3ttffMpvpLUnrZD5CHg7CTrk6wCtgB3DNJBks1Jdh49enQoBUqSFjbKJb67gQeA\nc5JMJ9lWVS8D1wH3AAeBW6vq0UH6dSQiSe0Z5eqsrQu07wX2jqoOSdLyafty1pJ5OUuS2tP5EPFy\nliS1p/Mh4khEktqTqmq7hmWR5Cjw+ACHnAosljxL3efNPjsdeHaRvsdJP7+LcfkeS+lnkGP73Xex\n/Zp+7jk03O/TtJ9Bj2v7b9E5VXXKIn0vrKqOi3/AzuXef6n7LPLZvrZ/Z8P8/bb5PZbSzyDH9rvv\nYvs1/dxzaDzPo5X2t6jzl7NmWfDxKEvYf6n7DFrTOBvFz7Jc32Mp/QxybL/7LrbfUj/vilH9HG2f\nRyvqb9Fxczmra5Lsq6qptutQd3kOaTks9Tw6nkYiXbOz7QLUeZ5DWg5LOo8ciUiSGnMkIklqzBCR\nJDVmiEiSGjNEJEmNGSJjIMnJSb6S5E+TXNV2PeqmJGcl+XKS29quRd2V5Ld7f4v+c5J/ttj+hsiQ\nJLklyTNJ9s9p35TksSSHkuzoNV8J3FZV1wAfHnmxGluDnEdV9URVbWunUo2zAc+jb/b+Fl0L/N5i\nfRsiw7ML2DS7IckJwE3ApcBGYGuSjcy8Fvip3m7HRlijxt8u+j+PpIXsYvDz6I97n78pQ2RIquo+\n4Pk5zRcBh3r/j/FFYA9wBTDNTJCA/5tolgHPI2leg5xHmXEDcFdVfX+xvv2DNVoTvDbigJnwmAC+\nAXwkyZ9w/DwnScMz73mU5LQkNwMXJvlMO6WpQxb6e/Qp4GLgo0muXayTkb0eVwurqr8B/qDtOtRt\nVfUcM9expcaq6gvAF/rd35HIaB0Gzpy1PdlrkwbheaTlsCznkSEyWg8BZydZn2QVsAW4o+Wa1D2e\nR1oOy3IeGSJDkmQ38ABwTpLpJNuq6mXgOuAe4CBwa1U92madGm+eR1oOwzyPfIqvJKkxRyKSpMYM\nEUlSY4aIJKkxQ0SS1JghIklqzBCRJDVmiEiSGjNEJEmNGSJSn5L80yT/aYD9b07yvibHSl1hiEj9\nuwB4ZID93wt8t+GxUicYIlL/LgD+QZL7kvzvJBcnOTXJX7+6Q5KHe23nAv+rqo4NeuxofyRpaQwR\nqX8XAEeq6gPAp4Grquoo8NYkr76b5wfA+cy8cvTuhsdKnWGISH1IchJwGvDvek0nAT/tff1XwD/s\nff1rve1L6IXIoMcmOTnJV5L8aZKrhvMTScvDEJH6cy7wg6p6pbd9PrC/9/XTwBlJPgo8y8yLfX65\nqp4e9Niqehy4Eritqq4BPjzMH0paKl+PK/XnAmYuN73qfOD23tdPA5cxcwnrMuA3gO80PBZm3jD3\no97Xx5DGmCMRqT8XAD+ctf3rvH408THgI1X1LPPPh/R7LMA0M0EC/jeqMedLqaRlluT7wHuq6qWG\nx58M3Ai8APxFVX1tOeuTlpMhIklqzKGyJKkxQ0SS1JghIklqzBCRJDVmiEiSGjNEJEmNGSKSpMYM\nEUlSY/8fzesjR7uQwkAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1a17ac0050>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fraction of transmitted photons: 0.858\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAERCAYAAAC6kZqPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEohJREFUeJzt3X+sX3V9x/HnS2jtZIixwJR7a1pWhhQHYq5o4tT9wFHR\nilG3wdgyR4VhxLjsn2HiP2Yx0cwl08DmKrLq1DLGFmmVHy6ZBJegUCRqaefocI4LU0qZaBYZUN/7\n436r18tt7/d77vf7Pd/TPh9Jk3s+3/Pjfa+H78vP+ZxzPqkqJElq4lltFyBJ6i5DRJLUmCEiSWrM\nEJEkNWaISJIaM0QkSY0ZIpKkxgwRSVJjhogkqbFj2y5gWE488cRau3Zt22VIUqfcc889j1bVSU23\nP2JCZO3atezcubPtMiSpU5J8ZznbezlLktSYISJJaqzzIZJkU5Itjz/+eNulSNJRp/NjIlW1A9gx\nMzNzWdu1SBq/p556itnZWZ544om2S5loq1atYnp6mhUrVgx1v50PEUlHt9nZWY4//njWrl1LkrbL\nmUhVxf79+5mdnWXdunVD3XfnL2dJOro98cQTrF692gA5jCSsXr16JL01Q0RS5xkgSxvV38gQkSQ1\n5phIS9Ze9YWfWf7PD76hpUokqTl7IpKkxuyJSDqiLOzlL1ebVwkuvfRSPv/5z3PyySeza9euQ653\n66238p73vIcDBw7wjne8g6uuumpsNXa+J+LDhpKOVG9/+9u59dZbD7vOgQMHeNe73sUtt9zC7t27\n2bZtG7t37x5ThUdAiFTVjqq6/IQTTmi7FElHqR/84Aecc845nHnmmTznOc/hpS99Ka985Sv58Y9/\nvKz9vuY1r+H5z3/+Yde56667WL9+PaeeeiorV67koosu4qabblrWcQfh5SxJWqbnPve53Hvvvdx1\n11184AMfeMaX+Ktf/Wp++MMfPmO7D3/4w5x33nnLOvZDDz3EmjVrfrI8PT3NV7/61WXtcxCGiCQN\nya5duzjzzDOf0f7lL3+5hWrGwxCRpCHZvXs3L3vZy57RPsqeyNTUFA8++OBPlmdnZ5mamlrWPgdh\niEjSkDz88MNccMEFz2gfZU/k5S9/Offffz/f/va3mZqa4vrrr+ezn/3syI63kCEi6YjS5i25559/\nPps3b2br1q289rWvXfb+Lr74Ym6//XYeffRRpqenef/738/mzZu54IILuPbaaznllFM49thjufrq\nqzn//PM5cOAAl1566aKX1EYlVTW2g43SzMxMdWl6XJ9Yl4Zjz549nHHGGW2X0QmL/a2S3FNVM033\n2flbfCVJ7TFEJEmNGSKSOu9IuSw/SqP6G01siCQ5LsnOJG9suxZJk2vVqlXs37/fIDmMgzMbrlq1\nauj7HtvdWUmuA94IPFJVL5nXvhH4CHAMcG1VfbD30Z8CN4yrPkndND09zezsLPv27Wu7lIl2cI71\nYRvnLb5bgauBTx1sSHIMcA3wOmAWuDvJdmAK2A0MPzYlHVFWrFgx9HnD1b+xhUhV3ZFk7YLmc4G9\nVfUAQJLrgQuBnweOAzYAP0pyc1U9401mSS4HLgd40YteNLriJUmLavthwyngwXnLs8ArqupKgCRv\nBx5dLEAAqmoLsAXmnhMZbamSpIXaDpHDqqqtbdcgSTq0tu/OeghYM295utfWNyelkqT2tB0idwOn\nJVmXZCVwEbB9kB04KZUktWdsIZJkG3AncHqS2SSbq+pp4ErgNmAPcENV3Tfgfu2JSFJLfAHjGCx8\n2eJifAGjpDb4AkZJUms6HyJezpKk9nQ+RBxYl6T2dD5EJEnt6XyIeDlLktrT+RDxcpYktafzISJJ\nas9EvzvraLLYsyQ+OyJp0nW+J+KYiCS1p/Mh4piIJLWn8yEiSWqPYyIj0M+7siTpSGBPRJLUWOdD\nxIF1SWpP50PEgXVJak/nQ0SS1B5DRJLUmCEiSWrMEJEkNdb5EPHuLElqT+dDxLuzJKk9nQ8RSVJ7\nDBFJUmOGiCSpMV/AOMEWvsjRSaokTRp7IpKkxgwRSVJjnQ8RnxORpPZ0PkR8TkSS2tP5EJEktce7\ns4bA6XAlHa3siUiSGjNEJEmNGSKSpMYMEUlSY4aIJKkxQ0SS1JghIklqzOdEOsS3+kqaNBPZE0ly\nRpKPJbkxyTvbrkeStLixhUiS65I8kmTXgvaNSb6VZG+SqwCqak9VXQH8NvCqcdUoSRrMOHsiW4GN\n8xuSHANcA7we2ABcnGRD77M3AV8Abh5jjZKkAYwtRKrqDuCxBc3nAnur6oGqehK4Hriwt/72qno9\ncMmh9pnk8iQ7k+zct2/fqEqXJB1C2wPrU8CD85ZngVck+VXgLcCzOUxPpKq2AFsAZmZmanRlSpIW\n03aILKqqbgdu72fdJJuATevXrx9lSZKkRbR9d9ZDwJp5y9O9tr45KZUktaftELkbOC3JuiQrgYuA\n7S3XJEnq0zhv8d0G3AmcnmQ2yeaqehq4ErgN2APcUFX3Dbhf51iXpJak6sgYj56ZmamdO3eO/DiT\nNIuhT6xLWq4k91TVTNPt276cJUnqsM6HiJezJKk9nQ8R786SpPZM5HMi6s9i4zOOk0gap873RLyc\nJUnt6XyIeDlLktrT+RCRJLXHEJEkNdb5EHFMRJLa0/kQcUxEktrT+RCRJLXHEJEkNdb5EHFMRJLa\n0/kQcUxEktozcIgkOS7JMaMoRpLULUuGSJJnJfndJF9I8gjwb8B/J9md5M+TOLm5JB2l+umJfAn4\nReC9wAuqak1VnQz8CvAV4ENJfm+ENUqSJlQ/b/E9r6qeWthYVY8B/wj8Y5IVQ69sQkzSTIaSNGmW\nDJGDAZLkeOAs4LtV9R+LraP2LQw9Xw0vaZQGmU/kS8DXgRcn2Qdc2uuNtCrJJmDT+vUOzUjSuPV9\nd1ZVzVTV5qp6FfBp4OYkp4yutL7r8hZfSWrJQLf49u7UmgK+A9zKXO9EknSU6vtyVpLvAj8HfA94\nuPdvx4jqkiR1wCBjIr9UVT8YWSWSpM4ZJERmklwCfB/YBXwD2FVV/zeSyiRJE2+QELkO+GNgBXO3\n+r4ZOBPwtihJOkoNEiLfqarP9X7+h1EUI0nqln7enfWpJH8MfCXJn4yhpoH4KnhJak8/t/huBQL8\nAvD7Sb6TZHuSP0vyWyOtrg8+JyJJ7enntSf/AvzLweUkxwJnAGcD5+KlLUk6ai0ZIklSVXVwuaqe\nBr7Z+/fpxdaRJB0d+noVfJJ3J3nR/MYkK5P8epJPAn8wmvIkSZOsn7uzNgKXAtuSnAr8D3NPrj8L\n+CLwl1V17+hKlCRNqn5C5Beq6q+Av+rNG3Ii8KOq+v5oS5MkTbp+QuSWJCczNy3uN5gbC/lGkm/6\nGhRJOrr1c3fWhiTPBjYAv8zc0+oXAmcl+b+qWjfiGiVJE6qvJ9Z778e6N8le4EfMXdL6JeZ6Jppg\ni03v62yHkoaln1t8TwfeALwROAn4Z+AzwOVV9eRoyxsv51OXpMH00xPZA9wLfAi4ybf2SpIO6idE\n3gm8BHgXcHWS/fz0YcNvznsp41AleTNzPaDnAp+oqi+O4jiSpOaWfNiwqv6mqt5dVa+tqpOB1wF/\nCzwJvHWQgyW5LskjSXYtaN+Y5FtJ9ia5qnfcz1XVZcAVwO8MchxJ0ngM8ip4AKpqFpgFbmlwvK3A\n1cCnDjYkOQa4hrlwmgXuTrK9qnb3Vnlf73NJ0oTp57UnQ1NVdwCPLWg+F9hbVQ/0BuqvBy7MnA8B\nt1TV1xbbX5LLk+xMsnPfvn2jLV6S9AxjDZFDmAIenLc822t7N3Ae8LYkVyy2YVVtqaqZqpo56aST\nRl+pJOlnDHw5a1yq6qPAR5daL8kmYNP69c7SK0njNgk9kYeANfOWp3ttfXFSKklqzySEyN3AaUnW\nJVkJXARs73djp8eVpPaMNUSSbAPuBE5PMptkc2+SqyuB25h7sPGGqrqv333aE5Gk9ox1TKSqLj5E\n+83AzeOsRZK0fJNwOWtZvJwlSe3pfIh4OUuS2tP5EJEktccQkSQ11vkQcUxEktrT+RBxTESS2jOx\nrz3R6CycwdHpciU11fmeiJezJKk9nQ8RL2dJUns6HyKSpPYc1WMiC8cGJEmDsSciSWqs8yHiwLok\ntafzIeLAuiS1p/MhIklqjyEiSWrMEJEkNWaISJIa63yIeHeWJLWn8w8bVtUOYMfMzMxlbdfSVYs9\ndOlLGSX1o/M9EUlSewwRSVJjhogkqTFDRJLUmCEiSWrMEJEkNdb5W3wH4fwhkjRcne+J+LChJLWn\n8yHiq+AlqT2dDxFJUnsMEUlSY4aIJKkxQ0SS1JghIklqzBCRJDVmiEiSGjNEJEmNGSKSpMYmMkSS\nnJrkE0lubLsWSdKhje0FjEmuA94IPFJVL5nXvhH4CHAMcG1VfbCqHgA2GyLtWfiySudcl7SYcfZE\ntgIb5zckOQa4Bng9sAG4OMmGMdYkSVqGsYVIVd0BPLag+Vxgb1U9UFVPAtcDF/a7zySXJ9mZZOe+\nffuGWK0kqR9tj4lMAQ/OW54FppKsTvIx4Jwk7z3UxlW1papmqmrmpJNOGnWtkqQFJnJSqqraD1zR\ndh2SpMNruyfyELBm3vJ0r61vTkolSe1pO0TuBk5Lsi7JSuAiYPsgO3BSKklqz9hCJMk24E7g9CSz\nSTZX1dPAlcBtwB7ghqq6b8D92hORpJaMbUykqi4+RPvNwM3L2O8OYMfMzMxlTfchSWqm7ctZkqQO\nm8i7swaRZBOwaf369T/TvvCJa0nS8HW+J+LAuiS1p/MhIklqT+dDxLuzJKk9nQ8RL2dJUns6HyKS\npPYYIpKkxo7YW3w1XE5SJWkxne+JOCYiSe3pfIhIktpjiEiSGjNEJEmNdT5EfNhQktrT+RBxYF2S\n2tP5EJEktccQkSQ1ZohIkhozRCRJjXU+RLw7S5La0/kQ8e4sSWpP50NEktQeQ0SS1JghIklqzBCR\nJDVmiEiSGjNEJEmNpararmFZDk6Pe+zzXnjZ1B99vO1ytASn1ZUmS5J7qmqm6fad74kcfE7kWauO\na7sUSTrqdD5EJEntMUQkSY0ZIpKkxgwRSVJjhogkqTFDRJLUmCEiSWrMEJEkNdb5J9YPSvI4cP8A\nm5wALDUd4nLXOdxnJwKPLrHvSdLP32JSjrGc/Qyybb/rLrVe0889h0Z7nKb7GXS7tr+LTq+q45fY\n96FV1RHxD9gy7PWXu84Sn+1s+282yr9vm8dYzn4G2bbfdZdar+nnnkOTeR4dbd9FR9LlrB0jWH+5\n6wxa0yQbx+8yrGMsZz+DbNvvukutt9zPu2Jcv0fb59FR9V10xFzO6pokO2sZLz2TPIc0DMs9j46k\nnkjXbGm7AHWe55CGYVnnkT0RSVJj9kQkSY0ZIpKkxgwRSVJjhogkqTFDZAIkOS7JJ5N8PMklbdej\nbkpyapJPJLmx7VrUXUne3Psu+vskv7nU+obIiCS5LskjSXYtaN+Y5FtJ9ia5qtf8FuDGqroMeNPY\ni9XEGuQ8qqoHqmpzO5Vqkg14Hn2u9110BfA7S+3bEBmdrcDG+Q1JjgGuAV4PbAAuTrIBmAYe7K12\nYIw1avJtpf/zSDqUrQx+Hr2v9/lhGSIjUlV3AI8taD4X2Nv7f4xPAtcDFwKzzAUJ+L+J5hnwPJIW\nNch5lDkfAm6pqq8ttW+/sMZrip/2OGAuPKaAfwLemuSvOXLek6TRWfQ8SrI6yceAc5K8t53S1CGH\n+j56N3Ae8LYkVyy1k2NHU5sGUVX/C/xh23Wo26pqP3PXsaXGquqjwEf7Xd+eyHg9BKyZtzzda5MG\n4XmkYRjKeWSIjNfdwGlJ1iVZCVwEbG+5JnWP55GGYSjnkSEyIkm2AXcCpyeZTbK5qp4GrgRuA/YA\nN1TVfW3WqcnmeaRhGOV55Ft8JUmN2RORJDVmiEiSGjNEJEmNGSKSpMYMEUlSY4aIJKkxQ0SS1Jgh\nIklqzBCR+pTkN5L83QDrfyzJq5psK3WFISL172zg3gHWfyXwlYbbSp1giEj9Oxt4QZI7kvxXkvOS\nnJDkewdXSHJPr+0M4N+r6sCg2473V5KWxxCR+nc2sK+qXgO8B7ikqh4HnpPk4Nw8XwfOYm7K0Vsb\nbit1hiEi9SHJCmA18Be9phXA93s/fxd4Ye/nF/eWz6cXIoNum+S4JJ9M8vEkl4zmN5KGwxCR+nMG\n8PWq+nFv+SxgV+/nh4FTkrwNeJS5iX2eV1UPD7ptVd0PvAW4saouA940yl9KWi6nx5X6czZzl5sO\nOgu4qffzw8AFzF3CugD4NeBLDbeFuRnmvtn7+QDSBLMnIvXnbOAb85Zfws/2Jn4XeGtVPcri4yH9\nbgswy1yQgP+NasI5KZU0ZEm+Bryiqp5quP1xwNXAE8C/VtVnhlmfNEyGiCSpMbvKkqTGDBFJUmOG\niCSpMUNEktSYISJJaswQkSQ1ZohIkhozRCRJjf0/smdtEsscLw8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1a181401d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fraction of transmitted photons: 0.483\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAERCAYAAAC6kZqPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAE6JJREFUeJzt3X+sX3Wd5/Hni1LswGJnhGFG7y2hTBmG4oKYK5q46uwu\nMxS0YtCZobKbdag0TsQ42X8Wk9lszMYEs26yY8BhK5LqrpZhmIlQLbCbjAQ3i0KRqKVd1wZ0e8vO\n8GtF4sgA5b1/3G+Ha3tv7/d7vvf7Pd/TPh9Jk3s+95zP932vh/vycz7nfE6qCkmSmjih7QIkSd1l\niEiSGjNEJEmNGSKSpMYMEUlSY4aIJKkxQ0SS1JghIklqzBCRJDV2YtsFLJfTTz+9zjrrrLbLkKRO\nefjhh5+uql9tevwxEyJnnXUWu3btarsMSeqUJD8e5ngvZ0mSGjNEJEmNGSKSpMY6PyeSZCOwcd26\ndW2XIqkFL730ErOzs7zwwgttlzLRVq1axfT0NCtXrlzWfnOsvE9kZmamnFiXjj+PP/44p556Kqed\ndhpJ2i5nIlUVzzzzDM8//zxr1679he8lebiqZpr27eUsSZ32wgsvGCBLSMJpp502ktGaISKp8wyQ\npY3qd2SISJIa6/zEeleddf3Xf2H7Rze8u6VKJKk5RyKSpMYciUg6phw+yh9Wm1cJrrnmGr72ta9x\nxhlnsHv37kX3u+eee/j4xz/OwYMH+fCHP8z1118/tho7PxJJsjHJ1ueee67tUiRpWX3oQx/innvu\nOeo+Bw8e5KMf/Sh33303e/bsYfv27ezZs2dMFR4DIVJVO6pqy+rVq9suRdJx6qc//SkXXXQR559/\nPieffDJvetObeNvb3sYrr7wyVL/vfOc7ed3rXnfUfR588EHWrVvH2WefzUknncRVV13FnXfeOdTn\nDsLLWZI0pNe+9rU88sgjPPjgg3zqU5864o/4O97xDp5//vkjjvvMZz7DJZdcMtRnHzhwgDVr1vzD\n9vT0NN/+9reH6nMQhogkLZPdu3dz/vnnH9H+zW9+s4VqxsMQkaRlsmfPHt785jcf0T7KkcjU1BT7\n9+//h+3Z2VmmpqaG6nMQhogkLZMnnniCyy+//Ij2UY5E3vKWt/DDH/6Qxx9/nKmpKW677Ta+8pWv\njOzzDmeISDqmtHlL7qWXXsrmzZvZtm0b73rXu4bub9OmTdx33308/fTTTE9P88lPfpLNmzdz+eWX\nc8stt/CGN7yBE088kRtvvJFLL72UgwcPcs011yx4SW1UXMW3JT6xLi2PvXv3ct5557VdRics9Lty\nFV9JUmsMEUlSYxMbIklOSbIryXvarkXSZDtWLsuP0qh+R2MLkSS3Jnkyye7D2jck+UGSfUnmL/jy\nb4Dbx1WfpG5atWoVzzzzjEFyFIfebLhq1apl73ucd2dtA24EvnSoIckK4Cbgd4BZ4KEkdwFTwB5g\n+X9iSceU6elpZmdneeqpp9ouZaIdesf6chtbiFTV/UnOOqz5YmBfVT0GkOQ24ArgHwGnAOuBnyfZ\nWVVHLEKTZAuwBeDMM88cXfGSJtbKlSuPeG+4xqft50SmgP3ztmeBt1bVdQBJPgQ8vVCAAFTVVmAr\nzN3iO9pSJUmHaztEjqqqti21T5KNwMZ169aNviBJ0i9o++6sA8CaedvTvba+uRS8JLWn7RB5CDgn\nydokJwFXAXe1XJMkqU9ju5yVZDvw28DpSWaBf1dVX0hyHXAvsAK4taoeHbDfib+ctdyv65SkSTHO\nu7M2LdK+E9g5RL87gB0zMzPXNu1DktRM25ezhuY71iWpPZ0PESfWJak9nQ8RSVJ7Oh8iXs6SpPZ0\nPkS8nCVJ7el8iEiS2mOISJIa63yIOCciSe2Z6AUY+3GsPGy40FPtP7rh3S1UIkn96/xIRJLUHkNE\nktRY50PEORFJao9zIiPgqr2SjhedH4lIktpjiEiSGjNEJEmNGSKSpMY6HyLenSVJ7el8iLiKryS1\np/MhIklqjyEiSWqs8w8bHssOf2jRBRklTRpHIpKkxgwRSVJjhogkqbHOh4jPiUhSezofIj4nIknt\n8e6sZeDS75KOV50fiUiS2mOISJIaM0QkSY0ZIpKkxgwRSVJj3p3VIa6lJWnSOBKRJDU2kSGS5Lwk\nNye5I8kftV2PJGlhYwuRJLcmeTLJ7sPaNyT5QZJ9Sa4HqKq9VfUR4PeBt4+rRknSYMY5EtkGbJjf\nkGQFcBNwGbAe2JRkfe977wW+DuwcY42SpAGMLUSq6n7g2cOaLwb2VdVjVfUicBtwRW//u6rqMuDq\nxfpMsiXJriS7nnrqqVGVLklaRNt3Z00B++dtzwJvTfLbwJXAazjKSKSqtgJbAWZmZmp0ZUqSFtJ2\niCyoqu4D7mu5DEnSEtq+O+sAsGbe9nSvrW++T0SS2tP2SOQh4Jwka5kLj6uADw7SQVXtAHbMzMxc\nO4L6juCy75L0qnHe4rsdeAA4N8lsks1V9TJwHXAvsBe4vaoeHbBfRyKS1JKxjUSqatMi7TsZ4jbe\ncY9EJEmvantOZGiORCSpPW3PiQzteB6JLDQ/46KMksap8yMRSVJ7DBFJUmOdDxHnRCSpPZ0Pkara\nUVVbVq9e3XYpknTc6XyISJLa0/kQ8XKWJLWn8yHi5SxJak/nQ0SS1B5DRJLUWOefWE+yEdi4bt26\ntkuZCIc/xe4T7JJGqfMjEedEJKk9nR+JjJrvD5GkxXV+JCJJao8hIklqzBCRJDXW+RDxiXVJak/n\nQ8S7sySpPQOHSJJTkqwYRTGSpG5ZMkSSnJDkg0m+nuRJ4H8B/zfJniT/IYlP+UnScaqfkcg3gN8A\nPgH8elWtqaozgH8CfAv4dJJ/McIaJUkTqp+HDS+pqpcOb6yqZ4G/BP4yycplr0ySNPGWHIkcCpAk\npyZ5e5LfWGwfSdLxZZBlT74BfBf4rSRPAdf0RiOSpONU3yFSVTOHvk7yAWBnkiur6omRVNYnV/GV\npPYMtABjkhOA1wM/Bu5hbnRy7gjq6ltV7QB2zMzMXNtmHZNqoQUkXR5e0nLpO0SS/A3wS8DfAk/0\n/u0YUV2SpA4YZCTym1X105FVIknqnEFCZCbJ1cBPgN3A94DdVfX3I6msBb47RJIGM0iI3Ar8MbAS\nuAB4H3A+4Iy2JB2nBgmRH1fVV3tf/8UoipEkdUs/a2d9KckfA99K8q/HUJMkqSP6WTtrGxDg14B/\nmeTHSe5K8u+T/N5Iq5MkTbQlL2dV1V8Df31oO8mJwHnAhcDFeGlLko5bS4ZIklRVHdquqpeB7/f+\n/deF9pEkHR/6Wgo+yceSnDm/MclJSf5Zki8C/2q5C0vyviSfT/LnSX53ufuXJA2vn7uzNgDXANuT\nnA38P+aeXD8B+G/Af6qqR/r5sCS3Au8BnqyqN85r3wD8KbACuKWqbujdCfbVJL8CfKb3WVoGhz8P\n4zIokprqJ0R+rao+B3yu996Q04GfV9VPGnzeNuBG4EuHGnqv2r0J+B1gFngoyV1Vtae3y5/0vi9J\nmjD9hMjdSc5g7rW432NuLuR7Sb4/6DIoVXV/krMOa74Y2FdVjwEkuQ24Isle4Abg7qr6zkL9JdkC\nbAE488wzF9pFkjRC/dydtT7Ja4D1wD9m7mn1K4ALkvx9Va0dsoYpYP+87VngrcDHgEuA1UnWVdXN\nC9S2FdgKMDMz48S+JI1ZX0+s99bHeiTJPuDnzF3S+k3mRiYjUVWfBT671H6+T0SS2tPPLb7nAu9m\nbkL8V4H/DnwZ2FJVLy5DDQeANfO2p3ttfRnmfSIuuChJw+lnJLIXeAT4NHDnCFbtfQg4J8la5sLj\nKuCD/R7sSESS2tPPcyJ/BPxP4KPA/iR7k9ye5N8med8gH5ZkO/AAcG6S2SSbew8vXgfcy1xg3V5V\nj/bbZ1XtqKotq1evHqQUSdIy6Gdi/T/P304yzasT7O8HvrrQcYv0tWmR9p3Azn77kSRNhoHesQ5Q\nVbPM3UF19/KXMzgvZ0lSe/q5nDXRvJwlSe3pfIhIktoz8OWsSePlrOEtdKuz62lJ6kfnRyJezpKk\n9nQ+RCRJ7el8iCTZmGTrc88913YpknTc6XyIeDlLktrT+RCRJLXHEJEkNdb5W3wH4aq9krS8Oj8S\ncWJdktrT+RBxYl2S2tP5EJEktccQkSQ1ZohIkhrrfIg4sS5J7el8iDixLknt6XyISJLac1w9bKj+\nLfVgpu8bkQSORCRJQzBEJEmNGSKSpMYMEUlSY50PEZ8TkaT2dD5EfE5EktrT+RCRJLXHEJEkNWaI\nSJIaM0QkSY0ZIpKkxo7ZtbOWWvtJkjQ8RyKSpMYMEUlSY4aIJKmxiQyRJGcn+UKSO9quRZK0uLFN\nrCe5FXgP8GRVvXFe+wbgT4EVwC1VdUNVPQZsHiREvn/gOSfTJWnMxjkS2QZsmN+QZAVwE3AZsB7Y\nlGT9GGuSJA1hbCFSVfcDzx7WfDGwr6oeq6oXgduAK/rtM8mWJLuS7Dr4d67iK0nj1vacyBSwf972\nLDCV5LQkNwMXJfnEYgdX1daqmqmqmRUnu4qvJI3bRD5sWFXPAB/pZ98kG4GNJ/7y60dblCTpCG2P\nRA4Aa+ZtT/fa+nbofSInrDplWQuTJC2t7RB5CDgnydokJwFXAXcN0sGhNxu+8sLPRlKgJGlxYwuR\nJNuBB4Bzk8wm2VxVLwPXAfcCe4Hbq+rRQfp1JCJJ7RnbnEhVbVqkfSewc1x1SJKWz0ROrA/CiXVJ\nak/bcyJD83KWJLWn8yEiSWqPl7MkSY11fiTi5SxJak/nQ0SS1B5DRJLUmHMiaqSfd7f86IZ3j6ES\nSW3q/EjEORFJak/nQ0SS1B5DRJLUmHMikqTGOj8ScU5EktrT+RCRJLXHEJEkNWaISJIaM0QkSY15\nd5YkqbHOj0S8O0uS2tP5EJEktccQkSQ1ZohIkhozRCRJjRkikqTGDBFJUmM+J6KR8e2H0rGv8yMR\nnxORpPZ0PkQkSe0xRCRJjRkikqTGDBFJUmOGiCSpMUNEktSYISJJaswQkSQ1NpFPrCc5Bfgc8CJw\nX1V9ueWSJEkLGNtIJMmtSZ5Msvuw9g1JfpBkX5Lre81XAndU1bXAe8dVoyRpMOO8nLUN2DC/IckK\n4CbgMmA9sCnJemAa2N/b7eAYa5QkDWBsIVJV9wPPHtZ8MbCvqh6rqheB24ArgFnmguSoNSbZkmRX\nkl0H/+65UZQtSTqKtifWp3h1xAFz4TEF/BXw/iR/BuxY7OCq2lpVM1U1s+Lk1aOtVJJ0hImcWK+q\nnwF/2M++LgUvSe1peyRyAFgzb3u619Y3l4KXpPa0HSIPAeckWZvkJOAq4K6Wa5Ik9Wmct/huBx4A\nzk0ym2RzVb0MXAfcC+wFbq+qRwfsd2OSra+88LPlL1qSdFRjmxOpqk2LtO8Edg7R7w5gx2tef861\nTfuQJDUzkRPrg3BiXZLa0/acyNCcWJek9nQ+RCRJ7UlVtV3DUA5dzgKuZm5yvl+rgaUecx92n6N9\n73Tg6SX6niT9/C4m5TOG6WeQY/vdd6n9mn7fc2i0n9O0n0GPa/tv0blVdeoSfS+uqo6Jf8DW5d5/\n2H2W+N6utn9no/z9tvkZw/QzyLH97rvUfk2/7zk0mefR8fa36Fi6nLXo8ihD7D/sPoPWNMnG8bMs\n12cM088gx/a771L7Dfv9rhjXz9H2eXRc/S3q/OWsrkqyq6pm2q5D3eU5pOUw7Hl0LI1EumZr2wWo\n8zyHtByGOo8ciUiSGnMkIklqzBCRJDVmiEiSGjNEJEmNGSITIMkpSb6Y5PNJrm67HnVTkrOTfCHJ\nHW3Xou5K8r7e36I/T/K7S+1viIxIkluTPJlk92HtG5L8IMm+JNf3mq8E7qiqa4H3jr1YTaxBzqOq\neqyqNrdTqSbZgOfRV3t/iz4C/MFSfRsio7MN2DC/IckK4CbgMmA9sCnJeuZeC7y/t9vBMdaoybeN\n/s8jaTHbGPw8+pPe94/KEBmRqrofePaw5ouBfb3/x/gicBtwBTDLXJCA/5tongHPI2lBg5xHmfNp\n4O6q+s5SffsHa7ymeHXEAXPhMQX8FfD+JH/GsbNOkkZnwfMoyWlJbgYuSvKJdkpThyz29+hjwCXA\nB5J8ZKlOOv9mw2NBVf0M+MO261C3VdUzzF3Hlhqrqs8Cn+13f0ci43UAWDNve7rXJg3C80jLYVnO\nI0NkvB4CzkmyNslJwFXAXS3XpO7xPNJyWJbzyBAZkSTbgQeAc5PMJtlcVS8D1wH3MvcWxtur6tE2\n69Rk8zzSchjleeQqvpKkxhyJSJIaM0QkSY0ZIpKkxgwRSVJjhogkqTFDRJLUmCEiSWrMEJEkNWaI\nSH1K8s+T/JcB9r85ydubHCt1hSEi9e9C4JEB9n8b8K2Gx0qdYIhI/bsQ+PUk9yf5P0kuSbI6yd8e\n2iHJw72284D/XVUHBz12vD+SNBxDROrfhcBTVfVO4OPA1VX1HHBykkPv5vkucAFzrxy9p+GxUmcY\nIlIfkqwETgP+Y69pJfCT3td/A7y+9/Vv9bYvpRcigx6b5JQkX0zy+SRXj+YnkpaHISL15zzgu1X1\nSm/7AmB37+sngDck+QDwNHMv9vnlqnpi0GOr6ofAlcAdVXUt8N5R/lDSsHw9rtSfC5m73HTIBcCd\nva+fAC5n7hLW5cA/Bb7R8FiYe8Pc93tfH0SaYI5EpP5cCHxv3vYb+cXRxAeB91fV0yw8H9LvsQCz\nzAUJ+N+oJpwvpZKWWZLvAG+tqpcaHn8KcCPwAvA/qurLy1mftJwMEUlSYw6VJUmNGSKSpMYMEUlS\nY4aIJKkxQ0SS1JghIklqzBCRJDVmiEiSGvv/xlKU1VhDnicAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1a0ea42b10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fraction of transmitted photons: 0.482\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAERCAYAAAC6kZqPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEy9JREFUeJzt3X+MXXWZx/H3Q2npgtAmLSjMlLS1iC0sFDMCiQtCFpZS\nKCiw2MpuiCBdEjFsNtlYE7PG7D+adROXgLpVat1daddFIzSUIi4QNLERkKKlXZYGxE6LUupSiKFY\nhmf/mKkM02nn3u+de889nfcraTL3e88595nxcD9+f5xzIjORJKnEEVUXIEmqL0NEklTMEJEkFTNE\nJEnFDBFJUjFDRJJUzBCRJBUzRCRJxQwRSVKxI6suYLzMnDkzZ8+eXXUZklQrTzzxxMuZeXzp/rUP\nkYhYAiyZN28ejz/+eNXlSFKtRMQLrexf++GszFyXmcunTZtWdSmSNOHUPkQkSdUxRCRJxQ6rORFJ\nE8++ffvo7+9n7969VZfS1aZOnUpvby+TJ08e1+PG4fI8kb6+vnRiXZp4nn/+eY499lhmzJhBRFRd\nTlfKTHbv3s1rr73GnDlz3vFeRDyRmX2lx3Y4S1Kt7d271wAZQ0QwY8aMtvTWDBFJtWeAjK1df6Pa\nh0hELImIlXv27Km6FEmacGo/sZ6Z64B1fX19N1VdSzNmr7jvHa9/9cXLKqpEksrVviciSapO7Xsi\nkjTcyF5+q6oaJdi7dy/nn38+b7zxBm+++SbXXHMNX/jCF0bddsOGDdx6660MDAzwyU9+khUrVnSs\nTnsiktSFjjrqKB566CGeeuopNm3axIYNG9i4ceMB2w0MDPCpT32K+++/ny1btrBmzRq2bNnSsTpr\nHyJOrEuq2quvvspZZ53FaaedxtFHH83ChQs599xzeeutt4qPGRG8613vAgYvqNy3b9+oK6x+9rOf\nMW/ePObOncuUKVNYunQp99xzT/HnNqv2IeINGCVV7bjjjuPJJ5/kW9/6FhdffDGbNm1i48aNHHHE\n4Ffseeedx8KFCw/496Mf/eiQxx0YGGDhwoWccMIJXHzxxZxzzjkHbLNjxw5mzZr1x9e9vb3s2LFj\nfH/BQ3BORJLGyebNmznttNMOaP/xj39cdLxJkyaxadMmXnnlFT760Y+yefNmTj/99FbLHFeGSAeM\n90SfpO60ZcsWPvCBDxzQft555/Haa68d0P7lL3+Ziy66aMzjTp8+nQsvvJANGzYcECI9PT1s3779\nj6/7+/vp6ekpqL6MISJJ42Tnzp0sXrz4gPaSnsiuXbuYPHky06dP5/XXX+fBBx/kM5/5zAHbffCD\nH+TZZ5/l+eefp6enh7Vr13LXXXcV1V/CEJF0WKnywt1LLrmEG2+8kdWrV/PhD3+4pWO9+OKLXH/9\n9QwMDPDWW29x7bXXcvnllwOwePFivvnNb3LSSSdx5JFHcvvtt3PJJZcwMDDADTfcMOqQWrt4F98O\naGQ4yyvWpTJbt25l/vz5VZdRC6P9rVq9i2/teyKHy/NERgsag0VSt3OJrySpWO1DRJIOl2H5dmrX\n38gQkVRrU6dOZffu3QbJIex/suHUqVPH/di1nxPpRl4XInVOb28v/f397Nq1q+pSutr+Z6yPN0NE\nUq1Nnjz5gOeGq3MczpIkFbMn0sV8+qGkbmdPRJJUrPYh4vNEJKk6tQ8RLzaUpOrUPkQkSdUxRCRJ\nxVydVSOu1pLUbeyJSJKK2RMZB97mRNJEZU9EklTMEJEkFTNEJEnFDBFJUjFDRJJUzBCRJBXryiW+\nETEfuBWYCfx3Zn6t4pK60mhLi70AUVIndawnEhGrIuKliNg8on1RRDwTEdsiYgVAZm7NzJuBa4EP\ndapGSVJzOjmctRpYNLwhIiYBdwCXAguAZRGxYOi9K4D7gPUdrFGS1ISODWdl5qMRMXtE89nAtsx8\nDiAi1gJXAlsy817g3oi4D7hrtGNGxHJgOcDJJ5/cpsrfyavTJeltVc+J9ADbh73uB86JiAuAq4Cj\nOERPJDNXAisB+vr6sn1lSpJGU3WIjCozHwEeqbgMSdIYql7iuwOYNex171Bbw3w8riRVp+oQeQw4\nJSLmRMQUYClwbzMH8PG4klSdjg1nRcQa4AJgZkT0A5/PzDsj4hbgAWASsCozn+5UTYcjH1wlqZM6\nuTpr2UHa19PCMt6IWAIsmTdvXukhJEmFqh7OapnDWZJUndqHiCSpOrUPEVdnSVJ1ah8iDmdJUnW6\n8mLDbuJtTiTp4GofIq7OOjRvFy+pnRzOkiQVq32ISJKqY4hIkooZIpKkYrUPEa8TkaTq1D5EnFiX\npOrUPkQkSdWp/XUiap63i5c0XgyRYbw6XZKaU/vhLCfWJak6tQ8RJ9YlqTq1DxFJUnUMEUlSMUNE\nklTMEJEkFat9iLg6S5KqU/sQcXWWJFWn9iEiSaqOISJJKuZtT+Rz2CUVsyciSSpmiEiSik3o4Szv\n2itJrbEnIkkqZohIkorVPkS8Yl2SqlP7EPGKdUmqzoSeWNfB+Rx2SY2ofU9EklQdQ0SSVMwQkSQV\nM0QkScUMEUlSMUNEklRsQi3x9V5ZkjS+7IlIkooZIpKkYoaIJKlY03MiEXEMsDczB9pQz/DP+Qhw\nGXAccGdm/rCdn6dD8zYokkYzZk8kIo6IiI9HxH0R8RLwP8CLEbElIv4pIuY1+mERsSoiXoqIzSPa\nF0XEMxGxLSJWAGTmDzLzJuBm4GPN/VqSpE5oZDjrYeC9wGeB92TmrMw8AfgzYCPwpYj4qwY/bzWw\naHhDREwC7gAuBRYAyyJiwbBNPjf0viSpyzQynHVRZu4b2ZiZvwO+B3wvIiY38mGZ+WhEzB7RfDaw\nLTOfA4iItcCVEbEV+CJwf2b+fLTjRcRyYDnAySef3EgJkqRxNGZPZH+ARMSxEfGhiHjvwbYp1ANs\nH/a6f6jt08BFwDURcfNBaluZmX2Z2Xf88ce3UIIkqUQzE+sPA08B74+IXcANQ72RtsjM24Db2nV8\nSVLrGg6RzOzb/3NEXAOsj4irMnNnizXsAGYNe9071NaQiFgCLJk3r+H5fUnSOGlqiW9EHAGcCLwA\nbGCwd3JqizU8BpwSEXMYDI+lwMcb3Tkz1wHr+vr6bmqxDjVhtFvIuOxXmngaDpGI+A3wJ8BvgZ1D\n/9Y182ERsQa4AJgZEf3A5zPzzoi4BXgAmASsysynmzimPRFJqkgzPZH3ZearrXxYZi47SPt6YH3h\nMe2JSFJFmgmRvoi4DngF2Az8AticmW+0pbIWecdeSWq/Zu6dtYrB4auNwFzgH4CGh53aJSKWRMTK\nPXv2VF2KJE04zfREXsjMHwz9/F/tKKaEw1mSVJ1G7p31bxHxt8DGiPi7DtQkSaqJRoazVgMBvBv4\n64h4ISLujYh/jIi/bGt1kqSuNuZwVmY+BDy0/3VEHAnMB85k8L5XlQ5tucRXkqrTyHBWDH+dmW9m\n5i8z8z8y8+9H26aTMnNdZi6fNm1aVSVI0oQVmXnoDSIeYfBuvfdk5q+HtU9h8Hbw1wMPZ+bq9pU5\ntqNOPCVPvP4rVZagEbyCXep+EfHE8NtaNauR1VmLgBuANRExF/g/Bq9cPwL4IfCVzHyytABJUn01\nEiLvzsyvAl8dem7ITOD1zHylvaU1Zv+cyJHTT6y6FEmacBpZnXV/RLwcET8B/gX4CHBaRBzX3tIa\ns39O5Iipx1RdiiRNOI2szloQEUcx+OjaPwXOAK4EzoiINzJzTptrlCR1qYauWB+6P9aTEbENeJ3B\nIa33MXj/LEnSBDVmiETEqcBlwOXA8cCDwHeA5Zn5h/aWJ0nqZo30RLYCTwJfYnCZb1fdtdeJdUmq\nTiPXifwNcDqDcyHzgd3AL/f/G3ZTxkp5nUg9eO2I1F3afp1IZv7riA/s5e0J9quBrggRSVLnNfWM\ndYDM7Af6gfvHvxxJUp0081AqSZLewRCRJBUzRCRJxZqeE+k2LvGVpOrUvifivbMkqTq1DxFJUnUM\nEUlSMUNEklTMEJEkFTNEJEnFar/EV/Uye8V973jtDRmlerMnIkkqZohIkorVfjjLK9YlqTq174l4\nxbokVaf2ISJJqo4hIkkqZohIkooZIpKkYoaIJKmYISJJKlb760RUb94GRao3eyKSpGKGiCSpmMNZ\n6iojh7fAIS6pm9kTkSQV68oQiYi5EXFnRNxddS2SpIPrWIhExKqIeCkiNo9oXxQRz0TEtohYAZCZ\nz2XmjZ2qTZJUppM9kdXAouENETEJuAO4FFgALIuIBR2sSZLUgo5NrGfmoxExe0Tz2cC2zHwOICLW\nAlcCWxo5ZkQsB5YDTDru+HGrVd3Fa0mk7lX1nEgPsH3Y636gJyJmRMTXgbMi4rMH2zkzV2ZmX2b2\nTTp6WrtrlSSN0JVLfDNzN3BzI9v6ZENJqk7VPZEdwKxhr3uH2hrmkw0lqTpVh8hjwCkRMScipgBL\ngXsrrkmS1KCODWdFxBrgAmBmRPQDn8/MOyPiFuABYBKwKjOfbvK4DmdNMI1c1e5kvNQZnVydtewg\n7euB9S0cdx2w7qgTT7mp9BiSpDJVD2dJkmqsK1dnNcPhLEmqTu17Iq7OkqTq1D5EJEnVMUQkScWc\nE5EkFat9T8Q5EUmqTu1DRJJUHUNEklTMORFJUrHa90ScE5Gk6tQ+RCRJ1TFEJEnFDBFJUjEn1jVh\n+cwRqXW174k4sS5J1al9iEiSqmOISJKKGSKSpGKGiCSpmCEiSSrmEl9JUrHa90Rc4itJ1al9iEiS\nqmOISJKKGSKSpGKGiCSpmCEiSSpmiEiSihkikqRiXmyow8LIZ4NI6oza90S82FCSqlP7EJEkVccQ\nkSQVM0QkScUMEUlSMUNEklTMEJEkFTNEJEnFDBFJUjFDRJJUzBCRJBUzRCRJxbryBowRcQzwVeAP\nwCOZ+Z2KS5IkjaJjPZGIWBURL0XE5hHtiyLimYjYFhErhpqvAu7OzJuAKzpVoySpOZ0czloNLBre\nEBGTgDuAS4EFwLKIWAD0AtuHNhvoYI2SpCZ0bDgrMx+NiNkjms8GtmXmcwARsRa4EuhnMEg2cYig\ni4jlwHKASccdP/5Fa0IpeSbJr754WRsqkeqj6on1Ht7uccBgePQA3weujoivAesOtnNmrszMvszs\nm3T0tPZWKkk6QFdOrGfm74FPVF2HJOnQqg6RHcCsYa97h9oa5uNxJak6VQ9nPQacEhFzImIKsBS4\nt5kD+HhcSapOJ5f4rgF+CpwaEf0RcWNmvgncAjwAbAW+m5lPd6omSVJrOrk6a9lB2tcD60uP63CW\nJFWn6uGsljmcJUnVqX2ISJKqU/XqrJY5nCVJ1YnMrLqGcRERe4Bnm9hlGrCnzdsc6r2ZwMtjHLub\nNPK36JbPaOU4zezb6LZjbVf6vudQez+n9DjN7lf1d9GpmXnsGMc+uMw8LP4BK8d7+1a3GeO9x6v+\nm7Xz71vlZ7RynGb2bXTbsbYrfd9zqDvPo4n2XXQ4zYkc9PYoLWzf6jbN1tTNOvG7jNdntHKcZvZt\ndNuxtmv1/bro1O9R9Xk0ob6LDpvhrLqJiMczs6/qOlRfnkMaD62eR4dTT6RuVlZdgGrPc0jjoaXz\nyJ6IJKmYPRFJUjFDRJJUzBCRJBUzRCRJxQyRLhARx0TEtyPiGxFxXdX1qJ4iYm5E3BkRd1ddi+or\nIj4y9F30nxHxF2Ntb4i0SUSsioiXImLziPZFEfFMRGyLiBVDzVcBd2fmTcAVHS9WXauZ8ygzn8vM\nG6upVN2syfPoB0PfRTcDHxvr2IZI+6wGFg1viIhJwB3ApcACYFlELGDwscDbhzYb6GCN6n6rafw8\nkg5mNc2fR58bev+QDJE2ycxHgd+NaD4b2Db0/xj/AKwFrgT6GQwS8H8TDdPkeSSNqpnzKAZ9Cbg/\nM38+1rH9wuqsHt7uccBgePQA3weujoivcfjcJ0ntM+p5FBEzIuLrwFkR8dlqSlONHOz76NPARcA1\nEXHzWAep/fNEDgeZ+XvgE1XXoXrLzN0MjmNLxTLzNuC2Rre3J9JZO4BZw173DrVJzfA80ngYl/PI\nEOmsx4BTImJOREwBlgL3VlyT6sfzSONhXM4jQ6RNImIN8FPg1Ijoj4gbM/NN4BbgAWAr8N3MfLrK\nOtXdPI80Htp5HnkXX0lSMXsikqRihogkqZghIkkqZohIkooZIpKkYoaIJKmYISJJKmaISJKKGSJS\ngyLizyPi35vY/usR8aGSfaW6MESkxp0JPNnE9ucCGwv3lWrBEJEadybwnoh4NCJ+HREXRcS0iPjt\n/g0i4omhtvnA/2bmQLP7dvZXklpjiEiNOxPYlZnnA7cC12XmHuDoiNj/bJ6ngDMYfOTohsJ9pdow\nRKQGRMRkYAbwz0NNk4FXhn7+DXDi0M/vH3p9CUMh0uy+EXFMRHw7Ir4REde15zeSxochIjVmPvBU\nZr419PoMYPPQzzuBkyLiGuBlBh/sMz0zdza7b2Y+C1wF3J2ZNwFXtPOXklrl43GlxpzJ4HDTfmcA\n9wz9vBNYzOAQ1mLgQuDhwn1h8Alzvxz6eQCpi9kTkRpzJvCLYa9P5529iY8DV2fmy4w+H9LovgD9\nDAYJ+N+oupwPpZLGWUT8HDgnM/cV7n8McDuwF/hJZn5nPOuTxpMhIkkqZldZklTMEJEkFTNEJEnF\nDBFJUjFDRJJUzBCRJBUzRCRJxQwRSVKx/wfzTFC79TPDEwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10b0d5210>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fraction of transmitted photons: 0.259\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAERCAYAAAC6kZqPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAExZJREFUeJzt3X+MXlWdx/HPh6GlC5Y2W4kLMzVtLWBbFkRHIDGi7FZb\nCuVnF1vJhkhtgymGzf5jTdzdmI0JZt1ECShbsam62sqiC21oy7IBg2ZBWmxxh3bZTlC3T4tSyjIY\nl0I7/e4f84wM05nOfc7z4zx35v1KSOaeufc83xku8+Gcc384IgQAQIpTchcAACgvQgQAkIwQAQAk\nI0QAAMkIEQBAMkIEAJCMEAEAJCNEAADJCBEAQLJTcxdQL9tLJS2dOnXqqvPOOy93OQBQKs8888zL\nEXFW6vEeL4896e7ujp07d+YuAwBKxfYzEdGdejzTWQCAZIQIACBZ6UPE9lLb6/r6+nKXAgATTukX\n1iNii6Qt3d3dq3LXAqD1jh49qkqloiNHjuQupa1NmTJFXV1dmjRpUkP7LX2IAJjYKpWKpk6dqlmz\nZsl27nLaUkTo8OHDqlQqmj17dkP7Lv10FoCJ7ciRI5oxYwYBchK2NWPGjKaM1kofIqyJACBAxtas\n31HpQyQitkTE6mnTpuUuBQAmHNZEMpm19uG3bf/qzqsyVQIA6Uo/EgEA5MNIBMC4MnyUX69cswRH\njhzR5ZdfrjfeeEPHjh3TsmXL9MUvfnHEfbdv36477rhD/f39+vSnP621a9e2rM7Sj0RYWAcwHp12\n2ml67LHH9Oyzz2r37t3avn27nnrqqRP26+/v15o1a7Rt2zbt2bNHGzdu1J49e1pWZ+lDhIV1ALm9\n9tpruvjii7VgwQKdfvrpet/73qfLLrtMx48fT+7Ttt7xjndIGrih8ujRoyNeYfX0009r7ty5mjNn\njiZPnqzly5froYceSv7cWjGdBQB1OvPMM7Vr1y49/fTT+tKXvnTCH/EPf/jD+t3vfnfCcV/5yle0\ncOHCUfvt7+/XBz7wAfX29mrNmjW69NJLT9jnwIEDmjlz5h+2u7q69LOf/ayOn6Y2hAgANEhPT48W\nLFhwQvtPfvKTpP46Ojq0e/duvfrqq7r++uvV09OjCy64oN4yG4oQaYFGL/QBaE979uzR+9///hPa\nU0cig6ZPn64rrrhC27dvPyFEOjs7tX///j9sVyoVdXZ2JlSfhhABgAY5ePCglixZckJ7ykjk0KFD\nmjRpkqZPn67XX39djz76qD73uc+dsN8HP/hB7du3T7/85S/V2dmpTZs26fvf/35S/SkIEQDjSs4b\ndxctWqSVK1dqw4YN+shHPlJXXy+++KJuueUW9ff36/jx47rpppt09dVXS5KWLFmi++67T+ecc45O\nPfVU3X333Vq0aJH6+/t16623jjil1iylfz3u4DvW586du2rfvn25yxlRkeks7lgH0uzdu1fz5s3L\nXUYpjPS7mvCvx+USXwDIh+msNjHSaIXRCYB2V/qRCACUfVq+FZr1OyJEAJTalClTdPjwYYLkJAbf\nbDhlypSG9810VhNwXwjQOl1dXapUKjp06FDuUtra4DvWG40QAVBqkyZNavh7w1EcIdLGeHEVgHbH\nmggAIBkhAgBIVvoQ4aVUAJBP6UOEO9YBIJ/ShwgAIB9CBACQjEt8S4RLfgG0G0KkAbhDHcBExXQW\nACAZIQIASEaIAACSESIAgGSECAAgGVdnlRiv1AWQGyMRAECytgwR2/Ns32v7AdufyV0PAGBkLQsR\n2+ttv2S7Z1j7YtvP2+61vVaSImJvRNwm6SZJH2pVjQCA2rRyTWSDpLslfWewwXaHpHskfUxSRdIO\n25sjYo/tayR9RtJ3W1jjmLg7HQDe0rKRSEQ8IemVYc2XSOqNiBci4k1JmyRdW91/c0RcKenm0fq0\nvdr2Tts7Dx061KzSAQCjyH11Vqek/UO2K5Iutf1RSTdIOk3S1tEOjoh1ktZJUnd3dzSvTADASHKH\nyIgi4seSflxkX9tLJS2dO3duM0sqDZ70C6CVcl+ddUDSzCHbXdW2wnizIQDkkztEdkg61/Zs25Ml\nLZe0OXNNAICCWnmJ70ZJT0o633bF9sqIOCbpdkmPSNor6f6IeK7GfpfaXtfX19f4ogEAJ9WyNZGI\nWDFK+1adZPG8QL9bJG3p7u5eldoHACBN7uksAECJteXVWbXg6qyT4yGNAJqp9CHS7Oks7lAHgNEx\nnQUASEaIAACSlT5EuMQXAPIpfYhwxzoA5FP6EAEA5EOIAACSlT5EWBMBgHxKHyKsiQBAPqUPEQBA\nPqW/Yx2148VVABqFEBmCR5wAQG1KP53FwjoA5FP6EGFhHQDyKX2IAADyIUQAAMkIEQBAMkIEAJCM\nS3zBK3QBJCv9SIRLfAEgn9KHCJf4AkA+E3o6izvUAaA+pR+JAADyIUQAAMkm9HQWRseTfgEUwUgE\nAJCMEAEAJCNEAADJSh8i3GwIAPmUPkS42RAA8il9iAAA8iFEAADJCBEAQLIJdbMhz8oCgMZiJAIA\nSDahRiJIx2NQAIyEkQgAIBkhAgBIRogAAJIRIgCAZIQIACAZIQIASFZziNg+w3ZHM4oZ9jnX2f6m\n7R/Y/nizPw8AULsxQ8T2KbY/afth2y9J+i9JL9reY/sfbM8t+mG219t+yXbPsPbFtp+33Wt7rSRF\nxIMRsUrSbZI+UduPBQBohSI3Gz4u6d8lfV5ST0QclyTbfyzpCklftv2vEfHPBfraIOluSd8ZbKiO\nau6R9DFJFUk7bG+OiD3VXb5Q/T7ayEiPkOEGRGDiKRIiCyPi6PDGiHhF0g8l/dD2pCIfFhFP2J41\nrPkSSb0R8YIk2d4k6VrbeyXdKWlbRPx8pP5sr5a0WpLe/e53FykBANBAY05nDQaI7am2P2T7PaPt\nk6hT0v4h25Vq22clLZS0zPZto9S2LiK6I6L7rLPOqqMEAECKWp6d9bikZyW91/YhSbdWRyNNERF3\nSbprrP1sL5W0dO7cwkszaBKerwVMPIVDJCK6B7+2vUzSVts3RMTBOms4IGnmkO2ualvRurZI2tLd\n3b1qaDuPfQeA5qvpEt/qlVqdkn4tabsGRif12iHpXNuzbU+WtFzS5gb0CwBossIjEdu/kfRHkn4r\n6WD1ny21fJjtjZI+KumdtiuS/i4ivmX7dkmPSOqQtD4inquhT6azACCTWtZEzouI1+r5sIhYMUr7\nVklbE/sccToLANB8tYRIt+2bJb0qqUfSLzRw38gbTakMAND2alkTWa+B6aunJM2R9LeSCk87NYvt\npbbX9fX15S4FACacWkYiv46IB6tf/0sziknBdBYA5FPk2Vnfsf1Xkp6y/dctqAkAUBJFprM2SLKk\nd0n6S9u/tr3Z9t/b/oumVgcAaGtjTmdFxGOSHhvctn2qpHmSLtLAc6+yTm1xiS8A5OOIOPkOtmOM\nnYrs02ynnX1unH3LV3OWgGF47AnQ/mw/M/SJJLUqMp31uO3P2n7bY3JtT7b9Z7a/LemW1AIAAOVV\n5OqsxZJulbTR9hxJ/6uBO9dPkfRvkr4aEbuaVyLKqsjzyxitAOVWJETeFRFfl/T16ntD3inp9Yh4\ntbmlFTO4JnLq9LNzlwIAE06R6axttl+2/VNJX5N0naQFts9sbmnFRMSWiFh9ypQzcpcCABNOkauz\n5ts+TdJ8SX8q6UJJ10q60PYbETG7yTUCANpUoTvWq8/H2mW7V9LrGpjSOk8Dz88CAExQY4aI7fMl\nXSXpaklnSXpU0vckrY6IN5tbHgCgnRUZieyVtEvSlyU91G5P7WVhHQDyKbKw/hlJ/yFpjaT9tvfa\nvt/239i+rrnljY2FdQDIp8jC+j8N3bbdpbcW2G+U9OBIxwEAxr9aHgUvSYqIiqSKpG2NLwcAUCa1\nvJQKAIC3IUQAAMlqns4Cmmmk523xfC2gfZU+RLjEFwDyKf10Fpf4AkA+pQ8RAEA+hAgAIBkhAgBI\nVvqFdZRbkbcfAmhfjEQAAMkIEQBAMkIEAJCs9Gsi3GwIAPmUfiTCzYYAkE/pQwQAkA8hAgBIRogA\nAJIRIgCAZIQIACAZIQIASEaIAACSESIAgGSECAAgGSECAEhGiAAAkrVliNieY/tbth/IXQsAYHQt\ne4qv7fWSrpb0UkRcMKR9saSvSeqQdF9E3BkRL0haSYhAOvHth7+686pMlQAYrpUjkQ2SFg9tsN0h\n6R5JV0qaL2mF7fktrAkAUIeWhUhEPCHplWHNl0jqjYgXIuJNSZskXVu0T9urbe+0vbP///oaWC0A\noIjcayKdkvYP2a5I6rQ9w/a9ki62/fnRDo6IdRHRHRHdHadPa3atAIBh2vLNhhFxWNJtuesAAJxc\n7hA5IGnmkO2ualthvB534hm+0D4SFt+B1sg9nbVD0rm2Z9ueLGm5pM21dMDrcQEgn5aFiO2Nkp6U\ndL7tiu2VEXFM0u2SHpG0V9L9EfFcq2oCANSnZdNZEbFilPatkram9st0FgDkk3s6q25MZwFAPqUP\nEQBAPrmvzqob01kAkE/pRyJMZwFAPqUPEQBAPkxnAQCSlX4kwnQWAORT+hABAORDiAAAkhEiAIBk\nLKwDAJKVfiTCwjoA5FP6EAEA5EOIAACSESIAgGQsrGNc4hW6QGuUfiTCwjoA5FP6EAEA5EOIAACS\nESIAgGSECAAgGSECAEjGJb4AgGSlH4lwiS8A5FP6EAEA5EOIAACSESIAgGSECAAgGSECAEhGiAAA\nkhEiAIBk3GwIAEhW+pEINxsCQD6lDxEAQD6ECAAgGSECAEhGiAAAkhEiAIBkhAgAIBkhAgBIRogA\nAJIRIgCAZIQIACAZIQIASNaWD2C0fYakr0t6U9KPI+J7mUsCAIygZSMR2+ttv2S7Z1j7YtvP2+61\nvbbafIOkByJilaRrWlUjAKA2rZzO2iBp8dAG2x2S7pF0paT5klbYni+pS9L+6m79LawRAFCDlk1n\nRcQTtmcNa75EUm9EvCBJtjdJulZSRQNBslsnCTrbqyWtlqSOM89qfNEY12atfXjMfX5151U19zHW\nMcB4knthvVNvjTikgfDolPQjSTfa/oakLaMdHBHrIqI7Iro7Tp/W3EoBACdoy4X1iPi9pE/lrgMA\ncHK5Q+SApJlDtruqbYXxelwAyCf3dNYOSefanm17sqTlkjbX0gGvxwWAfFp5ie9GSU9KOt92xfbK\niDgm6XZJj0jaK+n+iHiuVTUBAOrTyquzVozSvlXS1tR+mc4CgHxyT2fVjeksAMin9CECAMgn99VZ\ndWM6CwDycUTkrqEhbPdJ2lfDIdMk9TV5n5N9752SXh6j73ZS5HfRLp9RTz+1HFt037H2S/0+51Bz\nPye1n1qPy/236PyImDpG36OLiHHxj6R1jd6/3n3G+N7O3L+zZv5+c35GPf3UcmzRfcfaL/X7nEPt\neR5NtL9F42lNZNTHo9Sxf7371FpTO2vFz9Koz6inn1qOLbrvWPvV+/2yaNXPkfs8mlB/i8bNdFbZ\n2N4ZEd2560B5cQ6hEeo9j8bTSKRs1uUuAKXHOYRGqOs8YiQCAEjGSAQAkIwQAQAkI0QAAMkIEQBA\nMkKkDdg+w/a3bX/T9s2560E52Z5j+1u2H8hdC8rL9nXVv0U/sP3xsfYnRJrE9nrbL9nuGda+2Pbz\ntnttr6023yDpgYhYJemalheLtlXLeRQRL0TEyjyVop3VeB49WP1bdJukT4zVNyHSPBskLR7aYLtD\n0j2SrpQ0X9IK2/M18Frg/dXd+ltYI9rfBhU/j4DRbFDt59EXqt8/KUKkSSLiCUmvDGu+RFJv9f8Y\n35S0SdK1kioaCBKJfycYosbzCBhRLeeRB3xZ0raI+PlYffMHq7U69daIQxoIj05JP5J0o+1vaPw8\nJwnNM+J5ZHuG7XslXWz783lKQ4mM9vfos5IWSlpm+7axOin9+0TGg4j4vaRP5a4D5RYRhzUwjw0k\ni4i7JN1VdH9GIq11QNLMIdtd1TagFpxHaISGnEeESGvtkHSu7dm2J0taLmlz5ppQPpxHaISGnEeE\nSJPY3ijpSUnn267YXhkRxyTdLukRSXsl3R8Rz+WsE+2N8wiN0MzziKf4AgCSMRIBACQjRAAAyQgR\nAEAyQgQAkIwQAQAkI0QAAMkIEQBAMkIEAJCMEAEKsv3ntr9bw/732v5QyrFAWRAiQHEXSdpVw/6X\nSXoq8VigFAgRoLiLJP2J7Sds/4/thban2f7t4A62n6m2zZP03xHRX+uxrf2RgPoQIkBxF0k6FBGX\nS7pD0s0R0SfpdNuD7+Z5VtKFGnjl6PbEY4HSIESAAmxPkjRD0j9WmyZJerX69W8knV39+r3V7UWq\nhkitx9o+w/a3bX/T9s3N+YmAxiBEgGLmSXo2Io5Xty+U1FP9+qCkc2wvk/SyBl7sMz0iDtZ6bETs\nk3SDpAciYpWka5r5QwH14vW4QDEXaWC6adCFkh6qfn1Q0hINTGEtkXSFpMcTj5UG3jD3n9Wv+wW0\nMUYiQDEXSfrFkO0L9PbRxCcl3RgRL2vk9ZCix0pSRQNBIvHfKNocL6UCGsz2zyVdGhFHE48/Q9Ld\nko5I+mlEfK+R9QGNRIgAAJIxVAYAJCNEAADJCBEAQDJCBACQjBABACQjRAAAyQgRAEAyQgQAkOz/\nAcOySqYd2qtxAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1a18220790>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fraction of transmitted photons: 0.261\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAERCAYAAAC6kZqPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEktJREFUeJzt3X2MXNV9xvHnsbHZYoERNiRm19Ga2iHYKS/RmiC1IX2x\nYmO6EEGq2KHNix1bVDgiqlTJkSJVURUpqKkUIUioE5DTluCmEAFLDAQpIFIpwby4gQVKcSHEYzcB\nO42JUBxg+fWPGcN6veu9c2bu2873I1nauXPu3d8ud+fh3HPuuY4IAQCQYlbZBQAA6osQAQAkI0QA\nAMkIEQBAMkIEAJCMEAEAJCNEAADJCBEAQDJCBACQ7ISyC+iWhQsXxuDgYNllAECtPP744wci4vTU\n/WdMiAwODuqxxx4ruwwAqBXbL3WyP5ezAADJCBEAQLLah4jtYdvbDh06VHYpANBzaj8mEhEjkkaG\nhoY2lV0LgGK88cYbajQaOnz4cNml1EZfX58GBgY0Z86crh639iECoPc0Gg2dfPLJGhwclO2yy6m8\niNDBgwfVaDS0ZMmSrh679pezAPSew4cPa8GCBQRIRra1YMGCXHpuhAiAWiJA2pPX76v2IcLAOgCU\np/ZjIjNlYH1w6/eP2fazr1x63DYT3weAotW+JwIAKE/teyJ1NVnPI6XNdPvQW0EvSPlbOZ6i/242\nbNige+65R2eccYZGR0ff3n7ffffp2muv1djYmD772c9q69atk+6ftV0eCJEa6/YfDoByfPrTn9aW\nLVv0yU9+8u1tY2Njuuaaa/TAAw9oYGBAK1eu1GWXXably5cftW/WdnnhchYAJHj11Vd1wQUXaMWK\nFTrppJN0/vnn66KLLtJbb73V9rEuvvhinXbaaUdt27Vrl5YuXaqzzjpLc+fO1bp163TXXXcds2/W\ndnmhJ1KAqvUYuOQFdO6UU07R7t27tWvXLn35y18+5oP7Qx/6kH7zm98cs99Xv/pVrVq1atrj79u3\nT4sXL3779cDAgB555JHkdnmpfYjYHpY0vHTp0rJLqaTUsReCBchmdHRUK1asOGb7j370oxKqKV7t\nQ2SmTPEFUE/PPPOMPvCBDxyzvdOeSH9/v/bu3fv260ajof7+/uR2eal9iABAmfbv36+1a9ces73T\nnsjKlSv1/PPP68UXX1R/f7927Nih73znO8nt8kKIAKi9Mi+/rl69Whs3btT27dv14Q9/OOkY69ev\n10MPPaQDBw5oYGBAX/rSl7Rx40bdcMMNWr16tcbGxrRhw4ajLputXbtW3/rWt3TmmWcet13eHBGF\nfbM8DQ0NRVUej1u1gfRuYIwEVfLss8/qnHPOKbuM2pns92b78YgYSj0mU3wBAMkIEQBAMsZEkAn3\nlqBqIoLl4NuQ19AFIdIFM3EMBKiyvr4+HTx4kAdTZXTkyYZ9fX1dPzYh0iYCo4kbFFGmgYEBNRoN\nvfLKK2WXUhtHnrHebbUPEe5YB3rPnDlzuv6scKSp/cB6RIxExOb58+eXXQoA9JzahwgAoDyECAAg\nWe3HRFBdDL4DMx8hMg1mYwHA1AgRdA2BC/QexkQAAMkIEQBAMi5noVCswQXMLPREAADJCBEAQDJC\nBACQrPZjIt1cgJEpqgDQntqHSESMSBoZGhra1O6+hAYAdKb2IYJ6Y7YWUG+MiQAAkhEiAIBkhAgA\nIBkhAgBIRogAAJIxOwuVwoOsgHohRFA7BA1QHVzOAgAkoyeCymNlAaC6eipE+DACgO7ichYAIBkh\nAgBI1lOXszBzsZAjUI4ZGyKMfwBA/mZsiKC30TMBilHZELH9UUmXSjpF0s0R8YOSSwIATFDowLrt\nW2y/bHt0wvY1tp+zvcf2VkmKiDsjYpOkqyV9vMg6AQDZFD07a7ukNeM32J4t6UZJl0haLmm97eXj\nmnyx9T4AoGIKDZGIeFjSryZsvlDSnoh4ISJel7RD0uVuuk7SvRHxRJF1AgCyqcKYSL+kveNeNyR9\nUNLnJK2SNN/20oi4aeKOtjdL2ixJ73nPewooFXXFoo1APqoQIpOKiOslXT9Nm22StknS0NBQFFEX\nAOAdVbhjfZ+kxeNeD7S2AQAqrgoh8qikZbaX2J4raZ2ku7PubHvY9rZDhw7lViAAYHKFXs6yfZuk\nP5a00HZD0t9FxM22t0i6X9JsSbdExNNZjxkRI5JGTly0bBN3qQNAsQoNkYhYP8X2nZJ2FlkLAKBz\nlR1YB/LG0ihA52ofIraHJQ2fcOqisksBgJ5ThYH1jkTESERsntU3r+xSAKDn1D5EAADlqf3lLKBb\nuKsdaF/tQ4QxEQAoT+0vZzEmAgDlqX2IAADKU/vLWUCRuLcEOBohAhwHS+kAx1f7EGFgHQDKU/sx\nEQbWAaA8tQ8RAEB5CBEAQDJCBACQjBABACRjdhYAIFnteyLMzgKA8tS+JwKUiZV/0etq3xMBAJSH\nEAEAJCNEAADJCBEAQLLaD6wzxRcAyuOIKLuGrjhx0bJY9KmvlV0GkAkzuFAVth+PiKHU/bmcBQBI\nRogAAJIRIgCAZIQIACBZ7WdnAXU0cbkUBtpRV233RGzPsz07j2IAAPUybYjYnmX7E7a/b/tlSf8l\n6X9tP2P7H2wvzb9MAEAVZemJPCjp9yV9QdK7I2JxRJwh6Y8k/UTSdbb/Mscaj8v2sO1tbx1+rawS\nAKBnZRkTWRURb0zcGBG/knSHpDtsz+l6ZRlFxIikkRMXLdtUVg0A0KumDZEjAWL7ZEnnSvpFRPzP\nZG0ApOG5JKirdmZnPSjpp5LeZ/sVSRtavREABWBGF6ooc4iMX1vF9sck7bR9RUTsz6UyAEDltTXF\ntzVTq1/SS5LuU7N3AgDoUZl7IrZ/Ien3JP1S0v7Wv5Gc6gIA1EA7YyLvjYhXc6sEAFA7mZ8nYvtP\nJV0l6deSRiU9KWk0In6XX3nZ8TwRgMF2tK/T54m00xO5RdLnJc1Rc6rvRyWtkMQd6wDQo9oJkZci\n4s7W1/+eRzEAgHrJsnbWP9v+vKSf2P6bAmoCANRElim+2yVZ0rsk/ZXtl2zfbfvvbf9FrtUBACot\ny7InP5T0wyOvbZ8g6RxJ50m6UCVf2rI9LGn4hFMXlVkGAPSkLJezPP51RLwZEU9FxL9GxN9O1qZI\nETESEZtn9c0rqwQA6FlZBtYftH2HpLsi4udHNtqeq+Zy8J9S88717blUCCAz1tdC0bKEyBpJGyTd\nZvssSf+n5p3rsyT9QNLXImJ3fiUCAKoqS4i8KyK+LunrreeGLJT024j4db6lAQCqLkuI3Gv7DDUf\ni/ukpKckPWn7KZZBAYDelmV21nLbJ0paLukP1Lxb/XJJ59r+XUQsyblGAEBFZbpjvbU+1m7beyT9\nVs1LWu9Vs2cCAOhR04aI7bMlXSrpzyWdLukBSbdK2hwRr+dbHgCgyrL0RJ6VtFvSdWpO863Eqr0A\n0jANGN2UJUT+WtL7JV0j6QbbB9UcXH9K0lPjFmUEAPSYLAPr/zT+te0BvTPAfqUkQgQAelQ7S8FL\nkiKiIakh6d7ulwMAqJMsq/gCADCptnsiAOpj4iA60G30RAAAyeiJADjKZL0XpgFjKvREAADJCBEA\nQLJKhojts2zfbPv2smsBAEytsBCxfYvtl22PTti+xvZztvfY3ipJEfFCRGwsqjYAQJoieyLb1XxK\n4ttsz5Z0o6RL1Fxqfr3t5QXWBADoQGGzsyLiYduDEzZfKGlPRLwgSbZ3qPmskmeyHNP2ZkmbJWn2\nKad3rVagl2S5l4RFGzGVssdE+iXtHfe6Ianf9gLbN0m6wPYXpto5IrZFxFBEDM0+aX7etQIAJqjk\nfSIRcVDS1WXXAQA4vrJDZJ+kxeNeD7S2ZWZ7WNLwCacu6mZdAIAMyr6c9aikZbaX2J4raZ2ku9s5\nQESMRMTmWX3zcikQADC1Iqf43ibpx5LOtt2wvTEi3pS0RdL9aj5B8bsR8XRRNQEAOlPk7Kz1U2zf\nKWlnUXUAALqn7DGRjjEmAhQvyyKNWaYFM3W4/soeE+kYYyIAUJ7ahwgAoDyECAAgGWMiAAqR+qhe\nxk2qrfY9EcZEAKA8tQ8RAEB5CBEAQDJCBACQjIF1AF2ROnCOeqt9T4SBdQAoT+1DBABQHkIEAJCM\nEAEAJCNEAADJmJ0FoNayLEtf5HF6Te17IszOAoDy1D5EAADlIUQAAMkIEQBAMkIEAJCMEAEAJGOK\nLwAgWe17IkzxBYDy1D5EAADlIUQAAMkIEQBAMkIEAJCMEAEAJCNEAADJCBEAQDJuNgRQGZM90yOv\n4+b1rJCJ32umP5Ok9j0RbjYEgPLUPkQAAOUhRAAAyQgRAEAyQgQAkIwQAQAkI0QAAMkIEQBAMkIE\nAJCMEAEAJCNEAADJCBEAQDIWYARQK1kWacxrIccUWRZ/zLJoY1UXdqx9T4QFGAGgPLUPEQBAeQgR\nAEAyQgQAkIwQAQAkI0QAAMkIEQBAMkIEAJCMEAEAJCNEAADJCBEAQDJCBACQjBABACQjRAAAyQgR\nAEAyQgQAkIwQAQAkI0QAAMkq+Xhc2/MkfV3S65IeiohbSy4JADCJwnoitm+x/bLt0Qnb19h+zvYe\n21tbm6+QdHtEbJJ0WVE1AgDaU+TlrO2S1ozfYHu2pBslXSJpuaT1tpdLGpC0t9VsrMAaAQBtKOxy\nVkQ8bHtwwuYLJe2JiBckyfYOSZdLaqgZJP+p4wSd7c2SNkvS7FNO737RAGaswa3fb7vNz75yadvH\nSJHXcfNQ9sB6v97pcUjN8OiX9D1JV9r+hqSRqXaOiG0RMRQRQ7NPmp9vpQCAY1RyYD0iXpP0mbLr\nAAAcX9k9kX2SFo97PdDaBgCogbJD5FFJy2wvsT1X0jpJd7dzANvDtre9dfi1XAoEAEytyCm+t0n6\nsaSzbTdsb4yINyVtkXS/pGclfTcinm7nuBExEhGbZ/XN637RAIDjKnJ21voptu+UtLOoOgAA3VPJ\ngfV22B6WNHzCqYvKLgUAek7ZYyId43IWAJSn9iECACiPI6LsGrrC9iFJz7exy3xJh3Juc7z3Fko6\nMM2xqyTL76Iq36OT47Szb9a207VLfZ9zKN/vk3qcdvcr+7Po7Ig4eZpjTy0iZsQ/Sdu63b7TNtO8\n91jZv7M8f79lfo9OjtPOvlnbTtcu9X3OoWqeR732WTSTLmdNuTxKB+07bdNuTVVWxM/Sre/RyXHa\n2Tdr2+nadfp+XRT1c5R9HvXUZ9GMuZxVN7Yfi4ihsutAfXEOoRs6PY9mUk+kbraVXQBqj3MI3dDR\neURPBACQjJ4IACAZIQIASEaIAACSESIAgGSESAXYnmf727a/afuqsutBPdk+y/bNtm8vuxbUl+2P\ntj6L/s32R6ZrT4jkxPYttl+2PTph+xrbz9neY3tra/MVkm6PiE2SLiu8WFRWO+dRRLwQERvLqRRV\n1uZ5dGfrs+hqSR+f7tiESH62S1ozfoPt2ZJulHSJpOWS1tteruZjgfe2mo0VWCOqb7uyn0fAVLar\n/fPoi633j4sQyUlEPCzpVxM2XyhpT+v/GF+XtEPS5ZIaagaJxH8TjNPmeQRMqp3zyE3XSbo3Ip6Y\n7th8YBWrX+/0OKRmePRL+p6kK21/QzNnnSTkZ9LzyPYC2zdJusD2F8opDTUy1efR5yStkvQx21dP\nd5DaP9lwJoiI1yR9puw6UG8RcVDN69hAsoi4XtL1WdvTEynWPkmLx70eaG0D2sF5hG7oynlEiBTr\nUUnLbC+xPVfSOkl3l1wT6ofzCN3QlfOIEMmJ7dsk/VjS2bYbtjdGxJuStki6X9Kzkr4bEU+XWSeq\njfMI3ZDnecQqvgCAZPREAADJCBEAQDJCBACQjBABACQjRAAAyQgRAEAyQgQAkIwQAQAkI0SAjGz/\nme1/aaP9Tbb/MGVfoC4IESC78yTtbqP9RZJ+krgvUAuECJDdeZLebfth2z+3vcr2fNu/PNLA9uOt\nbedI+u+IGGt332J/JKAzhAiQ3XmSXomIiyVdK+mqiDgk6STbR57N81NJ56r5yNH7EvcFaoMQATKw\nPUfSAkn/2No0R9KvW1//QtKi1tfva71erVaItLuv7Xm2v237m7avyucnArqDEAGyOUfSTyPirdbr\ncyWNtr7eL+lM2x+TdEDNB/ucGhH72903Ip6XdIWk2yNik6TL8vyhgE7xeFwgm/PUvNx0xLmS7mp9\nvV/SWjUvYa2V9CeSHkzcV2o+Ye6p1tdjAiqMngiQzXmSnhz3+v06ujfxCUlXRsQBTT4eknVfSWqo\nGSQSf6OoOB5KBXSZ7SckfTAi3kjcf56kGyQdlvQfEXFrN+sDuokQAQAko6sMAEhGiAAAkhEiAIBk\nhAgAIBkhAgBIRogAAJIRIgCAZIQIACDZ/wN2wO+U7bJJ2wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1a186d5b50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fraction of transmitted photons: 0.102\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAERCAYAAAC6kZqPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEhZJREFUeJzt3X+MVXV+xvHnAdHpEsQEcRdnSIAFXNBK1oy2yVZpU1oQ\nO2rUZkHb1AUhNrpx078w2aQxjcma2qQ16FpWDbtbhVrdiCjgNlmNNv5AlKoj1EpElwvdFbCiscui\n46d/3Ds6wPy49zv3nHPP3PcrIZl77jlnPjMe7zOf8z3nexwRAgAgxbiiCwAAlBchAgBIRogAAJIR\nIgCAZIQIACAZIQIASEaIAACSESIAgGSECAAg2SlFFzBatnsk9UyaNGnV3Llziy4HAErllVdeORQR\nU1O391iZ9qS7uzt27NhRdBkAUCq2X4mI7tTtOZ0FAEhGiAAAkhEiAIBkpR9YB9B+Pv30U1UqFR09\nerToUkqjo6NDXV1dmjBhQlP3W/oQ6b86a/bs2UWXAiAnlUpFkyZN0owZM2S76HJaXkTo8OHDqlQq\nmjlzZlP3XfrTWRGxOSJWT548uehSAOTk6NGjmjJlCgFSJ9uaMmVKJp1b6UMEQHsiQBqT1e+LEAEA\nJCv9mEg7mbHmyRHXefcHl+VQCQBUESItrJ7QGGkbQgVAlgiRMY7uBe0g5Q+u4eT9/8SKFSv0xBNP\n6KyzzlJvb+8Xy7dt26ZbbrlFfX19uuGGG7RmzZpBt693vSyUfkzEdo/tdUeOHCm6lIbMWPPkcf8A\ntK/rr79e27ZtO25ZX1+fbrrpJm3dulW7du3Shg0btGvXrpO2rXe9rJS+E4mIzZI2d3d3ryq6ltEo\nMkgG+950J8DwPvroIy1cuFDHjh3T3r17NXfuXHV0dOj555/XuHGN/X1+ySWX6N133z1u2fbt2zV7\n9mzNmjVLkrRs2TJt2rRJ8+fPT1ovK6UPEQAowumnn66dO3dq+/btuv3227Vp06bj3r/44ov18ccf\nn7TdnXfeqUWLFo24//3792v69OlfvO7q6tJLL72UvF5WCBEMigF6oD69vb0699xzT1r+3HPPFVBN\n/giRHDDmAYxdu3bt0gUXXHDS8tF2Ip2dndq3b98XryuVijo7O5PXywoh0gTt8Fd7O/yMQIoDBw5o\n6dKlJy0fbSdy4YUX6u2339bevXvV2dmpjRs36qGHHkpeLyuECIDSK/KPmsWLF2vlypVav369Fi5c\nmLSP5cuX65lnntGhQ4fU1dWl2267TStXrtTatWu1ePFi9fX1acWKFcedNlu6dKnuu+8+nX322cOu\nlzUej9sE7Xi6ik4ERdq9e7fmzZtXdBmlM9jvjcfjAgAKU/rTWTxPpBjcWwJAGgOdCM8TAdrTWDkV\nn5esfl+lDxEA7aejo0OHDx8mSOrU/2TDjo6Opu+79KezALSfrq4uVSoVHTx4sOhSSqP/GevNRoig\nabiXBHmZMGFC058VjjSczgIAJCNEAADJCBEAQDJCBACQjIH1BrXjFCepuCERGPsIEeRqpBAmZIBy\n4XQWACAZnQhaCqfAgHIhRAbgAwwAGlP6EMl6Fl8G0gFgaKUPkYjYLGlzd3f3qka3JSAAYHQYWAcA\nJCNEAADJCBEAQDJCBACQrPQD641gIB0AmqutQgTlxMOugNZFiKB0uCkUaB2MiQAAktGJYEzglBdQ\nDDoRAEAyOhGMSXQmQD7oRAAAyQgRAEAyTmehLXBZMJCNMRsi3J0OANnjdBYAIBkhAgBI1rKns2xf\nKekySadLuj8ifl5wSQCAE+Taidh+wPb7tntPWL7E9lu299heI0kR8VhErJJ0o6Rv51knAKA+eXci\n6yWtlfST/gW2x0u6W9KfSKpIetn24xGxq7bK92vvA03FDYnA6OXaiUTEs5I+OGHxRZL2RMQ7EXFM\n0kZJV7jqDklbI+LVwfZne7XtHbZ3HDx4MNviAQAnaYWB9U5J+wa8rtSWfVfSIknX2L5xsA0jYl1E\ndEdE99SpU7OvFABwnJYdWI+IuyTdVXQdAIChtUInsl/S9AGvu2rL6mK7x/a6I0eONL0wAMDwWqET\neVnSHNszVQ2PZZKurXfjiNgsaXN3d/eqjOpDm2BqFKBxuYaI7Q2S/lDSmbYrkv42Iu63fbOkpySN\nl/RARLyZZ13AULiCCxheriESEcuHWL5F0pY8awEAjF4rjImMCmMiAFCc0odIRGyOiNWTJ08uuhQA\naDulDxEAQHEIEQBAMkIEAJCs9CHCwDoAFKcVbjYclf6bDU+bNmcVj8QFgHyVPkSAInGXO9odIQI0\ngG4XOF7px0QAAMUpfSdiu0dSzylnTCu6FABoO6XvRPrvWB/XMbHoUgCg7ZS+EwFaDTP/op2UvhMB\nABSHTgQoAN0KxorShwgD6xirCBqUQelPZzGwDgDFKX0nAowF3MSIsip9JwIAKA4hAgBIRogAAJIx\nJgJkjPEOjGV0IgCAZKXvRLhPBO2C+0bQikrfiXCfCAAUp/QhAgAoDiECAEhGiAAAkpV+YB3Alxh8\nR97oRAAAyQgRAEAyQgQAkKzhELE90fb4LIpJYbvH9rrPj35SdCkA0HZGDBHb42xfa/tJ2+9L+i9J\n/2N7l+2/tz07+zKHxs2GAFCcejqRpyV9XdKtkr4WEdMj4ixJfyDpRUl32P6LDGsEALSoei7xXRQR\nn564MCI+kPSopEdtT2h6ZQCAljdiJ9IfILYn2f6W7a8PtQ4AoL00crPh05Jek/QN2wclrah1IwAK\nwHNK0ArqDpGI6O7/2vY1krbYvioiDmRSGYBRGyxouIsdzdTQJb61K7U6Jb0naZuq3QkAoE3V3YnY\n/pWk35H0a0kHav82Z1QXAKAEGhkTmRsRH2VWCQCgdBoJkW7b10n6UFKvpNcl9UbEbzOpDEAhGEdB\nIxoJkQckfU/SBEnnS7pS0rmSCr1jHQBQnEZC5L2IeKz29b9lUQwAoFzqmTvrJ7a/J+lF23+TQ00N\nYQJGAChOPZf4rpdkSV+V9Je237P9uO2/s/3nmVZXByZgBIDijHg6KyJ+IekX/a9tnyJpnqQFki4S\np7aAUuERumimEUPEtiMi+l9HxGeS3qj9+5fB1gEAtId6Btaftv2opE0R8cv+hbZPVXU6+L9S9c71\n9ZlUCCBTzMGF0agnRJZIWiFpg+1Zkv5X1TvXx0n6uaR/jIid2ZUIAGhV9YTIVyPiHkn31J4bcqak\n30TEh9mWBgBodfWEyFbbZ6n6WNzXVR0Led32G0yDAgDtrZ6rs+bbPk3SfEm/q+rd6ldIOt/2byNi\nZsY1AgBaVF13rNfmx9ppe4+k36h6Smuuqp0JAKBN1XOJ7zmSLpP0Z5KmSvp3SQ9KWh0Rx7ItDwDQ\nyurpRHZL2inpDlUv82XWXgCApPpC5K8lnSfpJklrbR/WlzcbvjFgUkYAY1Qz7nJnivmxqZ6B9X8e\n+Np2l74cYL9aEiECAG2qkangJUkRUZFUkbS1+eUAAMqknll8AQAYFCECAEhGiAAAkhEiAIBkLRki\ntmfZvt/2I0XXAgAYWsNXZ6Wy/YCqd72/HxHnDVi+RNI/SRov6b6I+EFEvCNpJSEClEfKvSQ8ZbH8\n8uxE1qv6bJIv2B4v6W5Jl6o6weNy2/NzrAkAMAq5hUhEPCvpgxMWXyRpT0S8U5uHa6OqMwQDAEqg\n6DGRTkn7BryuSOq0PcX2vZK+afvWoTa2vdr2Dts7+v7vSNa1AgBOkNuYSCMi4rCkG+tYb52kdZJ0\n2rQ5kXVdAIDjFd2J7Jc0fcDrrtoyAEAJFN2JvCxpju2ZqobHMknXNrID2z2Sek45Y1oG5QHIU6vN\n9MvVYyPLrROxvUHSC5LOsV2xvTIiPpN0s6SnVH1uycMR8WYj+42IzRGxelzHxOYXDQAYVm6dSEQs\nH2L5Fklb8qoDANA8RY+JAABKrOgxkVFjTATI32BjF2hPpe9EGBMBgOKUPkQAAMUhRAAAyRgTAZCJ\nIsdNuL8jP6XvRBgTAYDilD5EAADFIUQAAMkIEQBAMgbWAbS0kQbJ6xnAb7WJHceS0nciDKwDQHFK\nHyIAgOIQIgCAZIQIACAZIQIASMbVWQBKhWnoW0vpOxGuzgKA4pQ+RAAAxSFEAADJCBEAQDJCBACQ\njBABACTjEl8AaDFlmjCy9J0Il/gCQHFKHyIAgOIQIgCAZIQIACAZIQIASEaIAACSESIAgGSECAAg\nGSECAEjGHesA2lLKw61O3KZV7yLPU+k7Ee5YB4DilD5EAADFIUQAAMkIEQBAMkIEAJCMEAEAJCNE\nAADJCBEAQDJCBACQjBABACQjRAAAyQgRAEAyJmAEgCaqZ5LGlMkf69nHSBNCNuP7nqj0nQgTMAJA\ncUofIgCA4hAiAIBkhAgAIBkhAgBIRogAAJIRIgCAZIQIACAZIQIASEaIAACSESIAgGSECAAgGSEC\nAEhGiAAAkhEiAIBkhAgAIBkhAgBIRogAAJK15ONxbU+UdI+kY5KeiYgHCy4JADCI3DoR2w/Yft92\n7wnLl9h+y/Ye22tqi6+S9EhErJJ0eV41AgAak+fprPWSlgxcYHu8pLslXSppvqTltudL6pK0r7Za\nX441AgAakNvprIh41vaMExZfJGlPRLwjSbY3SrpCUkXVIPlPDRN0tldLWi1J40+f2vyiAWAYM9Y8\n2ZR1spLH9y56YL1TX3YcUjU8OiX9TNLVtn8oafNQG0fEuojojoju8V+ZnG2lAICTtOTAekR8Iuk7\nRdcBABhe0Z3IfknTB7zuqi0DAJRA0SHysqQ5tmfaPlXSMkmPN7ID2z22131+9JNMCgQADC3PS3w3\nSHpB0jm2K7ZXRsRnkm6W9JSk3ZIejog3G9lvRGyOiNXjOiY2v2gAwLDyvDpr+RDLt0jaklcdAIDm\nKfp0FgCgxFry6qxG2O6R1HPKGdOKLgUA2k7pOxHGRACgOI6IomtoCttHJL3dwCaTJR3JeJ3h3jtT\n0qER9t1K6vldtMr3GM1+Gtm23nVHWi/1fY6hbL9P6n4a3a7oz6JzImLSCPseWkSMiX+S1jV7/dGu\nM8J7O4r+nWX5+y3ye4xmP41sW++6I62X+j7HUGseR+32WVT601kDDDk9yijWH+06jdbUyvL4WZr1\nPUazn0a2rXfdkdYb7ftlkdfPUfRx1FafRWPmdFbZ2N4REd1F14Hy4hhCM4z2OBpLnUjZrCu6AJQe\nxxCaYVTHEZ0IACAZnQgAIBkhAgBIRogAAJIRIgCAZIRIC7A90faPbf/I9nVF14Nysj3L9v22Hym6\nFpSX7Strn0X/avtPR1qfEMmI7Qdsv2+794TlS2y/ZXuP7TW1xVdJeiQiVkm6PPdi0bIaOY4i4p2I\nWFlMpWhlDR5Hj9U+i26U9O2R9k2IZGe9pCUDF9geL+luSZdKmi9pue35qj4WeF9ttb4ca0TrW6/6\njyNgKOvV+HH0/dr7wyJEMhIRz0r64ITFF0naU/uL8ZikjZKukFRRNUgk/ptggAaPI2BQjRxHrrpD\n0taIeHWkffOBla9OfdlxSNXw6JT0M0lX2/6hxs48ScjOoMeR7Sm275X0Tdu3FlMaSmSoz6PvSlok\n6RrbN460k9I/lGosiIhPJH2n6DpQbhFxWNXz2ECyiLhL0l31rk8nkq/9kqYPeN1VWwY0guMIzdCU\n44gQydfLkubYnmn7VEnLJD1ecE0oH44jNENTjiNCJCO2N0h6QdI5tiu2V0bEZ5JulvSUpN2SHo6I\nN4usE62N4wjNkOVxxCy+AIBkdCIAgGSECAAgGSECAEhGiAAAkhEiAIBkhAgAIBkhAgBIRogAAJIR\nIkCdbP+x7Z82sP69tr+Vsi1QFoQIUL8FknY2sP7vS3oxcVugFAgRoH4LJH3N9rO2f2l7ke3Jtn/d\nv4LtV2rL5kn674joa3TbfH8kYHQIEaB+CyQdjIhLJN0i6bqIOCLpK7b7n83zmqTzVX3k6LbEbYHS\nIESAOtieIGmKpH+oLZog6cPa17+SNK329TdqrxerFiKNbmt7ou0f2/6R7euy+YmA5iBEgPrMk/Ra\nRHxee32+pN7a1wcknW37GkmHVH2wzxkRcaDRbSPibUlXSXokIlZJujzLHwoYLR6PC9Rngaqnm/qd\nL2lT7esDkpaqegprqaQ/kvR04rZS9Qlzb9S+7hPQwuhEgPoskPT6gNfn6fhu4lpJV0fEIQ0+HlLv\ntpJUUTVIJP4fRYvjoVRAk9l+VdLvRcSnidtPlLRW0lFJ/xERDzazPqCZCBEAQDJaZQBAMkIEAJCM\nEAEAJCNEAADJCBEAQDJCBACQjBABACQjRAAAyf4f7CcbvkxBT3UAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10b0bd290>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fraction of transmitted photons: 0.101\n"
     ]
    }
   ],
   "source": [
    "# Run four Monte Carlo simulations:\n",
    "\n",
    "n_phot=100000 # number of photons\n",
    "H=1e7 # corona height\n",
    "v=0.1*c # electron velocity\n",
    "#v = np.sqrt(2.*k*T/m)\n",
    "Te = pi * me * v**2 /(8.*k) # Temperature of a Maxwellian distribution with mean velocity v \n",
    "h_nu_zero=1.6e-9 # 1 keV seed energy\n",
    "taus=np.array([0.1,1.,3.,10.]) # optical depths\n",
    "for tau in taus:\n",
    "    plot_mc(tau,H,v,h_nu_zero)\n",
    "    plot_mc_maxwellian(tau,H,Te,h_nu_zero)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}