{ "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 }