''' Virtual imager of the ENA from Io torus.
The :class:`SinglePointImager` provides a virtual ENA imager
that measures ENA (differential) flux from multiple directions.
'''
import os
import sys
import time, datetime
from io import StringIO
import numpy as np
from argparse import ArgumentParser
from . import simulator
from . import observer
from pyana.util.future import OrderedDict
import pyana.util.datafile
from . import parameters as p
from . import f2j
def singlepointimager_juice_at(t, vena_sc, mass=16, primary_axis=None, secondary_axis=None):
''' Return the imager on JUICE at the time given.
'''
import pyana.juice.jspice as js
js.init()
posvel = js.get_posvel(t, frame='JEqS')
pos = posvel[:3] / p.rj
vel = posvel[3:]
io_jeqs = js.get_posvel(t, frame='JEqS', target='IO', origin='Jupiter')
iophase = np.rad2deg(np.arctan2(io_jeqs[1], io_jeqs[0]))
sun = js.sun_longitude_system3(t)
return SinglePointImager(pos, vel, vena_sc, iophase, mass=mass,
primary_axis=primary_axis, secondary_axis=secondary_axis,
sunsys3long_t0=sun)
class SinglePointImager:
''' An imager class to calculate the flux map.
'''
def __init__(self, scpos, scvel, vena_sc, iophase, mass=16,
primary_axis=None,
secondary_axis=None,
sunsys3long_t0=0):
# SC info in JEqS
self.scpos = np.array(scpos) * p.rj # in km.
self.scvel = np.array(scvel) # in km/s
if secondary_axis == None and (self.scvel ** 2).sum() == 0:
raise ValueError('Velocity is 0, so that secondary_axis should be specified')
self.sc = observer.VirtualSpacecraft(self.scpos, self.scvel, primary_axis=primary_axis, secondary_axis=secondary_axis)
self.vena_sc = vena_sc
self.mass = mass
self.pax = self.sc.primary_axis
self.sax = self.sc.secondary_axis
self.tax = self.sc.third_axis
# Io position at t=0.
self.iophase_t0 = iophase
self.iopos_t0 = 5.9 * np.array([np.cos(np.deg2rad(self.iophase_t0)), np.sin(np.deg2rad(self.iophase_t0)), 0]) * p.rj
self.sunsys3long_t0 = sunsys3long_t0
def do_calculation(self, xlist, ylist, output=None):
datf = pyana.util.datafile.Datafile()
# xlist = np.linspace(-maxangle, maxangle, resolution)
# ylist = np.linspace(-maxangle, maxangle, resolution)
xlist = np.array(xlist)
ylist = np.array(ylist)
maxangle = max([np.abs(xlist).max(), np.abs(ylist).max()])
xytot = len(xlist) * len(ylist)
xydone = 0
t0 = time.time()
datf.add_header('SCPOS', '%f %f %f' %
(self.scpos[0], self.scpos[1], self.scpos[2]))
datf.add_header('SCVEL', '%f %f %f' %
(self.scvel[0], self.scvel[1], self.scvel[2]))
datf.add_header('MASS, ENAVEL(SCf)', '%f %f' %
(self.mass, self.vena_sc))
datf.add_header('SUN_SYS3, IO', '%f %f %f %f' %
(self.sunsys3long_t0,
self.iopos_t0[0], self.iopos_t0[1], self.iopos_t0[2]))
datf.add_header('PAX', '%f % f %f' %
(self.pax[0], self.pax[1], self.pax[2]))
datf.add_header('SAX', '%f % f %f' %
(self.sax[0], self.sax[1], self.sax[2]))
datf.add_header('TAX', '%f % f %f' %
(self.tax[0], self.tax[1], self.tax[2]))
datf.add_header('MAXANGLE', '%f' % maxangle)
datf.add_header('DATAFORMAT',
'Ximg Yimg theta phi Nx Ny Nz f[s3/m6] j[/cm2_sr_eV_s]')
sim = simulator.SimulatorCore()
sim.set_iophase(self.iophase_t0)
sim.set_jupiterphase(self.sunsys3long_t0)
dl = 0.05 * p.rj # Resolution in L. Used to define tlist.
dt = dl / self.vena_sc
totl = 50 * p.rj # LOS length.
tott = totl / self.vena_sc
tlist = np.arange(0, -tott, -dt)
print('Dt = %f / TotT = %f / len(tlist) = %d' % (dt, tott, len(tlist)), file=sys.stderr)
### Image
for x in xlist:
for y in ylist:
### To angle.
theta = np.deg2rad(np.sqrt(x * x + y * y))
phi = np.arctan2(y, x)
### To vector
n = self.sc.get_vector(theta, phi)
print(x, y, np.rad2deg(theta), np.rad2deg(phi), n, end=' ', file=sys.stderr)
obsrvr = observer.VirtualObserver(
self.sc,
observer.VirtualSensor(n, self.vena_sc))
sim.set_observer(obsrvr)
# sim.set_tlist(np.arange(0, -20001, -100))
sim.set_tlist(tlist)
sim.enamass = self.mass
f, j = sim.run()
datf.add_data('DATA',
'%f %f %f %f %f %f %f %e %e' %
(x, y, np.rad2deg(theta), np.rad2deg(phi),
n[0], n[1], n[2], f, j))
print('f = %e s^3/m^6' % f, end=' ', file=sys.stderr)
print('j = %e /cm2 sr eV s' % j, file=sys.stderr)
xydone += 1
t1 = time.time()
print('%d/%d %f s ellapsed. Exp left=%fs' % (xydone, xytot, t1 - t0, (t1 - t0) / xydone * (xytot - xydone)), file=sys.stderr)
if output == None:
output = 'main_image_singlepoint_%s.dat' % datetime.datetime.now().strftime('%y%m%d-%H%M%S')
fp = open(output, 'w')
datf.dump(fp)
fp.close()