apps120925_h_ena.hena

H-ENA production at Ganymede

''' H-ENA production at Ganymede
'''

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 irfpy.util.unitq as u
from irfpy.util.sputtering import ThompsonSigmund as tsmodel
from irfpy.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  # Hydrogen

    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 gch_bs(estep):
    return bs(estep, 80, .25, 1.5e7)

def gco_sp(estep):
    return sp(estep, 25., 3, 4.5e7, 3000)

def gph_bs(estep):
    return bs(estep, 9, .25, 2.5e7)

def ghh_bs(estep):
    return bs(estep, 3000, .25, 5e6)

def ghh_sp(estep):
    return sp(estep, 1., 3., 5e6, 100e3)

def gho_sp(estep):
    return sp(estep, 25., 50., 2e6, 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(gch_bs, estep)
    je_ch_bs2, ce_ch_bs2 = calc_jc(gch_bs, estep2)
#    axj.plot(estep, je_ch_bs, ':', label='H+ (corotation) BS')
#    axc.plot(estep2, ce_ch_bs2, ':', label='H+ (corotation) BS')

    je_co_sp, ce_co_sp = calc_jc(gco_sp, estep)
    je_co_sp2, ce_co_sp2 = calc_jc(gco_sp, estep2)
#    axj.plot(estep, je_co_sp, ':', label='O+ (corotation) ISP')
#    axc.plot(estep2, ce_co_sp2, ':', label='O+ (corotation) ISP')

    je_ph_bs, ce_ph_bs = calc_jc(gph_bs, estep)
    je_ph_bs2, ce_ph_bs2 = calc_jc(gph_bs, estep2)
    axj.plot(estep, je_ph_bs, 'b', label='H+ (cusp) BS')
    axc.plot(estep2, ce_ph_bs2, 'b:', label='H+ (cusp) BS')

    je_hh_bs, ce_hh_bs = calc_jc(ghh_bs, estep)
    je_hh_bs2, ce_hh_bs2 = calc_jc(ghh_bs, estep2)
    axj.plot(estep, je_hh_bs, 'r', label='H+ (energetic) BS')
    axc.plot(estep2, ce_hh_bs2, 'r:', label='H+ (energetic) BS')

    # Sputtering.

    je_hh_sp, ce_hh_sp = calc_jc(ghh_sp, estep)
    je_hh_sp2, ce_hh_sp2 = calc_jc(ghh_sp, estep2)
    axj.plot(estep, je_hh_sp, 'm', label='H+ (energetic) ISP')
    axc.plot(estep2, ce_hh_sp2, 'm:', label='H+ (energetic) ISP')

    je_ho_sp, ce_ho_sp = calc_jc(gho_sp, estep)
    je_ho_sp2, ce_ho_sp2 = calc_jc(gho_sp, estep2)
    axj.plot(estep, je_ho_sp, 'c', label='Heavy (energetic) ISP')
    axc.plot(estep2, ce_ho_sp2, 'c:', label='Heavy (energetic) ISP')

#    axc.plot(estep, ce_ph_bs+ce_hh_bs+ce_hh_sp+ce_ho_sp, 'k', label='Total')
    histtypeplot(axc, estep2, ce_ph_bs2+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-3, 1e6)
    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-2, 1e2)
    axc.set_xlim(5, 5000)
    axc.grid(True)
    axc.legend(loc='best')
    axc.set_ylabel('Count rate [cps]')

    axj.set_title('Ganymede')
    axc.set_title('Ganymede')

    axc.axhline(np.sqrt(1.97/(15./8.)), color='orange')
    axc.text(5000, np.sqrt(1.97/(15./8.)), 'BG (1 scan=15s)', horizontalalignment='right', color='orange')
    axc.axhline(np.sqrt(1.97/(60/8.)), linestyle='--', color='orange')  # 60 sec
    axc.text(5000, np.sqrt(1.97/(60./8.)), 'BG (4 scan=60s)', horizontalalignment='right', color='orange')
    axc.axhline(np.sqrt(1.97/(300/8.)), linestyle='-.', color='orange') # 300 sec
    axc.text(5000, np.sqrt(1.97/(300./8.)), 'BG (20 scan=5m)', horizontalalignment='right', color='orange')

    print('S/N ratio max', np.max(ce_ph_bs2+ce_hh_bs2+ce_hh_sp2+ce_ho_sp2) / np.sqrt(1.97/(60./8.)))

    fig.savefig('hena.eps')

if __name__ == "__main__":
    main()