%pylab inline
from scipy import interpolate
from matplotlib.colors import LogNorm
from datetime import datetime, timedelta
from pyglow import pyglow
from scipy.stats import rv_discrete
import sys
import time
from IPython.display import display, clear_output

# Change matplotlib defaults
matplotlib.rcParams['savefig.dpi'] = 150
matplotlib.rcParams['figure.figsize'] = (4,3)
matplotlib.rcParams['font.size'] = 8
matplotlib.rcParams['savefig.bbox'] = 'tight'
matplotlib.rcParams['ytick.labelsize'] = 'small'
matplotlib.rcParams['xtick.labelsize'] = 'small'
matplotlib.rcParams['legend.fancybox'] = 'True'
matplotlib.rcParams['legend.fontsize'] = 'small'
matplotlib.rcParams['legend.numpoints'] = 1

Populating the interactive namespace from numpy and matplotlib

WARNING: pylab import has clobbered these variables: ['gamma', 'datetime', 'e']
`%matplotlib` prevents importing * from pylab and numpy

Oxygen aurora model v0.0

What's included:

• Initial preciptating neutral atom source at 1000 km
• Elastic collisions of fast $O$ on atmospheric $O$
  • differential scattering cross section from Ishimoto et al. (1992).
• Inelastic collisions:
  • Excitation to $O(^1D)$
  • excitation cross section from Ishimoto et al. (1994).
• MSIS neutral density (assumed to be entirely $O$)

What's not included:

• Charge exchange of precipitating $O^+$ to create fast $O$
  • Assumes $\sigma_{\textrm{chg ex}} >> \sigma_{\textrm{other things}}$ (Kozyra et al. 1982)
• Ionization collisions
• Stripping collisions
• Collisions with ambient $N_2$, $O_2$, $He$
• Thermal velocities of ambient particles (assumed to be at rest)
• Loss of $O(^1D)$
• Long lifetime of $O(^1D)$
• Excitation of other particles and states (e.g., $O(^1S)$ or $N_2$-band states)
• Gravity
• Magnetic field
• Momentum transfer in inelastic collisions (forward-scattering assumed)

O-O elastic collision cross section

from scipy.special import gamma

Semi-classical Landau-Schiff formula, given in Kharchenko et al (2000) and used for O+ precipitation modeling in Shematovich et al (2005)

mO_const = 15.9994 * 1.660538921e-27 # oxygen mass in kg
e_const = 1.602e-19 # elementary charge

def sigma_LS(energy_ev):
    '''
    Return Landau-Shiff cross section at energy given in eV
    '''
    n = 6.
    Cn = 16.72 # erg cm^6 1e60
    Cn_SI = 1e-60*1e-7*1e-12*Cn # J m^6
    hbar = 1.0545718e-34
    energy_J = energy_ev*e_const
    #v = sqrt(2*energy_J/mO_const)
    #print v
    #nu = energy_ev/6.582119e-16
    #print nu
    c0 = 2*np.pi**(n/(n-1)) * np.sin(np.pi/2*(n-3)/(n-1)) * gamma((n-3)/(n-1)) * \
            (gamma((n-1)/2)/gamma(n/2))**(2/(n-1))
    #x = sqrt(8e77*energy_ev) # Why.
    v = sqrt(2*energy_J/mO_const)
    sigma = c0 * (Cn_SI/(v*sqrt(np.pi)*hbar))**(2/(n-1))
    return sigma*1e4

e = linspace(0.1,5,100)
factor = 1
sig = sigma_LS(e)
print '%e should be close to 7.5e-15'% sigma_LS(1)


plot(e,sig*1e16,'k')
xlabel('Energy [eV]')
ylabel('Cross section [$10^{-16}$ cm$^2$]')
#ylim(50,120)

7.446066e-15 should be close to 7.5e-15

<matplotlib.text.Text at 0x7fdb16eb5590>

This doesn't really make a whole lot of sense for high energies, which is where it matters, I think.

The DSCS used in Ishimoto 1992

def p_i92(tau):
    if tau < 1:
        p = 6e-17
    elif tau < 2e2:
        p = exp(log(6e-17) + (log(1e-18)-log(6e-17))/(log(2e2)-log(1)) * (log(tau)-log(1)))
    else:
        p = exp(log(1e-21) + (log(1e-18)-log(1e-21))/(log(2e2)-log(5e3)) * (log(tau)-log(5e3)))
    return p

def dscs(e,th):
    '''
    Return the DSCS*sin(th), which should be interpreted as a pdf of the scattering angle,
    multiplied by the total cross section.
    INPUTS:
        e - incident energy, eV
        th - scattering angle, radians
    OUTPUT:
        s - see above. cm^2/rad
    '''
    if th==0.0: # Hack to avoid singularity at th=0
        th = 1e-16
    
    tau = e*th*180./pi
    p = p_i92(tau)
    s = p/th
    return s
    

tau = logspace(-2, 5, 100)
p = zeros(len(tau))
for i in range(len(tau)):
    p[i] = p_i92(tau[i])

loglog(tau,p,'k-')
grid()
xlabel('tau')
ylabel('p')
title('Ishimito 92 O-O DSCS')

<matplotlib.text.Text at 0x7f053d2b2ed0>

Pre-calculate total elastic cross section as function of energy

e_tr_min = 1e-9 # eV. energy transfer below which collision is ignored.
#e_min = 2*e_tr_min # eV. energy below which particle is no longer tracked
e_min = 0.1
e_vec = logspace(log10(e_min),7,101)
N = 10000 # number of angular bins for integration
elastic_xsec_vec = zeros(len(e_vec))
cdf_mat = zeros((len(e_vec),N))


th_edge = logspace(-9,log10(pi),N)
th_cdf = th_edge.copy() # better name, to indicate it's the grid for the CDF
th_mid  = (th_edge[1:] + th_edge[:-1])/2
for i in range(len(e_vec)):
    e = e_vec[i]
    if e >= e_tr_min: # Use minimum scattering angle
        th_min = 2*arcsin(sqrt(e_tr_min/e))
    else: # do full calculation
        th_min = th_edge[0]
    jstart = sum(th_edge < th_min)
    
    s_cdf = zeros(len(th_edge))
    for j in range(jstart,len(th_mid)):
        s = dscs(e,th_mid[j])
        dx = th_edge[j+1] - th_edge[j]
        s_cdf[j+1] = s_cdf[j] + s*dx
    scatt_angle_cdf = s_cdf / s_cdf[-1]
    
    elastic_xsec_vec[i] = s_cdf[-1]
    cdf_mat[i,:] = scatt_angle_cdf

Total elastic scattering cross section

loglog(e_vec,elastic_xsec_vec,'k-')

xlabel('Energy [eV]');
ylabel('Elastic Scattering Cross section [cm$^2$]');

X,Y = meshgrid(th_cdf, e_vec)
pcolormesh(X,Y,cdf_mat,cmap='jet')
xscale('log')
yscale('log')
colorbar()
title('CDF of scattering angle')
xlabel('Scattering angle [rad]')
ylabel('Incident energy [eV]')

<matplotlib.text.Text at 0x7f04d3fc5bd0>

Differential scattering cross section @ 100 eV

ei = 37 # eV
#print e_vec[ei]
s_pdf = diff(cdf_mat[ei,:])
loglog(th_cdf[1:],s_pdf,'k-')
xlabel('Scattering angle [rad]')
ylabel('PDF')
xlim((1e-5,pi));

Create interpolation variables to be used in model

elastic_xsec = interpolate.InterpolatedUnivariateSpline(e_vec,elastic_xsec_vec, ext=3)

For each energy, calculate inverse CDF and reconfigure into a matrix for 2D interpolation

u = linspace(0,1,101)
icdf_mat = zeros((len(e_vec),len(u)))
for i in range(len(e_vec)):
    cdf_i = cdf_mat[i,:]
    
    idx = cdf_i > 0.0
    cdf_i_new = concatenate(([0.0], cdf_i[idx]))
    th_cdf_new = concatenate(([0.0], th_cdf[idx]))
    
    th_interp = interpolate.interp1d(cdf_i_new,th_cdf_new)
    thu = th_interp(u)

    icdf_mat[i,:] = thu

X,Y = meshgrid(u, e_vec)
pcolormesh(X,Y,icdf_mat,cmap='jet',norm=LogNorm())
yscale('log')
colorbar()
title('Inverse CDF of scattering angle')
xlabel('Value of CDF')
ylabel('Incident energy [eV]')

<matplotlib.text.Text at 0x7f053a24b310>

elastic_scatt_angle_icdf = interpolate.RectBivariateSpline(e_vec,u,icdf_mat)

O-O excitation cross section to O(1D)

Grean-McNeal formula (see Ishimoto 1994)

def excitation_xsec(e):
    '''
    The excitation cross section of the O + O --> O + O(1D) reaction,
    as a function of energy, e
    INPUT:
        e - incident energy, eV
    OUTPUT:
        sigma - cross section, cm^2
    '''
    sig_max = 3.8e-18
    nu = 1.0
    om = 1.0
    J = 5.4e3   

    pre = J**nu/(2*J**(om+nu))
    sigma =  sig_max/pre * (e**nu/(J**(om+nu) + e**(om+nu)))
    
    # But if e < 1.96 eV, it can't excite
    if e < 1.96:
        sigma = 0.0
    
    return sigma

Excitation cross section

e = logspace(0, 6, 1000)

sigma = zeros(len(e))
for i in range(len(e)):
    sigma[i] = excitation_xsec(e[i])

loglog(e,sigma,'k-')


xlabel('Energy [eV]');
ylabel('Excitation Cross section [cm$^2$]');

Neutral density

lat = 45
lon = 0
t = datetime(2004,1,1,0) # Solar medium

altvec = linspace(30,1200,1000)
n_O = zeros(len(altvec))
n_N2 = zeros(len(altvec))
n_He = zeros(len(altvec))
n_tot = zeros(len(altvec))
for i in range(len(altvec)):
    pt = pyglow.Point(t,lat,lon,altvec[i])
    pt.run_msis()
    n_O[i] = pt.nn['O']
    n_N2[i] = pt.nn['N2']
    n_He[i] = pt.nn['HE']
    n_tot[i] = sum([pt.nn[x] for x in ['AR','H','HE','N','N2','O','O2']])

/usr/local/lib/python2.7/dist-packages/scipy-0.15.1-py2.7-linux-x86_64.egg/scipy/stats/stats.py:310: RuntimeWarning: invalid value encountered in double_scalars
  return np.mean(x, axis) / factor

semilogx(n_tot,altvec,'k-',linewidth=2,label='total')
semilogx(n_O, altvec,'b-',label='O')
semilogx(n_N2, altvec,'r-', label='N2')
semilogx(n_He, altvec,'g-', label='He')
legend(loc='best')
xlim((1e5,1e12))
xlabel('Neutral density [cm^-3]')
ylabel('Altitude [km]')

<matplotlib.text.Text at 0x7f053804f710>

Use total density, but assume that it's all atomic oxygen

density_O = interpolate.InterpolatedUnivariateSpline(altvec*1e3, n_tot)

altvec = linspace(60e3,1000e3,100000)
nO_interp = zeros(len(altvec))
for i in range(len(altvec)):
    nO_interp[i] = density_O(altvec[i])
    
semilogx(nO_interp, altvec, 'k-')

[<matplotlib.lines.Line2D at 0x7f0538189a90>]

Model, version 0

Following notes in notebook on 4 Nov 2015.

Variables/functions used from above:

• density_O(altitude_m)
• excitation_xsec(energy_ev)
• elastic_xsec(energy_ev)
• elastic_scatt_angle_icdf(energy_ev, u)

# Constants
m_O = 16*1.67e-27 # kg
q = 1.6e-19 # fundamental charge
k = 8.617e-5 # ev/K
#e_thermal = 1.5*k*200 # thermal energy at 200 K
#e_min = 5*e_thermal # stop when particle has less than this energy. Convert to heat.
e_min = 1.96 # 

# Initialize variables for this particle
e = 1e3 # eV
ru = rand() # randomized initial direction as per:
rv = rand() # http://mathworld.wolfram.com/SpherePointPicking.html
th = arccos(rv-1) # isotropic over downward hemisphere
ph = 2*pi*ru
v_mag = sqrt(2*e*q/m_O)
v = v_mag * array([sin(th)*sin(ph), sin(th)*cos(ph), cos(th)])
r = array([0, 0, 1000e3]) # m/s [x,y,z]
state = 0 # 0 = 3P, 1 = 1D

# Variables to record secondary particles 
secondary_count = 0
secondary_v = []
secondary_r = []
secondary_state = []

# For local, thermalized primaries and secondaries, record their energy and altitude
final_e = [] # eV
final_alt = [] # m

# Temporary variables for tracking particles
r_history = []
r_elastic_collision = []
r_excitation_collision = []
t_history = []

steps = 0
t = 0
steps_max = 1e5
stop = False

while(not stop):

    # Density and cross sections
    e = 0.5*m_O/q*norm(v)**2 # eV
    sigma_elastic    = elastic_xsec(e)    # cm^2
    sigma_excitation = excitation_xsec(e) # cm^2
    n_O = density_O(r[2]) # cm^-3

    # Mean free path
    mfp = 1e-2/(n_O*(sigma_elastic + sigma_excitation)) # m

    # Step size
    H = 50e3 # scale size, approximate for now
    ds = 0.1 * min(H,mfp) # so that step is shorter than the scale size and collisional scale
    dt = ds/norm(v) # sec

    # Determine if collision happened
    pc = 1 - exp(-ds/mfp)
    collision = rand() < pc

    if collision:
        # Determine type of collision
        collision_probs = array([sigma_elastic, sigma_excitation])
        collision_pmf = collision_probs/sum(collision_probs)
        # Randomly sample from the pmf
        collision_cdf = cumsum(collision_pmf)
        u = rand()
        collision_type = sum(u>=collision_cdf)
        ################## Elastic collision ########################
        if collision_type == 0:
            # Random scattering angle
            chi = elastic_scatt_angle_icdf.ev(e, rand())
            # Rotate to collision frame
            v1 = v.copy()
            v2 = 0.0 # assume ambient particle is at rest
            ez = v1/norm(v1)
            ey_tmp = cross(v1, array([0,1,0]))
            ey = ey_tmp/norm(ey_tmp)
            ex = cross(ey,ez)
            A = zeros((3,3))
            A[:,0] = ex
            A[:,1] = ey
            A[:,2] = ez
            v1p = A.T.dot(v1-v2)
            v2p = 0.0
            # Velocity magnitudes after collision
            m1 = m_O
            m2 = m_O
            v1p_mag_new = sqrt(m1**2 + 2*m1*m2*cos(chi) + m2**2)/(m1+m2) * norm(v1p)
            v2p_mag_new = 2*m1/(m1+m2) * norm(v1p) * sin(chi/2)
            # Scattering angles
            th1p = arctan2(m2*sin(chi), m1 + m2*cos(chi))
            th2p = (pi-chi)/2
            phi1p = 2*pi*rand()
            # New velocity vectors
            v1p_new = v1p_mag_new * array([sin(th1p)*sin(phi1p),    sin(th1p)*cos(phi1p),    cos(th1p)])
            v2p_new = v2p_mag_new * array([sin(th2p)*sin(phi1p+pi), sin(th2p)*cos(phi1p+pi), cos(th2p)])
            v1_new = A.dot(v1p_new) + v2
            v2_new = A.dot(v2p_new) + v2
            # Record secondary particle for later
            e2 = 0.5*m_O/q*norm(v2_new)**2 
            if e2 < e_min: # No need to trace secondary. Just record its energy.
                final_e.append(e2)
                final_alt.append(r[2])
            else:
                secondary_count += 1
                secondary_v.append(v2_new)
                secondary_r.append(r.copy())
                secondary_state.append(0)
            # Adjust this particle
            v = v1_new

            r_elastic_collision.append(r)

        ################# Excitation collision #######################
        elif collision_type == 1:
            # TODO: treat this as an elastic collision instead of a forward-scatter collision?
            # Determine which particle gets excited
            if rand() > 0.5: # Ambient particle
                e2 = 0.0 # It is left undisturbed (for now)
                if e2 < e_min: # Don't need to trace particle later.
                    # TODO: record O(1D) local production to be used later in emission calculation
                    pass
                else: # This can't occur, for now
                    pass

            else: # Incident particle
                state = 1 # TODO: record information for O(1D) production rate throughout
                # No need to record anything about secondary particle
            # Subtract excitation energy from incident particle (not ambient particle?)
            e -= 1.96 # eV
            v_mag_new = sqrt(2*e*q/m_O)
            v = v_mag_new * v / norm(v)
            
            r_excitation_collision.append(r)

    else:
        # Step forward in time
        if ds*v[2]/norm(v) > 20e3:
            print ds, v, e
            raise Exception('blah')
        r = r + ds*v/norm(v)
        t = t + dt
    r_history.append(r)
    t_history.append(t)

    # Check halting conditions
    steps += 1
    if r[2] > 1100e3: # Escape flux
        stop = True
        print 'Particle exited above 1100 km with energy %.1f eV' % e
    if e < e_min:
        stop = True
        print 'Particle thermalized at %.1f km with energy %.2f eV' % (r[2]/1e3, e)
    if r[2] < 20e3: # This probably shouldn't happen
        stop = True
        print 'Particle got to 20 km'
    if steps > steps_max:
        stop = True
        print 'Max steps reached (TEMPORARY)'
        
print '%i steps' % steps
print '%i elastic collisions' % len(r_elastic_collision)
print '%i excitation collisions' % len(r_excitation_collision)
print '%i secondaries to trace' % secondary_count

Particle thermalized at 111.6 km with energy 0.17 eV
11884 steps
949 elastic collisions
1 excitation collisions
11 secondaries to trace

rh = array(r_history)
plot(rh[:,0]/1e3,rh[:,2]/1e3,'k-',linewidth=0.5) # Particle track
plot(rh[0,0]/1e3,rh[0,2]/1e3,'ko') # Particle start
rhel = array(r_elastic_collision)
plot(rhel[:,0]/1e3,rhel[:,2]/1e3,'k.',markersize=2)
if len(r_excitation_collision) > 0:
    rhex = array(r_excitation_collision)
    plot(rhex[:,0]/1e3,rhex[:,2]/1e3,'r.')

xlabel('x [km]');
ylabel('z [km]');
xlim((min(rh[:,0]/1e3)-10, max(rh[:,0]/1e3)+10));
ylim((50,1100));
xlim((-3500,100));
title('Path of incident particle');

Run entire tree of particles

def trace_particle(r, v, state, depth, max_depth = np.inf):
    '''
    Recursive function to trace a single particle. It calls itself to trace its secondaries.
    '''
    
    # For local, thermalized primaries and secondaries, record their energy and altitude
    final_e = [] # eV
    final_alt = [] # m
    
    # For tracking O(1D) production
    excited_alt = []
    excited_v = []

    # Temporary variables for tracking this particle
    r_history = []
    r_elastic_collision = []
    r_excitation_collision = []
    t_history = []
    
    # For tracking secondary particle locations:
    r_histories = []
    secondary_particles = 0

    steps = 0
    t = 0
    steps_max = 1e5
    stop = False

    while(not stop):

        # Density and cross sections
        e = 0.5*m_O/q*norm(v)**2 # eV
        sigma_elastic    = elastic_xsec(e)    # cm^2
        sigma_excitation = excitation_xsec(e) # cm^2
        n_O = density_O(r[2]) # cm^-3

        # Mean free path
        mfp = 1e-2/(n_O*(sigma_elastic + sigma_excitation)) # m

        # Step size
        H = 50e3 # scale size, approximate for now
        ds = 0.1 * min(H,mfp) # so that step is shorter than the scale size and collisional scale
        dt = ds/norm(v) # sec

        # Determine if collision happened
        pc = 1 - exp(-ds/mfp)
        collision = rand() < pc

        if collision:
            # Determine type of collision
            collision_probs = array([sigma_elastic, sigma_excitation])
            collision_pmf = collision_probs/sum(collision_probs)
            # Randomly sample from the pmf
            collision_cdf = cumsum(collision_pmf)
            u = rand()
            collision_type = sum(u>=collision_cdf)
            ################## Elastic collision ########################
            if collision_type == 0:
                # Random scattering angle
                chi = elastic_scatt_angle_icdf.ev(e, rand())
                # Rotate to collision frame
                v1 = v.copy()
                v2 = 0.0 # assume ambient particle is at rest
                ez = v1/norm(v1)
                ey_tmp = cross(v1, array([0,1,0]))
                ey = ey_tmp/norm(ey_tmp)
                ex = cross(ey,ez)
                A = zeros((3,3))
                A[:,0] = ex
                A[:,1] = ey
                A[:,2] = ez
                v1p = A.T.dot(v1-v2)
                v2p = 0.0
                # Velocity magnitudes after collision
                m1 = m_O
                m2 = m_O
                v1p_mag_new = sqrt(m1**2 + 2*m1*m2*cos(chi) + m2**2)/(m1+m2) * norm(v1p)
                v2p_mag_new = 2*m1/(m1+m2) * norm(v1p) * sin(chi/2)
                # Scattering angles
                th1p = arctan2(m2*sin(chi), m1 + m2*cos(chi))
                th2p = (pi-chi)/2
                phi1p = 2*pi*rand()
                # New velocity vectors
                v1p_new = v1p_mag_new * array([sin(th1p)*sin(phi1p),    sin(th1p)*cos(phi1p),    cos(th1p)])
                v2p_new = v2p_mag_new * array([sin(th2p)*sin(phi1p+pi), sin(th2p)*cos(phi1p+pi), cos(th2p)])
                v1_new = A.dot(v1p_new) + v2
                v2_new = A.dot(v2p_new) + v2                # Record secondary particle for later
                e2 = 0.5*m_O/q*norm(v2_new)**2                
                if e2 < e_min: # No need to trace secondary. Just record its energy.
                    final_e.append(e2)
                    final_alt.append(r[2])
                elif depth >= max_depth: # Just record final quantities and move on
                    final_e.append(e2)
                    final_alt.append(r[2])
                else: # Trace all secondaries
                    fe,fa,ea,ev,rh = trace_particle(r, v2_new, 0, depth+1, max_depth)
                    r_histories.extend(rh)
                    final_e.extend(fe)
                    final_alt.extend(fa)
                    excited_alt.extend(ea)
                    excited_v.extend(ev)
                    secondary_particles += 1
                    
                # Adjust this particle
                v = v1_new

                r_elastic_collision.append(r)

            ################# Excitation collision #######################
            elif collision_type == 1:
                # TODO: treat this as an elastic collision instead of a forward-scatter collision?
                # Determine which particle gets excited
                if rand() > 0.5: # Ambient particle
                    e2 = 0.0 # It is left undisturbed (for now)
                    if e2 < e_min: # Don't need to trace particle later.
                        # TODO: record O(1D) local production to be used later in emission calculation
                        excited_alt.append(r[2])
                        excited_v.append(array([0,0,0]))
                    else: # This can't occur, for now
                        pass

                else: # Incident particle
                    # Record information for O(1D) production rate. No need to record secondary particle
                    state = 1 # TODO: record information for O(1D) production rate throughout
                    excited_alt.append(r[2])
                    excited_v.append(v)

                # Subtract excitation energy from incident particle (not ambient particle?)
                e -= 1.96 # eV
                v_mag_new = sqrt(2*e*q/m_O)
                v = v_mag_new * v / norm(v)

                r_excitation_collision.append(r)

        else:
            # Step forward in time
            r = r + ds*v/norm(v)
            t = t + dt
        r_history.append(r)
        t_history.append(t)

        # Check halting conditions
        steps += 1
        if r[2] > 1100e3: # Escape flux
            stop = True
            #print 'Particle exited above 1100 km with energy %.1f eV' % e
        if e < e_min:
            stop = True
            #print 'Particle thermalized at %.1f km with energy %.2f eV' % (r[2]/1e3, e)
        if r[2] < 20e3: # This probably shouldn't happen
            stop = True
            #print 'Particle got to 20 km'
        if steps > steps_max:
            stop = True
            #print 'Max steps reached (TEMPORARY)'
        
    # Concatenate this path with the secondaries
    r_histories.append(r_history)
    final_e.append(e)
    final_alt.append(r[2])
    
#    print '\t%i Depth level' % depth
#    print '\t%i steps' % steps
#    print '\t%i elastic collisions' % len(r_elastic_collision)
#    print '\t%i excitation collisions' % len(excited_alt)
#    print '\t%i secondaries to trace' % secondary_count
    
    return final_e, final_alt, excited_alt, excited_v, r_histories

# Use monoenergetic flux
e_0 = 20e3 # eV
N_particles = 100

final_e_all = []
final_alt_all = []
excited_alt_all = []
excited_v_all = []

t0 = datetime.now()
for i in range(N_particles):

    # Initialize variables for this particle
    e = e_0
    ru = rand() # randomized initial direction as per:
    rv = rand() # http://mathworld.wolfram.com/SpherePointPicking.html
    th = arccos(rv-1) # isotropic over downward hemisphere
    ph = 2*pi*ru
    v_mag = sqrt(2*e*q/m_O)
    v = v_mag * array([sin(th)*sin(ph), sin(th)*cos(ph), cos(th)])
    r = array([0, 0, 1000e3]) # m/s [x,y,z]
    state = 0 # 0 = 3P, 1 = 1D

    final_e, final_alt, excited_alt, excited_v, rh_all = trace_particle(r,v,state,0, max_depth=inf)
    
    final_e_all.extend(final_e)
    final_alt_all.extend(final_alt)
    excited_alt_all.extend(excited_alt)
    excited_v_all.extend(excited_v)
    
    clear_output(wait=True)
    time.sleep(0.1)
    print '%03i/%03i' % (i+1, N_particles)
    sys.stdout.flush()
t1 = datetime.now()

029/100

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-156-e19cce5ce64c> in <module>()
     22     state = 0 # 0 = 3P, 1 = 1D
     23 
---> 24     final_e, final_alt, excited_alt, excited_v, rh_all = trace_particle(r,v,state,0, max_depth=inf)
     25 
     26     final_e_all.extend(final_e)

<ipython-input-126-6855b0906604> in trace_particle(r, v, state, depth, max_depth)
     94                     final_alt.append(r[2])
     95                 else: # Trace all secondaries
---> 96                     fe,fa,ea,ev,rh = trace_particle(r, v2_new, 0, depth+1, max_depth)
     97                     r_histories.extend(rh)
     98                     final_e.extend(fe)

<ipython-input-126-6855b0906604> in trace_particle(r, v, state, depth, max_depth)
     94                     final_alt.append(r[2])
     95                 else: # Trace all secondaries
---> 96                     fe,fa,ea,ev,rh = trace_particle(r, v2_new, 0, depth+1, max_depth)
     97                     r_histories.extend(rh)
     98                     final_e.extend(fe)

<ipython-input-126-6855b0906604> in trace_particle(r, v, state, depth, max_depth)
     94                     final_alt.append(r[2])
     95                 else: # Trace all secondaries
---> 96                     fe,fa,ea,ev,rh = trace_particle(r, v2_new, 0, depth+1, max_depth)
     97                     r_histories.extend(rh)
     98                     final_e.extend(fe)

<ipython-input-126-6855b0906604> in trace_particle(r, v, state, depth, max_depth)
     94                     final_alt.append(r[2])
     95                 else: # Trace all secondaries
---> 96                     fe,fa,ea,ev,rh = trace_particle(r, v2_new, 0, depth+1, max_depth)
     97                     r_histories.extend(rh)
     98                     final_e.extend(fe)

<ipython-input-126-6855b0906604> in trace_particle(r, v, state, depth, max_depth)
     94                     final_alt.append(r[2])
     95                 else: # Trace all secondaries
---> 96                     fe,fa,ea,ev,rh = trace_particle(r, v2_new, 0, depth+1, max_depth)
     97                     r_histories.extend(rh)
     98                     final_e.extend(fe)

<ipython-input-126-6855b0906604> in trace_particle(r, v, state, depth, max_depth)
     94                     final_alt.append(r[2])
     95                 else: # Trace all secondaries
---> 96                     fe,fa,ea,ev,rh = trace_particle(r, v2_new, 0, depth+1, max_depth)
     97                     r_histories.extend(rh)
     98                     final_e.extend(fe)

<ipython-input-126-6855b0906604> in trace_particle(r, v, state, depth, max_depth)
     94                     final_alt.append(r[2])
     95                 else: # Trace all secondaries
---> 96                     fe,fa,ea,ev,rh = trace_particle(r, v2_new, 0, depth+1, max_depth)
     97                     r_histories.extend(rh)
     98                     final_e.extend(fe)

<ipython-input-126-6855b0906604> in trace_particle(r, v, state, depth, max_depth)
     94                     final_alt.append(r[2])
     95                 else: # Trace all secondaries
---> 96                     fe,fa,ea,ev,rh = trace_particle(r, v2_new, 0, depth+1, max_depth)
     97                     r_histories.extend(rh)
     98                     final_e.extend(fe)

<ipython-input-126-6855b0906604> in trace_particle(r, v, state, depth, max_depth)
     94                     final_alt.append(r[2])
     95                 else: # Trace all secondaries
---> 96                     fe,fa,ea,ev,rh = trace_particle(r, v2_new, 0, depth+1, max_depth)
     97                     r_histories.extend(rh)
     98                     final_e.extend(fe)

<ipython-input-126-6855b0906604> in trace_particle(r, v, state, depth, max_depth)
     30 
     31         # Density and cross sections
---> 32         e = 0.5*m_O/q*norm(v)**2 # eV
     33         sigma_elastic    = elastic_xsec(e)    # cm^2
     34         sigma_excitation = excitation_xsec(e) # cm^2

/usr/local/lib/python2.7/dist-packages/numpy-1.9.0-py2.7-linux-x86_64.egg/numpy/linalg/linalg.pyc in norm(x, ord, axis)
   2050 
   2051     """
-> 2052     x = asarray(x)
   2053 
   2054     # Check the default case first and handle it immediately.

/usr/local/lib/python2.7/dist-packages/numpy-1.9.0-py2.7-linux-x86_64.egg/numpy/core/numeric.pyc in asarray(a, dtype, order)
    460 
    461     """
--> 462     return array(a, dtype, copy=False, order=order)
    463 
    464 def asanyarray(a, dtype=None, order=None):

KeyboardInterrupt: 

t1 = datetime.now()

N_particles = i # However many were actually run

t = (t1-t0).total_seconds()
t_per = t/N_particles
print '%f sec per particle' % t_per
print '%f hours for N=1e2' % (t_per*1e2/3600)
print '%f days for N=1e6' % (t_per*1e6/86400)

1775.628535 sec per particle
49.323015 hours for N=1e2
20551.256193 days for N=1e6

final_e = array(final_e_all)
final_alt = array(final_alt_all)
excited_alt = array(excited_alt_all)
excited_v = array(excited_v_all)
rh_all = array(rh_all)

Path Visualizations

figure(figsize=(10,10))

rh = array(rh_all[-1])
plot(rh[:,0]/1e3,rh[:,2]/1e3,'k-',linewidth=2) # Particle track

for i in range(len(rh_all)-1):
    rh = array(rh_all[i])
    if len(rh) > 0:
        plot(rh[:,0]/1e3,rh[:,2]/1e3,'-',linewidth=0.5) # Particle track


xlabel('x [km]')
ylabel('z [km]')
#xlim((min(rh[:,0]/1e3)-10, max(rh[:,0]/1e3)+10))
#ylim((50,400))
#xlim((-2500,-1500))

<matplotlib.text.Text at 0x7f0497660a10>

figure(figsize=(10,10))

rh = array(rh_all[-1])
plot(rh[:,0]/1e3,rh[:,2]/1e3,'k-',linewidth=2) # Particle track
xfinal = rh[-1,0]

for i in range(len(rh_all)-1):
    rh = array(rh_all[i])
    if len(rh) > 0:
        plot(rh[:,0]/1e3,rh[:,2]/1e3,'-',linewidth=0.5) # Particle track

rh = array(rh_all[-1])
plot(rh[:,0]/1e3,rh[:,2]/1e3,'k-',linewidth=2) # Particle track

xlabel('x [km]');
ylabel('z [km]');
#xlim((min(rh[:,0]/1e3)-10, max(rh[:,0]/1e3)+10))
ylim((90,150));
xlim((xfinal/1e3-30,xfinal/1e3+30));

print '%i secondaries' % len(rh_all)

7647 secondaries

idx_esc = (final_alt >= 1100e3)
e_tot = sum(final_e) + 1.96*len(excited_alt)
e_esc = sum(final_e[idx_esc])
e_heat = sum(final_e[~idx_esc])
e_exc = 1.96*len(excited_alt)
print 'Escape energy: %f eV (%.2f %%)' % (e_esc, 100*e_esc/e_tot)
print 'Heat energy: %f eV (%.2f %%)' % (e_heat, 100*e_heat/e_tot)
print 'Excitation energy: %f eV (%.2f %%)' % (e_exc, 100*e_exc/e_tot)
print 'Total energy: %f eV' % (e_tot)

Escape energy: 198928.330998 eV (34.29 %)
Heat energy: 370947.318173 eV (63.95 %)
Excitation energy: 10207.680000 eV (1.76 %)
Total energy: 580083.329172 eV

Heating

alt_bin_edge = linspace(50e3,1100e3,100)
alt  = (alt_bin_edge[1:] + alt_bin_edge[:-1])/2
e_dep = zeros(len(alt))
for i in range(len(alt)):
    idx = (final_alt < alt_bin_edge[i+1]) & (final_alt >= alt_bin_edge[i])
    e_dep[i] = sum(final_e[idx])

# Scale to a 1 erg/cm^3/s source
e_model = e_0*N_particles # particle energy * number of particles
d_alt = alt[2]-alt[1]
heat_rate = e_dep/d_alt*6.2415e11/e_model*1e-2 # eV/cm^3/s. See notes.
    
heat_rate[heat_rate==0.0] = 1e-12 # because

semilogx(heat_rate,alt/1e3,'k.-')

xlabel('Heating rate [eV/cm^3/s]')
ylabel('Altitude [km]')

<matplotlib.text.Text at 0x7f04e3337510>

alt_i92 =       [450, 400, 300, 200 , 151, 141, 131, 125, 120,  110]
heat_rate_i92 = [20,  70 ,1000, 15e3, 50e3, 60e3, 50e3, 30e3, 10e3, 1e1 ]

semilogx(heat_rate,alt/1e3,'k.-', label='Me')
semilogx(heat_rate_i92, alt_i92, 'g-', linewidth=2, label='Ishimoto92')
legend(loc='best');
xlim((1e-3, 1e6));
ylim((0,600));
xlabel('Heating rate [eV/cm^3/s]');
ylabel('Altitude [km]');

Excitation

alt_bin_edge = linspace(50e3,1100e3,100)
alt  = (alt_bin_edge[1:] + alt_bin_edge[:-1])/2
excitation_counts = zeros(len(alt))
for i in range(len(alt)):
    idx = (excited_alt < alt_bin_edge[i+1]) & (excited_alt >= alt_bin_edge[i])
    excitation_counts[i] = sum(idx)

# Scale to a 1 erg/cm^3/s source
e_model = e_0*N_particles # particle energy * number of particles
d_alt = alt[2]-alt[1]
prod_rate = excitation_counts/d_alt*6.2415e11/e_model*1e-2 # cm^-3 s^-1. See notes.
    
#heat_rate[heat_rate==0.0] = 1e-12 # because

semilogx(prod_rate,alt/1e3,'k.-')

xlabel('O(1D) production rate [$cm^{-3} s^{-1}$]');
ylabel('Altitude [km]');

brightness_R = sum(prod_rate)*d_alt*1e2*1e4*1e-10
print 'Column brightness if all of those O(1D) immediately emitted: %f R' % brightness_R

Column brightness if all of those O(1D) immediately emitted: 5604.436552 R

import subprocess
import shutil
nb_name = 'Model_v0.ipynb' # There's probably a way to get this programmatically
cmd = ['ipython', 'nbconvert', '--to', 'slides', nb_name, \
       '--reveal-prefix', '"http://cdn.jsdelivr.net/reveal.js/2.6.2"',\
      '--template','output_toggle']
# The --template command points to a custom template which opens and closes the input code cells with a click
subprocess.call(cmd)
html_name = '%s.slides.html' % (nb_name.split('.')[0])
shutil.copy(html_name, '/home/bhardin2/public_html/slideshows/')