''' H-ENA production at Europa
'''
import os
import sys
import logging
logging.basicConfig()
import datetime
import math
import matplotlib.pyplot as plt
import numpy as np
import scipy as sp
import pyana.util.units as u
from pyana.util.sputtering import ThompsonSigmund as tsmodel
from pyana.cena import empirical as emp
def sp(estep, m1, yld, upstflx, e0):
''' Sputtering.
m1 in amu
upstflx in /cm2 s.
e0 in eV.
'''
eb = 1.0
m2 = 1 # Water.
vupst = 13.8e5 * np.sqrt(e0 / m1) # in cm/s
n0 = upstflx / vupst
print(n0, vupst)
ts = tsmodel(n0, yld, eb, e0, m1, m2)
fe = ts.f(estep) # It is energy spectra, in /cm3 eV
vlist = 13.8e5 * np.sqrt(estep / m1)
return np.clip(fe * vlist, 0, np.inf)
def bs(estep, temp_ev, refrate, upstflx):
''' Return the backscattered flux. J.
'''
corresponding_vsw = emp.guess_vsw(temp_ev * u.eV_asT) # Corresponding velocity
corresponding_vsw = corresponding_vsw.asNumber(u.km / u.s)
j = emp.j2(upstflx, corresponding_vsw)
jval = [j(e) for e in estep]
return np.array(jval)
def ech_bs(estep):
return bs(estep, 55, .25, 2.9e8)
def eco_sp(estep):
return sp(estep, 25., 1, 8.8e8, 1000)
def ehh_bs(estep):
return bs(estep, 3000, .25, 1.5e7)
def ehh_sp(estep):
return sp(estep, 1., 3., 1.5e7, 100e3)
def eho_sp(estep):
return sp(estep, 25., 50., 1e7, 100e3)
def calc_jc(func, estep):
je = func(estep)
c = je * 8e-6 * estep
return je, c
def histtypeplot(ax, xlist, ylist, x0, x1, *args, **kwds):
newx = []
newx.append(x0)
for ix in range(len(xlist) - 1):
newx.append(np.sqrt(xlist[ix] * xlist[ix + 1]))
newx.append(np.sqrt(xlist[ix] * xlist[ix + 1]))
newx.append(x1)
newy = []
for y in ylist:
newy.append(y)
newy.append(y) # Do twice.
newx = np.array(newx)
newy = np.array(newy)
ax.plot(newx, newy, *args, **kwds)
def main():
'''Main script'''
estep = np.logspace(0, 4, 1000)
estep2 = np.array([7, 11, 17, 25, 38, 57, 86, 129, 193, 290, 435, 652, 978, 1467, 2200, 3300])
e0 = 5.
e1 = 3800
fig = plt.figure(figsize=(15, 8))
axj = fig.add_subplot(121)
axc = fig.add_subplot(122)
je_ch_bs, ce_ch_bs = calc_jc(ech_bs, estep)
je_ch_bs2, ce_ch_bs2 = calc_jc(ech_bs, estep2)
axj.plot(estep, je_ch_bs, 'b', label='H+ (corotate) BS')
axc.plot(estep2, ce_ch_bs2, 'b:', label='H+ (corotate) BS')
# histtypeplot(axc, estep2, ce_ch_bs2, e0, e1, 'b', label='H+(corotate) BS')
je_co_sp, ce_co_sp = calc_jc(eco_sp, estep)
je_co_sp2, ce_co_sp2 = calc_jc(eco_sp, estep2)
axj.plot(estep, je_co_sp, 'g', label='Heavy (corotate) ISP')
axc.plot(estep2, ce_co_sp2, 'g:', label='Heavy (corotate) ISP')
# histtypeplot(axc, estep2, ce_co_sp2, e0, e1, 'g', label='Heavy(corotate) ISP')
je_hh_bs, ce_hh_bs = calc_jc(ehh_bs, estep)
je_hh_bs2, ce_hh_bs2 = calc_jc(ehh_bs, estep2)
axj.plot(estep, je_hh_bs, 'r', label='H+ (energetic) BS')
axc.plot(estep2, ce_hh_bs2, 'r:', label='H+ (energetic) BS')
# histtypeplot(axc, estep2, ce_hh_bs2, e0, e1, 'r', label='H+(energetic) BS')
# Sputtering.
je_hh_sp, ce_hh_sp = calc_jc(ehh_sp, estep)
je_hh_sp2, ce_hh_sp2 = calc_jc(ehh_sp, estep2)
axj.plot(estep, je_hh_sp, 'm', label='H+ (energetic) ISP')
axc.plot(estep2, ce_hh_sp2, 'm:', label='H+ (energetic) ISP')
# histtypeplot(axc, estep2, ce_hh_sp2, e0, e1, 'm', label='H+(energetic) ISP')
je_ho_sp, ce_ho_sp = calc_jc(eho_sp, estep)
je_ho_sp2, ce_ho_sp2 = calc_jc(eho_sp, estep2)
axj.plot(estep, je_ho_sp, 'c', label='Heavy (energetic) ISP')
axc.plot(estep2, ce_ho_sp2, 'c:', label='Heavy (energetic) ISP')
# histtypeplot(axc, estep2, ce_ho_sp2, e0, e1, 'c', label='Heavy(energetic) ISP')
# axc.plot(estep2, ce_ch_bs2+ce_co_sp2+ce_hh_bs2+ce_hh_sp2+ce_ho_sp2, 'ko', label='Total')
histtypeplot(axc, estep2, ce_ch_bs2+ce_co_sp2+ce_hh_bs2+ce_hh_sp2+ce_ho_sp2, e0, e1, 'k', label='Total')
axj.set_xscale('log')
axj.set_yscale('log')
axj.set_ylim(1e-1, 1e8)
axj.grid(True)
axj.legend(loc='best')
axj.set_ylabel('Differential Flux [/cm2 eV sr s]')
axc.set_xscale('log')
axc.set_yscale('log')
axc.set_ylim(1e-1, 1e3)
axc.set_xlim(5, 5000)
axc.grid(True)
axc.legend(loc='best')
axc.set_ylabel('Count rate [cps]')
axj.set_title('Europa')
axc.set_title('Europa')
axc.axhline(np.sqrt(450./(15/8.)), color='orange') # 15 sec
axc.text(5000, np.sqrt(450./(15./8.)), 'BG (1 scan=15s)', horizontalalignment='right', color='orange')
axc.axhline(np.sqrt(450./(60/8.)), linestyle='--', color='orange') # 60 sec
axc.text(5000, np.sqrt(450./(60./8.)), 'BG (4 scan=60s)', horizontalalignment='right', color='orange')
axc.axhline(np.sqrt(450./(300/8.)), linestyle='-.', color='orange') # 300 sec
axc.text(5000, np.sqrt(450./(300./8.)), 'BG (20 scan=5m)', horizontalalignment='right', color='orange')
print('S/N', np.max(ce_ch_bs2+ce_co_sp2+ce_hh_bs2+ce_hh_sp2+ce_ho_sp2) / np.sqrt(450./(60./8.)))
fig.savefig('henae.eps')
if __name__ == "__main__":
main()