{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 12,
   "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",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def lorentz_matrix(v):\n",
    "    #lorentz_matrix returns a 4x4 matrix of the lorentz transformation moving to a new reverence frame of velocity 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",
    "\n",
    "    ###Fill out code here\n",
    "    #matrix[0,0] = XXX\n",
    "\n",
    "    #Note: Instead of writing each element separately, one can write fewer lines of code with the following vector expressions\n",
    "\n",
    "\n",
    "    \n",
    "    return(matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# c = 3.e10 #cm/s \n",
    "# print(lorentz_matrix(np.asarray([0.1*c,0.2*c,0.3*c])))\n",
    "\n",
    "# [[ 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 ]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def lorentz_transform(p,v):\n",
    "    #lorentz_transform returns the four momentum vector in the new frame of reference\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",
    "\n",
    "    ###Fill out code here\n",
    "    \n",
    "    ###Return the appropriate quantity\n",
    "    return "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "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",
    "    # note: random orientation=even distribution in cos(theta)\n",
    "    ###Fill out code here\n",
    "\n",
    "\n",
    "    return(vector)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def scatter(p_in,v):\n",
    "    #p_in: The initial four momentum vector of the photon in the observer's frame\n",
    "    #v: the magnitude of the electron's velocity\n",
    "    \n",
    "    #choose a random direction for the velocity of the electron\n",
    "    \n",
    "    # first, transform to electron frame\n",
    "    \n",
    "    # neglect recoil, so energy unchanged\n",
    "    \n",
    "    # assume isotropic scattering, draw random new direction for the photon\n",
    "    \n",
    "    # new four momentum of photon in electron frame\n",
    "    \n",
    "    # transform back to observer frame\n",
    "    \n",
    "    #return the final four momentum vector after the scattering  in the observer's frame\n",
    "    return(p_out)\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "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",
    "    \n",
    "    ###Fill out code here\n",
    "    ### Take the following steps\n",
    "    # draw random probability\n",
    "\n",
    "    \n",
    "    # and find corresponding scattering depth\n",
    "\n",
    "    \n",
    "    # figure out how far and where the photon went\n",
    "\n",
    "    \n",
    "    # determine if it escaped or got absorbed\n",
    "\n",
    "    \n",
    "    # if it neither, execute inverse Compton scatter\n",
    "\n",
    "    \n",
    "    \n",
    "    #return the final four momentum vector\n",
    "    \n",
    "    return(p_out)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "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",
    "\n",
    "        # Define the initial photon four-momentum\n",
    "\n",
    "        \n",
    "        # define the initial position of photons\n",
    "\n",
    "        \n",
    "        # keep scattering until absorbed or escaped\n",
    "\n",
    "    \n",
    "    #return an array with the energies of all the photons that escaped\n",
    "    return"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "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": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "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",
    "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)"
   ]
  },
  {
   "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
}