%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
import mpld3

# 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

Comparison of cross sections from various sources

Ishimoto 1992

def p_i92(tau):
    '''
    tau: eV-deg
    p: cm^2/rad
    '''
    p0 = 6e-17
    tau0 = 1e3
    tau1 = 200e3
    b1 = 0.75
    b2 = 2.
    
    if tau < tau0:
        p = p0
    elif tau < tau1:
        p = p0 * tau0**b1 * tau**(-b1)
    else:
        p = p0 * tau0**b1 * tau1**(-b1+b2) * tau**(-b2)
        
    return p

def dscs_i92(e,th):
    '''
    Return DSCS*sin(th)
    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(-1, 8, 100)
p = zeros(len(tau))
for i in range(len(tau)):
    p[i] = p_i92(tau[i])

loglog(tau,p,'k-')
grid()
xlabel('tau [eV-deg]')
ylabel('p [cm^2/rad]')
title('Ishimito 92 O-O DSCS')

<matplotlib.text.Text at 0x7f37ca0ca590>

def F(th_rad, e, th_min_rad=1e-9):
    '''
    un-normalized CDF of scattering angle at energy e, calculated analytically (see notebook).
    F(inf, e) is the total cross section.
    F(th, e)/F(inf, e) is the CDF of the scattering angle.
    INPUTS:
        th: angle (rad)
        e: energy (eV)
    OUTPUTS:
        F: un-normalized CDF of scattering angle.
    '''
    p0 = 6e-17
    tau0 = 1e3
    tau1 = 200e3
    b1 = 0.75
    b2 = 2.0
    
    # Calculate intermediate variables
    p1 = p0 * tau0**b1 * e**(-b1)
    p2 = p0 * tau0**b1 * tau1**(-b1+b2) * e**(-b2)
    th0 = tau0/e
    th1 = tau1/e
    
    th = th_rad * 180/pi # convert to degrees
    th_min = th_min_rad * 180/pi # convert to degrees
    
    if th_min > th0:
        raise Exception('Error: th_min > th0, violating assumption')
    
    if th < th0:
        F = p0*log(th/th_min)
    elif th < th1:
        F = p0*log(th0/th_min) + p1/b1*(th0**(-b1) - th**(-b1))
    else: # th >= th1
        F = p0*log(th0/th_min) + p1/b1*(th0**(-b1) - th1**(-b1)) + p2/b2*(th1**(-b2) - th**(-b2))
    
    # Don't forget factor of 2*pi for azimuthal integral
    return 2*pi*F
    

def elastic_xsec(e, th_min = 1e-9):
    return F(pi ,e, th_min)

th_min = 1e-6
e_min = 1e-1
e_i92 = logspace(log10(e_min),6,101)
elastic_xsec_i92 = zeros(len(e_i92))

for i in range(len(e_i92)):
    e = e_i92[i]
    elastic_xsec_i92[i] = elastic_xsec(e, th_min)

loglog(e_i92,elastic_xsec_i92,'k-')

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

def dscs_i92(e,th):
    '''
    Return 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
    

Lewkow and Kharchenko 2014

def p_lk14(tau):
    '''
    Return the reduced coordinates from Lewknow and Kharchenko 2014.
    Note p = th*sin(th)*dsig/dOmega [cm^2/rad]
         tau = e*th [eV-deg]
    '''
    mO = 15.9994 # amu
    mu = mO*mO/(mO+mO) # for O-O collisions
    gamma = 1. # for atom-atom collisions
    
    # Create variable x in the units used in LK14 [ev-deg]
    x = tau/mu
    
    # Coefficients from the model 
    C1 = -0.13
    C2 = 1.0
    C3 = 2.7
    C4 = 10.0
    C5 = 2.04
    C6 = -0.03
    C7 = 32.3
    x0 = 50.12
    
    logx = log(x)
    if x >= x0:
        p = gamma * exp(C1*logx**2 + C2*logx + C3)
    else:
        p = gamma * C4*exp(C5 + C6*logx) + C7
    
    # dsig/dOmega is in units of a0^2. Change it to cm^2
    a0 = 5.292e-11 # m
    p = p*a0**2*1e4
    # th is in units of deg. Change it to rad
    p = p*pi/180.
    
    return p

p2 = zeros(len(tau))
for i in range(len(tau)):
    p2[i] = p_lk14(tau[i])
    
loglog(tau,p,'k-',label='I92')
loglog(tau,p2,'g-',label='LK14')
#grid()
xlabel('tau [eV-deg]')
ylabel('p [cm^2/rad]')
title('DSCS in reduced coordinates')
legend(loc='best', framealpha=1)

ylim((1e-20,1e-16))
#mpld3.display()

(1e-20, 1e-16)

def dscs_lk14(e,th):
    '''
    Return DSCS*sin(th)
    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_lk14(tau)
    s = p/th
    return s
    

# Calculate total cross section numerically
th_min = 1e-6
e_min = 1e-1
e_lk14 = logspace(log10(e_min),6,101)
elastic_xsec_lk14 = zeros(len(e_i92))
N = 10000 # number of angular bins for integration

th_edge = logspace(log10(th_min),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_lk14)):
    e = e_lk14[i]
    jstart = sum(th_edge < th_min)
    
    s_cdf = zeros(len(th_edge))
    for j in range(jstart,len(th_mid)):
        s = 2*pi*dscs_lk14(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_lk14[i] = s_cdf[-1]

Kharchenko et al 2000

mO_const = 15.9994 * 1.660538921e-27 # oxygen mass in kg
e_const = 1.602e-19 # elementary charge
from scipy.special import gamma

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
    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))
    v = sqrt(2*energy_J/mO_const)
    sigma = c0 * (Cn_SI/(v*sqrt(np.pi)*hbar))**(2/(n-1))
    return sigma*1e4

e_k00 = linspace(0.1,5,100)
elastic_xsec_k00 = zeros(len(e_k00))
for i in range(len(e_k00)):
    elastic_xsec_k00[i] = sigma_LS(e_k00[i])

Shematovich 2005 extrapolation

def s05_extrap(e):
    
    # Point to extrapolate from
    e1 = 5
    f1 = sigma_LS(e1)
    
    # extrapolation exponent
    b = 0.5
    return f1 * e1**b * e**(-b)

e_s05 = logspace(log10(5), 6)
elastic_xsec_s05 = zeros(len(e_s05))
for i in range(len(e_s05)):
    elastic_xsec_s05[i] = s05_extrap(e_s05[i])

Tully and Johnson 2001

xy = array([
        9.91625e-1, 2.35300e-15,
        1.27251e+0, 2.33477e-15,
        1.51681e+0, 2.28091e-15,
        1.83437e+0, 2.22830e-15,
        2.11526e+0, 2.21103e-15,
        2.48172e+0, 2.12668e-15,
        3.01905e+0, 2.11020e-15,
        3.42210e+0, 2.06153e-15,
        4.06940e+0, 1.99837e-15,
        4.66783e+0, 1.95227e-15,
        5.24182e+0, 1.92213e-15,
        5.96239e+0, 1.86325e-15,
        6.59744e+0, 1.83448e-15,
        7.30576e+0, 1.80616e-15,
        8.71018e+0, 1.75083e-15,
        9.45524e+0, 1.65804e-15,
        1.04933e+1, 1.64519e-15,
        ])


e_tj01 = xy[0::2]
elastic_xsec_tj01 = xy[1::2]

Kozyra et al 1982

e_k82 = [100, 400, 1000, 2000]
elastic_xsec_k82 = [3.6e-15, 2.26e-15, 1.68e-15, 1.36e-15]

Including all of the hard-sphere values used by various sources

hs_vals = [1.6e-16, 7.0e-16, 4.5e-16, 2e-15, 2.2e-15, 1.7e-15]

loglog(e_k00, elastic_xsec_k00, 'g-', label='Kharchenko 2000')
loglog(e_s05, elastic_xsec_s05,'g--', dashes=(3,3), label='Shematovich 2005 extrap.')
loglog(e_i92, elastic_xsec_i92, 'k-', label='Ishimoto 1992')
loglog(e_lk14, elastic_xsec_lk14, 'c-', label='Lewkow and Kharchenko 2014')
loglog(e_k82, elastic_xsec_k82, 'b-', label='Kozyra 1982')
loglog(e_tj01,elastic_xsec_tj01,'r-', label='Tully and Johnson 2001')
for h in hs_vals:
    loglog([e_i92[0], e_i92[-1]], [h, h],'-',color=[0.5, 0.5, 0.5], linewidth=0.5)
loglog([],[],'-',color=[0.5, 0.5, 0.5], linewidth=0.5, label='Hard Sphere (various sources)')

xlabel('Energy [eV]');
ylabel('Elastic Scattering Cross section [cm$^2$]');
legend(loc='best',framealpha=1,prop={'size':6});

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

p2 = zeros(len(tau))
for i in range(len(tau)):
    p2[i] = p_lk14(tau[i])
    
loglog(tau,p,'k-',label='I92')
loglog(tau,p2,'g-',label='LK14')
#grid()
xlabel(r'$\tau = \Theta E\,\,[eV\,deg]$')
ylabel('$p = \Theta\,\sin\Theta\,d\sigma/d\Omega\,\,[cm^2/rad]$')
title('DSCS in reduced coordinates')
legend(loc='best', framealpha=1)

ylim((1e-20,1e-16))
#mpld3.display()

(1e-20, 1e-16)

Angular dependence

Kharchenko 2000

# Using http://arohatgi.info/WebPlotDigitizer/
xy = array([
        1.2246666941841795, 1.000855744891543e-13,
1.529904651834734, 7.593522900975015e-14,
1.842572336649198, 5.48706879583845e-14,
2.4183761918520688, 3.965513365753816e-14,
3.0065629256614557, 2.6422093596124545e-14,
3.3316134890824287, 1.7602450061832085e-14,
4.170553514673713, 1.2723132440618299e-14,
5.016923267428908, 8.758712127215306e-15,
6.126429190572502, 6.030431263942718e-15,
6.96041606472118, 4.5028323787052665e-15,
8.591241177198995, 3.2035699389923968e-15,
10.75576835761754, 2.1713595621314264e-15,
12.371734015767537, 1.7030587102930312e-15,
13.714657200643913, 1.4252745913272884e-15,
15.343005737400418, 1.0306357744207715e-15,
15.898996986832877, 8.482622737801527e-16,
17.72361414950263, 9.514177263037084e-16,
19.615098856647542, 6.880814079921529e-16,
22.539315639575676, 5.670508319995892e-16,
24.140421843397874, 4.903083745120934e-16,
25.220208857885822, 4.1027619565840217e-16,
27.325298220993105, 4.1074438475599955e-16,
30.25199157964255, 3.330393596194494e-16,
31.81099599620258, 3.796059190585858e-16,
33.918561935031164, 3.739125567899173e-16,
35.82738267222521, 2.4134091166187313e-16,
37.63961695628845, 2.9360485948753154e-16,
40.58859949642961, 2.0566623602848476e-16,
42.9593016056466, 2.026105188435897e-16,
44.82354397985719, 1.7521497547158951e-16,
46.92615676724314, 1.782890940338306e-16,
49.05848846328476, 1.4927237558457384e-16,
52.23593511371611, 1.3129751022070246e-16,
54.06302885210715, 1.4489030821818414e-16,
55.92727122631774, 1.2529927836628348e-16,
58.03731374086763, 1.2143039376360466e-16,
59.11214760391299, 1.0496654284998435e-16,
60.92685846369753, 1.2563913660209578e-16,
62.76881165641639, 1.2576457976915665e-16,
66.48739010195237, 1.0037135202128896e-16,
68.32934329467123, 1.0047156681598037e-16,
71.48202418788955, 1.0396870373631528e-16,
76.53857266685928, 7.173638114704152e-17,
82.57832170718619, 7.806470076078735e-17,
84.96883642217361, 6.752860241773791e-17,
87.3494448342758, 6.233813241041422e-17,
93.13101085565691, 6.565845000647364e-17,
98.42097659635942, 5.5066600339123645e-17,
105.5182234696826, 5.804097190944397e-17,
113.4345977628266, 5.035793822351166e-17,
121.33611260164281, 4.816711670267666e-17,
130.5384488380732, 5.0826754172089004e-17,
138.96375944194492, 4.942619238112e-17,
147.13088702687082, 4.6520410486820796e-17,
155.03487844140832, 4.377921870472862e-17,
163.1821934205638, 4.6926732581915626e-17,
167.14904858216042, 4.1293634139776956e-17,
169.72840219589716, 5.817351508655542e-17,
171.2576877038015, 8.058637286882847e-17,
173.09964089652038, 8.066683354460272e-17,
175.00103190655054, 5.4667678429544413e-17,
176.27708754695175, 7.09532112711193e-17,
176.76621125190903, 9.056689879921133e-17,
177.50113509720558, 1.2954973762322165e-16,
177.99025880216286, 1.653613383607691e-16,
])
th_deg = xy[0::2]
dcs = xy[1::2]
semilogy(th_deg,dcs,'k.-')
xlabel('Scattering angle [deg]')
ylabel('Differential Cross Section [cm$^2$]')
title('Direct from paper @ 3 eV')

<matplotlib.text.Text at 0x7f37c8410710>

# See if I understand this enough to recreate total cross section from the paper
th_deg_edge = th_deg
th_deg_mid = (th_deg_edge[1:] + th_deg_edge[:-1])/2
dscs_mid = (dcs[1:] + dcs[:-1])/2
th_deg_width = th_deg_edge[1:] - th_deg_edge[:-1]
#sig = 360*sum(dscs_mid*sin(th_deg_mid*pi/180)*th_deg_width)
sig = 2*pi*sum(dscs_mid*sin(th_deg_mid*pi/180)*th_deg_width*pi/180)
print 'From angular distr: %e' % sig
print 'From LS approx: %e' % sigma_LS(3)

From angular distr: 2.485879e-15
From LS approx: 5.977266e-15

th_k00 = th_deg*pi/180
dcs_k00 = dcs

Try extrapolating lower values to see if that eliminates discrepancy.

# Extrapolate so it's linear in log-log scale, which is what it looks like.
xe = log10(th_deg_edge[0:20])
ye = log10(dcs[0:20])

p = polyfit(xe,ye,1)
print p
new_x = logspace(-9,0,1000)
new_logy = p[0]*log10(new_x) + p[1]
new_y = 10**new_logy
th_deg = concatenate((new_x, th_deg))
dcs = concatenate((new_y, dcs))

loglog(th_deg,dcs,'k.-')
loglog(new_x,new_y,'r.')
xlabel('Scattering angle [deg]')
ylabel('Differential Cross Section [cm$^2$]')
title('Direct from paper @ 3 eV, extrapolated')
grid()

[ -1.84397639 -12.76750146]

# See if I understand this enough to recreate total cross section from the paper
th_deg_edge = th_deg
th_deg_mid = (th_deg_edge[1:] + th_deg_edge[:-1])/2
dscs_mid = (dcs[1:] + dcs[:-1])/2
th_deg_width = th_deg_edge[1:] - th_deg_edge[:-1]
#sig = 360*sum(dscs_mid*sin(th_deg_mid*pi/180)*th_deg_width)
sig = 2*pi*sum(dscs_mid*sin(th_deg_mid*pi/180)*th_deg_width*pi/180)
print 'From angular distr: %e' % sig
print 'From LS approx: %e' % sigma_LS(3)

From angular distr: 4.563837e-15
From LS approx: 5.977266e-15

th_k00_extrap = new_x*pi/180
dcs_k00_extrap = new_y

Tully and Johnson

xy = array([
0, 4.0400150105702295e-15,
0.06798866855524088, 3.7548157629342325e-15,
0.14447592067988668, 3.1256260230707142e-15,
0.21246458923512757, 2.3302256951101475e-15,
0.2719546742209632, 1.69514008355028e-15,
0.3356940509915015, 1.1599118609927864e-15,
0.42067988668555245, 5.772769064372794e-16,
0.47167138810198295, 3.9017714689928886e-16,
0.518413597733711, 2.838225075042674e-16,
0.5651558073654392, 2.1949649073319997e-16,
0.6288951841359773, 1.895443341574557e-16,
0.7478753541076486, 1.7405548841754669e-16,
0.8116147308781869, 1.6176547663945786e-16,
0.8796033994334278, 1.4315776175066576e-16,
1.0028328611898016, 1.0545196290356566e-16,
1.0920679886685551, 8.461906176116692e-17,
1.1643059490084984, 7.580955488962529e-17,
1.2407932011331444, 7.046032934033348e-17,
1.3597733711048159, 6.795136414213667e-17,
1.4745042492917844, 6.473289164039662e-17,
1.5509915014164302, 6.165723251009937e-17,
1.6189801699716713, 5.801079671756118e-17,
1.7167138810198301, 5.2617344385747585e-17,
1.8441926345609063, 4.4895831078520774e-17,
1.9546742209631727, 4.173371164731576e-17,
2.073654390934844, 3.9757711049878946e-17,
2.1756373937677056, 3.9290091935407075e-17,
2.2776203966005664, 3.979083152683433e-17,
2.3668555240793197, 3.980533041031056e-17,
2.456090651558073, 3.933510398321281e-17,
2.5368271954674215, 3.7928529531312376e-17,
2.61756373937677, 3.6127054963805535e-17,
2.694050991501416, 3.4410547271434234e-17,
2.76628895184136, 3.277502755727129e-17,
2.838526912181303, 3.2385603768825294e-17,
2.927762039660056, 3.2003027701935424e-17,
3.0424929178470252, 3.2812006496186546e-17,
3.1274787535410766, 3.282339302284171e-17  ,
    ])
th = xy[0::2]
dcs = xy[1::2]
semilogy(th,dcs,'k.-')
xlabel('Scattering angle [rad]')
ylabel('Cross Section [cm$^2$]')
title('Direct from paper @ 5 eV [units issue?]')

<matplotlib.text.Text at 0x7f37c8771610>

# See if I understand this enough to recreate total cross section from the paper
th_edge = th
th_mid = (th_edge[1:] + th_edge[:-1])/2
dscs_mid = (dcs[1:] + dcs[:-1])/2
th_width = th_edge[1:] - th_edge[:-1]
#sig = 360*sum(dscs_mid*sin(th_deg_mid*pi/180)*th_deg_width)
#sig = 2*pi*sum(dscs_mid*th_width)
sig = 2*pi*sum(dscs_mid*sin(th_mid)*th_width)
print 'From angular distr: %e' % sig
print 'From paper: %e' % 1.9e-15

From angular distr: 1.935262e-15
From paper: 1.900000e-15

# Save for later
th_tj01 = th
#th_tj01[0] = 1e-9 # Because 0.0 won't work on a log-log plot
dcs_tj01 = dcs

Compare three papers

e = 5
th_edge = logspace(-4,log10(pi),1000)
th_mid = (th_edge[1:]+th_edge[:-1])/2
th_width = th_edge[1:] - th_edge[:-1]
dcs_i92_5ev = zeros(len(th_mid))
dcs_lk14_5ev = zeros(len(th_mid))
for i in range(len(th_mid)):
    dcs_i92_5ev[i] = dscs_i92(e,th_mid[i])/sin(th_mid[i])
    dcs_lk14_5ev[i] = dscs_lk14(e,th_mid[i])/sin(th_mid[i])

th_i92 = th_mid

e = 3
#th_edge = logspace(-9,log10(pi),1000)
#th_mid = (th_edge[1:]+th_edge[:-1])/2
#th_width = th_edge[1:] - th_edge[:-1]
dcs_i92_3ev = zeros(len(th_mid))
dcs_lk14_3ev = zeros(len(th_mid))
for i in range(len(th_mid)):
    dcs_i92_3ev[i] = dscs_i92(e,th_mid[i])/sin(th_mid[i])
    dcs_lk14_3ev[i] = dscs_lk14(e,th_mid[i])/sin(th_mid[i])

th_i92 = th_mid

dcs_i92_3ev_w = 2*pi * dcs_i92_3ev * sin(th_i92) # This is closer to a pdf of scattering angle
dcs_i92_5ev_w = 2*pi * dcs_i92_5ev * sin(th_i92) # This is closer to a pdf of scattering angle
dcs_lk14_3ev_w = 2*pi * dcs_lk14_3ev * sin(th_i92) # This is closer to a pdf of scattering angle
dcs_lk14_5ev_w = 2*pi * dcs_lk14_5ev * sin(th_i92) # This is closer to a pdf of scattering angle
dcs_k00_w = 2*pi * dcs_k00 * sin(th_k00)
dcs_tj01_w = 2*pi * dcs_tj01 * sin(th_tj01)

semilogy(th_k00, dcs_k00_w, 'g-', label='Kharchenko 2000 3eV')
semilogy(th_i92, dcs_i92_3ev_w, 'g--',dashes=(3,3), label='Ishimoto 1992 3eV')
semilogy(th_i92, dcs_lk14_3ev_w, 'g--',dashes=(1,1), label='Lewkow and Kharchenko 2014 3eV')
#semilogy(th_tj01, dcs_tj01_w, 'r-', label='Tully and Johnson 2001 5eV')
#semilogy(th_i92, dcs_i92_5ev_w, 'r--',dashes=(3,3), label='Ishimoto 1992 5eV')
#semilogy(th_i92, dcs_lk14_5ev_w, 'r--',dashes=(1,1), label='Lewkow and Kharchenko 2014 5eV')
xlabel('Scattering angle [rad]')
ylabel('DCS $d\sigma / d \Theta$ [cm$^2$/rad]')
grid()
legend(loc='best', framealpha=1)
ylim((1e-17,1e-13))

(1e-17, 1e-13)

dcs_i92_3ev_w = 2*pi * dcs_i92_3ev * sin(th_i92) # This is closer to a pdf of scattering angle
dcs_i92_5ev_w = 2*pi * dcs_i92_5ev * sin(th_i92) # This is closer to a pdf of scattering angle
dcs_lk14_3ev_w = 2*pi * dcs_lk14_3ev * sin(th_i92) # This is closer to a pdf of scattering angle
dcs_lk14_5ev_w = 2*pi * dcs_lk14_5ev * sin(th_i92) # This is closer to a pdf of scattering angle
dcs_k00_w = 2*pi * dcs_k00 * sin(th_k00)
dcs_tj01_w = 2*pi * dcs_tj01 * sin(th_tj01)

#semilogy(th_k00, dcs_k00_w, 'g-', label='Kharchenko 2000 3eV')
#semilogy(th_i92, dcs_i92_3ev_w, 'g--',dashes=(3,3), label='Ishimoto 1992 3eV')
#semilogy(th_i92, dcs_lk14_3ev_w, 'g--',dashes=(1,1), label='Lewkow and Kharchenko 2014 3eV')
semilogy(th_tj01, dcs_tj01_w, 'r-', label='Tully and Johnson 2001 5eV')
semilogy(th_i92, dcs_i92_5ev_w, 'r--',dashes=(3,3), label='Ishimoto 1992 5eV')
semilogy(th_i92, dcs_lk14_5ev_w, 'r--',dashes=(1,1), label='Lewkow and Kharchenko 2014 5eV')
xlabel('Scattering angle [rad]')
ylabel('DCS $d\sigma / d \Theta$ [cm$^2$/rad]')
grid()
legend(loc='best', framealpha=1)
ylim((1e-17,1e-13))

(1e-17, 1e-13)

dcs_i92_3ev_w = dcs_i92_3ev * sin(th_i92) # This is closer to a pdf of scattering angle
dcs_i92_5ev_w = dcs_i92_5ev * sin(th_i92) # This is closer to a pdf of scattering angle
dcs_k00_w = dcs_k00 * sin(th_k00)
dcs_k00_extrap_w = dcs_k00_extrap * sin(th_k00_extrap)
dcs_tj01_w = dcs_tj01 * sin(th_tj01)

loglog(th_k00, dcs_k00_w, 'g-', label='Kharchenko 2000 3eV')
loglog(th_k00_extrap, dcs_k00_extrap_w, 'g-',linewidth=0.5, label='Kharchenko 2000 3eV (extrap)')
loglog(th_i92, dcs_i92_3ev_w, 'g--',dashes=(3,3), label='Ishimoto 3eV')
loglog(th_tj01, dcs_tj01_w, 'r-', label='Tully and Johnson 2001 5eV')
loglog(th_i92, dcs_i92_5ev_w, 'r--',dashes=(3,3), label='Ishimoto 1992 5eV')
xlabel('Scattering angle [rad]')
ylabel('DCS $d\sigma / d \Theta$ [cm$^2$/rad]')
grid()
legend(loc='best', framealpha=1)
xlim((1e-4,1e1))
ylim((1e-19,1e-12))

(1e-19, 1e-12)

Show Ishimoto $d \sigma / d\theta$ at various energies

e_vec = [1,10,100,1000,10000]
th_edge = logspace(-9,log10(pi),1000)
th_mid = (th_edge[1:]+th_edge[:-1])/2
th_width = th_edge[1:] - th_edge[:-1]

for j in range(len(e_vec)):
    e = e_vec[j]

    dcs_i92 = zeros(len(th_mid))
    for i in range(len(th_mid)):
        dcs_i92[i] = dscs_i92(e,th_mid[i])/sin(th_mid[i])

    dcs_i92_w = 2*pi * dcs_i92 * sin(th_mid) # This is closer to a pdf of scattering angle

    semilogy(th_mid, dcs_i92_w, '-', label='%i eV' % e)
xlabel('Scattering angle [rad]')
ylabel('DCS $d\sigma / d \Theta$ [cm$^2$/rad]')
grid()
legend(loc='best', framealpha=1)
ylim((0.5*min(dcs_i92_w),1e-13))
title('Ishimoto 92 model')
clear_output()

e_vec = [1,10,100,1000,10000]
th_edge = logspace(-4,log10(pi),1000)
th_mid = (th_edge[1:]+th_edge[:-1])/2
th_width = th_edge[1:] - th_edge[:-1]

for j in range(len(e_vec)):
    e = e_vec[j]

    dcs_i92 = zeros(len(th_mid))
    for i in range(len(th_mid)):
        dcs_i92[i] = dscs_i92(e,th_mid[i])/sin(th_mid[i])

    dcs_i92_w = 2*pi * dcs_i92 * sin(th_mid) # This is closer to a pdf of scattering angle

    loglog(th_mid, dcs_i92_w, '-', label='%i eV' % e)
xlabel('Scattering angle [rad]')
ylabel('DCS $d\sigma / d \Theta$ [cm$^2$/rad]')
grid()
legend(loc='best', framealpha=1)
#ylim((0.5*min(dcs_i92_w),1e-13))

<matplotlib.legend.Legend at 0x7f37c7f4f110>

Compare Ishimoto 1992 and Lewkow 2014

e_vec = [1,1000,10000]
th_edge = logspace(-9,log10(pi),1000)
th_mid = (th_edge[1:]+th_edge[:-1])/2
th_width = th_edge[1:] - th_edge[:-1]
cols = ['k','g','r','b','m','y']

for j in range(len(e_vec)):
    e = e_vec[j]

    dcs_i92 = zeros(len(th_mid))
    dcs_lk14 = zeros(len(th_mid))
    for i in range(len(th_mid)):
        dcs_i92[i] = dscs_i92(e,th_mid[i])/sin(th_mid[i])
        dcs_lk14[i] = dscs_lk14(e,th_mid[i])/sin(th_mid[i])

    dcs_i92_w = 2*pi * dcs_i92 * sin(th_mid) # This is closer to a pdf of scattering angle
    dcs_lk14_w = 2*pi * dcs_lk14 * sin(th_mid) # This is closer to a pdf of scattering angle

    semilogy(th_mid, dcs_i92_w, '%s--'%cols[j], label='I92 %i eV' % e, dashes=(2,2))
    semilogy(th_mid, dcs_lk14_w, '%s-'%cols[j], label='LK14 %i eV' % e)
    
xlabel('Scattering angle [rad]')
ylabel('DCS $d\sigma / d \Theta$ [cm$^2$/rad]')
grid()
legend(loc='best', framealpha=1)
ylim((0.5*min(dcs_lk14_w),1e-13))

(4.5108787486253287e-21, 1e-13)

e_vec = [1,1000,10000]
th_edge = logspace(-4,log10(pi),1000)
th_mid = (th_edge[1:]+th_edge[:-1])/2
th_width = th_edge[1:] - th_edge[:-1]
cols = ['k','g','r','b','m','y']

for j in range(len(e_vec)):
    e = e_vec[j]

    dcs_i92 = zeros(len(th_mid))
    dcs_lk14 = zeros(len(th_mid))
    for i in range(len(th_mid)):
        dcs_i92[i] = dscs_i92(e,th_mid[i])/sin(th_mid[i])
        dcs_lk14[i] = dscs_lk14(e,th_mid[i])/sin(th_mid[i])

    dcs_i92_w = 2*pi * dcs_i92 * sin(th_mid) # This is closer to a pdf of scattering angle
    dcs_lk14_w = 2*pi * dcs_lk14 * sin(th_mid) # This is closer to a pdf of scattering angle

    loglog(th_mid, dcs_i92_w, '%s--'%cols[j], label='I92 %i eV' % e, dashes=(2,2))
    loglog(th_mid, dcs_lk14_w, '%s-'%cols[j], label='LK14 %i eV' % e)
    
xlabel('Scattering angle [rad]')
ylabel('DCS $d\sigma / d \Theta$ [cm$^2$/rad]')
grid()
legend(loc='best', framealpha=1)
#ylim((0.5*min(dcs_lk14_w),1e-13))

<matplotlib.legend.Legend at 0x7f37c9f51850>

import subprocess
import shutil
nb_name = 'Model_v0_cross_section_comparisons.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/')