snippet.jdc_dflux_simple

A very simple calculation of JDC energy spectra expected

A small concern about multiple charged particle, particularly where j to c conversion. (Does the conversion same as the singly charged? Where E/q plays into game? j(E) or j(E/q)?)

''' A very simple calculation of JDC energy spectra expected

A small concern about multiple charged particle,
particularly where j to c conversion.
(Does the conversion same as the singly charged?
Where E/q plays into game?
j(E) or j(E/q)?)

'''
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


def j(ene, n0=1., temp=100, m=1, q=1, vcr=150):
    '''

    ene to be energy in eV
    n0 is density in cm^-3
    temp is temperature in eV
    m is the mass in amu
    q is the charge in the unit of 1.6e-19
    '''
    k = 1.38e-23    # J/K
    mp = 1.67e-27  # kg
    qe = 1.6e-19

    m = m * mp
    q = q * qe

    vcr = vcr * 1e3    # in m/s

    ene = ene * qe   # in Joule
    n0 = n0 * 1e6         # in m^-3
    temp = temp * 1.16e4  # in K
    vth = np.sqrt(k * temp / m)

    j0 = 2 * ene / m / m * n0
    j1 = (0.5 / np.pi / (vth ** 2)) ** 1.5
    j2 = (np.sqrt(2 * ene / m) - vcr) ** 2
    j3 = 2 * (vth ** 2)

    jj = j0 * j1 * np.exp(-j2 / j3)   # in # / m2 J sr s
    jj = jj / 1e4  # in /cm2 J sr s
    jj = jj * qe   # in /cm2 eV sr s

    return jj

    
def main():
    '''Main script'''

    # Take a simple energy table.
    etbl = np.logspace(1, 4, 512)

    # For Ganymede
#    nlist = [.31, .12, .1, .1, .43, .1, .084]
    vcr = 150.

    # For Europa
    nlist = [3.20,    3.00,   2.60,   2.60,   11.00,  2.60,   2.10]
    vcr = 90.

    lablist = ['H+', 'O+', 'O++', 'S+', 'S++', 'S+++', 'Na+']
    mlist = [1, 16, 16, 32, 32, 32, 23]
    qlist = [1, 1, 2, 1, 2, 3, 1]


    for n, lab, m, q in zip(nlist, lablist, mlist, qlist):
        jp = np.zeros_like(etbl)
        for i, e in enumerate(etbl):
            jp[i] = j(e, n0=n, temp=100, m=m, q=q, vcr=vcr)
#        plt.plot(etbl / q, jp, '-', label=lab)
        ### Count rate.  J * G * E
        c = jp * 5e-4 * etbl
        plt.plot(etbl / q, c, '-', label=lab)

        print('#', lab)
        for cc, ee in zip(c, etbl):
            print(ee / q, cc)


    plt.xscale('log')
    plt.yscale('log')
    plt.ylim(1e+0, 1e6)
    plt.legend(loc='best')
    plt.xlabel('Energy [eV/q]')
#    plt.ylabel('Diff flux [#/cm2 eV sr s]')
    plt.ylabel('Count rate [cps]')


    plt.savefig('jdc_dflux_simple.eps')

if __name__ == "__main__":
    main()