''' Step by step version of ENA spectra calculation
'''
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 scipy.interpolate
import pyana.pep.iotorus
def get_density(iotorus, x, y, z):
''' An adapater of the torus density getter.
As iotorus.get_density will get Rj distance returning value with unit of per cm3, so
this wraps to the MKSA system.
'''
xrj = x / (7.15e4 * 1e3)
yrj = y / (7.15e4 * 1e3)
zrj = z / (7.15e4 * 1e3)
dens_per_cm3 = iotorus.get_density(xrj, yrj, zrj)
return dens_per_cm3 * 1e6 # Return the density in the unit of per m3.
def get_corotational_density(r):
''' Corotational plasma density model.
r to be m. Returned to be m^-3
Electron density is used to represent the "single fluid".
'''
r_rj = np.array([2, 4, 5., 6, 7, 9, 11, 13, 15, 20, 25]) * (7.15e4 * 1e3) # in m
n_cp_rj = np.log10(np.array([1e2/17/17/17, 1e2 / 17, 1e2, 1.7e3, 7.2e2, 4.3e1, 9.9, 3.9, 1.7, .23, .13])) + 6 # in m^-3
n_cp = scipy.interpolate.interp1d(r_rj, n_cp_rj)
return 10 ** n_cp(r)
rj2m = 7.15e4 * 1e3
m2rj = 1 / 7.15e4 / 1e3
def getenaflux(sc=[5.9, 13, 0], los=[0, -1, 0], ene0=100, ene1=200, dl_rj=0.01, rj1=50, iotorus_model=None):
'''
'''
logger = logging.getLogger('getenaflux')
logger.setLevel(logging.INFO)
# Energy range. 100--200 eV, oxygen assumed.
v0 = 13.8 * np.sqrt(ene0) / 4. * 1e3
v1 = 13.8 * np.sqrt(ene1) / 4. * 1e3
vc = (v0 + v1) / 2.
# First, the neutral torus is instanced, if needed.
if iotorus_model == None:
iotorus = pyana.pep.iotorus.TwoPeakModel()
else:
iotorus = iotorus_model
# Then, the spacecraft position is given.
# The frame used is the io fixed frame.
# This means that the Jupiter is at the center
# and the Io is at (5.9Rj, 0, 0).
# Let's use MKSA system internally always.
# scpos = np.array([15.3, 0, 0]) # At ganymede orbit.
scpos = np.array(sc) # Around ganymede, direct looking Io along -y.
scpos = scpos * 7.15e4 * 1e3 # in meter.
logger.debug('scpos: %s' % str(scpos))
# Looking direction
# lookdir = np.array([-1, 0, 0]) # Should be unit vector
lookdir = np.array(los) # Should be unit vector
lookdir = lookdir / np.sqrt((lookdir ** 2).sum())
logger.debug('lookdir %s' % str(lookdir))
# Density along the looking direction.
dl = dl_rj * rj2m
ls_rj = np.arange(0, rj1, dl_rj) #ls_rj is the distance from the spacecraft in Rj, (N, )
ls = ls_rj * (7.15e4 * 1e3) #ls is the distance from the spacecraft in m, (N, )
n_los = np.zeros_like(ls) # n_los stores the density at the ls along los., (N, )
n_corpla = np.zeros_like(ls) # n_corpla is the density of the plasma
enaflux = 0.
# For each segment along the LOS
for il, l in enumerate(ls):
lpos = scpos + lookdir * l # lpos is the vector of the LOS=l.
n_at_l = get_density(iotorus, lpos[0], lpos[1], lpos[2])
n_los[il] = n_at_l
lpos_r = np.sqrt((lpos ** 2).sum()) # Radius in m, from Jupietr
if lpos_r < 1 * rj2m: break
# Corotational plasma density depending on the radius
if lpos_r < 2 * rj2m or lpos_r > 25 * rj2m:
# Too close to Jupiter or too away.
n_corpla[il] = 0.
else:
n_corpla[il] = get_corotational_density(lpos_r)
# Corotational plasma velocity
sys3per = 9.9250 * 3600 # sec
sys3omega = 2 * np.pi / sys3per
vspeed = sys3omega * np.sqrt(lpos[0] * lpos[0] + lpos[1] * lpos[1])
v_corpla = np.array([-lpos[1], lpos[0], 0])
v_corpla = v_corpla / np.sqrt((v_corpla ** 2).sum()) * vspeed
# Angle between corotational plasma flow and sensor direction
tosensor = -lookdir
corpla_norm = v_corpla / vspeed
theta = np.arccos(np.clip((tosensor * corpla_norm).sum(), -1, 1))
# Corotational plasma temperature (100eV assumed)
vth = np.sqrt(100. * 1.6e-19 / 1.67e-27 / 16.)
# Maxwell distribution (polar coordinate)
maxwell = lambda v: n_at_l * ((2 * np.pi * (vth**2)) ** -1.5) * v * v * np.exp(- v * v * (np.sin(theta) ** 2) / (vth **2)) * np.exp(-((v * np.cos(theta) - vspeed) / vth) ** 2)
fmaxwell = maxwell(vc) * (v1 - v0) # This is "equivalent" to energy integral from E0 to E1.
dflux = fmaxwell * vc # Now the differential flux contribution of this segment
xsect = 1e-15 * 1e-4 # in m^2
dflux_ena = dflux * n_at_l * xsect * dl
enaflux += dflux_ena
# print >> sys.stderr, lpos[0], lpos[1], lpos[2], l, '(', l * m2rj, 'Rj)', n_at_l, n_corpla[il], v_corpla, vth, np.rad2deg(theta), fmaxwell, dflux, dflux_ena
# ENA flux is in the unit of #/m2 sr s.
logger.debug('%e' % (enaflux * 1e-4))
# plt.plot(ls_rj, n_los)
# plt.plot(ls_rj, n_corpla)
return enaflux * 1e-4
def main():
scpos = [0, 15.9 * np.cos(np.deg2rad(20)), 15.9 * np.sin(np.deg2rad(20))]
print('POS', scpos[0], scpos[1], scpos[2])
thetalist = list(range(0, 91, 3))
philist = list(range(0, 360, 3))
nall = len(thetalist) * len(philist)
ndone = 0
logger = logging.getLogger('app03_enaflux.py')
logger.setLevel(logging.INFO)
import time
t0 = time.time()
for theta in thetalist:
for phi in philist:
t = np.deg2rad(theta)
p = np.deg2rad(phi)
lookdir = [np.sin(t) * np.cos(p), -np.cos(t), np.sin(t) * np.sin(p)]
flx = getenaflux(scpos, lookdir)
print(theta, phi, lookdir[0], lookdir[1], lookdir[2], flx)
ndone += 1.0
t1 = time.time()
logger.info('%d/%d, %ds. Exp %ds' % (
ndone, nall, t1 - t0, (t1 - t0) / ndone * (nall - ndone)))
if __name__ == "__main__":
main()