swim_momintegΒΆ
This is a script to calculate the moment by integration in the LSE coordinate system.
#!/usr/bin/env python
''' This is a script to calculate the moment by integration in the LSE coordinate system.
'''
import matplotlib
from math import atan2,acos,sin
from pylab import *
import os,sys
import datetime
import logging
logging.basicConfig(level=logging.INFO)
logger=logging.getLogger(__file__)
logger.setLevel(logging.DEBUG)
import dateutil
import numpy
import irfpy.swim.swim_start
from irfpy.swim.swim_start import filter
from irfpy.util.vector3d import Vector3d
from irfpy.cy1orb.pvat import getlsevec
na = 21 # Division in azimuth angle
nb = 11 # Division in elevation angle
def ___pstring(str, fp=sys.stderr):
''' A helper function to print.
'''
if fp!=None:
print(str, file=fp)
def _swim_observation_table(t0, fp=None):
''' Returns and print a table of the SWIM observation table.
The input is time. The output is a 2-D double array [N,11].
One line consists of one measurement, and the column corresponds to:
0: time (julday)
1: energy step [0-15]
2: direction step [0-15]
3: count rate [0-]
4: energy in eV
5: velocity in km/s
6-8: velocity in LSE frame (x,y,z) in km/s
9: J(E)
10: f(v)
'''
## get the count data
dat=irfpy.swim.swim_start.getdata(t0, filter=[filter.remove_isolated_one_count, filter.remove_background_5numsum])
fv_=irfpy.swim.swim_start.getvdf(t0, filter=[filter.remove_isolated_one_count, filter.remove_background_5numsum]).getData()
je_=irfpy.swim.swim_start.getdeflux(t0, filter=[filter.remove_isolated_one_count, filter.remove_background_5numsum]).getData()
jd = dat.getJd()
dat = dat.getData()
print(jd, file=sys.stderr)
print(dat, file=sys.stderr)
### For nominal energy and velocity table
etbl=irfpy.swim.SwimEnergy.EnergyTable.table().getTable()
vtbl=irfpy.swim.SwimEnergy.eV2kms(etbl)
print(etbl, file=sys.stderr)
print(vtbl, file=sys.stderr)
### For FOV of SWIM
import irfpy.swim.SwimFov
swimfov = irfpy.swim.SwimFov.SwimFov.getInstance(coords='SC')
### SWIM duty time
import irfpy.swim.SwimTime
dt=irfpy.swim.SwimTime.duty_time()
print('DT=', dt, file=sys.stderr)
### SWIM G-factor
import irfpy.swim.SwimGfact
gfact=irfpy.swim.SwimGfact.gfactor()
gf=gfact.getAbsGProtonStart()[0] ### 0 is for 0.5 keV proton
print('GF@0.5keV=', gf, file=sys.stderr)
for ie in range(16):
for id in range(16):
### Here one can get the vector of the flux in SC.
vsc=swimfov.getDirection(id)
vsc.scale(vtbl[ie])
print(id, vsc, file=sys.stderr)
### The SC vector should be converted to LSE.
vlse = getlsevec(jd, vsc)
print(vlse, file=sys.stderr)
### The J(E) and f(V) is calculated
je = je_[ie][id]
fv = fv_[ie][id]
__pstring('%.5f %2d %2d %6d %4d %7.2f %7.2f %7.2f %7.2f %10.3e %10.3e'%(jd.juld(), ie, id, dat[ie,id], etbl[ie], vtbl[ie], vlse.x, vlse.y, vlse.z, je, fv), fp=sys.stdout)
__solids={}
def swim_solid_fov(id):
''' A first attempt of SWIM solid FoV.
Returns the FoV vectors (unit vectors) of the SWIM in the SC frame correpsonding to id.
SWIM FoV is decided by a shifted normal distribution, but here
we assumed a box with sigma in the cal rep.
'''
global __solids
if id in __solids:
return __solids[id]
cent_A=array([-63.1, -54.78, -46.04, -35.46, -27.27, -17.08, -9.74, -1.32, 7.09, 14.43, 24.62, 32.81, 43.40, 52.13, 62.84, 65.41])
sigma_A=array([7.49, 6.66, 5.74, 4.66, 3.92, 3.20, 2.87, 2.70, 2.77, 3.02, 3.62, 4.28, 5.27, 6.12, 7.09, 7.30])
# #The sigma, s, is exp(-t^2/2s^2). If it is converted to HWHM, multiply 2ln(2)
# hwhm_A = sigma_A* 2*log(2)
cent_B = 0.00
sigma_B = 3.25
# hwhm_B = sigma_B * 2 * log(2)
arr_A = arange(na)
# arr_A = (2 * arr_A / (na-1.) -1 )* hwhm_A[id]
arr_A = (2 * arr_A / (na-1.) -1 )* sigma_A[id]
arr_A = arr_A + cent_A[id]
arr_A = arr_A * pi / 180.
arr_B = arange(nb)
# arr_B = (2 * arr_B / (nb-1.) -1 )* hwhm_B
arr_B = (2 * arr_B / (nb-1.) -1 )* sigma_B
arr_B = arr_B + cent_B
arr_B = arr_B * pi / 180.
### Looking direction in the SC frame is
# xsc = -cos(B) * sin(A)
# ysc = sin(B)
# zsc = cos(B) * cos(A)
fovs=[]
for ia in range(na):
for ib in range(nb):
v=Vector3d( -cos(arr_B[ib])*sin(arr_A[ia]), sin(arr_B[ib]), cos(arr_B[ib])*cos(arr_A[ia]))
fovs.append(v)
__solids[id] = fovs
return fovs
def swim_mominteg(t0, t1):
''' Integrate the ion flux over time
'''
# First, prepare the matrixes of the integration for each energy bin
# Assumption is that the s/c ram velocity is ignored.
### Matrix of the reconvolved distribution function: now fixed to be '360x180' x 16
reconfv = []
for ie in range(16):
aarr=[]
for ia in range(360):
barr=[]
for ib in range(180):
barr.append([])
aarr.append(barr)
reconfv.append(aarr)
### Obtain the observation time
import irfpy.swim.swim_start
obstime=irfpy.swim.swim_start.getobstime(timerange=[t0,t1])
### Rotate the time
for t in obstime:
logger.info('T=%s (started at %s)' % (t, datetime.datetime.now()))
### Get the corresponding data (f(v) in MKSA)
fv_=irfpy.swim.swim_start.getvdf(t, filter=[filter.remove_isolated_one_count, filter.remove_background_5numsum]).getData()
### Rotate the deflection step
for id in range(16):
logger.debug('D=%d/16 started at %s' % (id, datetime.datetime.now()))
### Velocity vector of othe observed particles in the SC frame
fov = swim_solid_fov(id)
### Rotate the energy step
for ie in range(16):
### Now we can specify the direction of the coming ion and the energy step
fv = fv_[ie][id]
for vsc in fov:
# Convert SC to LSE
v=getlsevec(t, vsc)
v.normalize()
# theta is measured from z axis
theta = acos(v.z) * 180./pi
# phi is measured from x axis. z-+x plane is 0, z-+y plane is 90
phi = ( atan2(v.y,v.x) * 180./pi + 360 ) % 360
ia = int(phi)
ib = int(theta)
reconfv[ie][ia][ib].append(fv)
# print ie, ia, ib, reconfv[ie][ia][ib]
# Then, integrate over.
# velocity step and its separation
etbl=irfpy.swim.SwimEnergy.EnergyTable.table().getTable()
vtbl=irfpy.swim.SwimEnergy.eV2kms(etbl) * 1000. ### Order in MKSA
dvtbl=zeros(16)
dvtbl[0]=(vtbl[1]-vtbl[0])/2
for i in range(1,15):
dvtbl[i]=(vtbl[i+1]-vtbl[i-1])/2
dvtbl[15]=(vtbl[15]-vtbl[14])/2
## Initialize the density for three methods.
density_by_max=0.
density_by_ave=0.
density_by_med=0.
fv_avelist = zeros([16, 360, 180]) - 1e30
for ia in range(360):
for ib in range(180):
for ie in range(16):
if len(reconfv[ie][ia][ib]) != 0:
# Averaging f(v)
fvarr = array(reconfv[ie][ia][ib])
fv_ave = fvarr.sum()/len(fvarr)
fv_max = max(fvarr)
fv_med = median(fvarr)
partial_ave = vtbl[ie] * vtbl[ie] * fv_ave * dvtbl[ie] * sin(ib*pi/180.) * (pi/180) * (pi/180)
partial_max = vtbl[ie] * vtbl[ie] * fv_max * dvtbl[ie] * sin(ib*pi/180.) * (pi/180) * (pi/180)
partial_med = vtbl[ie] * vtbl[ie] * fv_med * dvtbl[ie] * sin(ib*pi/180.) * (pi/180) * (pi/180)
# logger.debug('%d %d %d %e %e %e %e' % (ia,ib,ie,fv,partial_ave,partial_max, partial_med))
print('%d %d %d %e %e %e %d' % (ie,ia,ib,fv_ave,fv_max,fv_med,len(fvarr)))
density_by_ave=density_by_ave+partial_ave
density_by_max=density_by_max+partial_max
density_by_med=density_by_med+partial_med
fv_avelist[ie, ia, ib] = fv_ave
print('# Density (fv_ave) = %e'% density_by_ave)
print('# Density (fv_max) = %e'% density_by_max)
print('# Density (fv_med) = %e'% density_by_med)
return (density_by_ave, fv_avelist)
if __name__=='__main__':
''' Main routine to investigate the argument.
'''
import getopt
if len(sys.argv) != 3:
print('***\nUSAGE: %s t0 t1\n***' % (sys.argv[0]))
sys.exit(-1)
t0 = dateutil.parser.parse(sys.argv[1])
t1 = dateutil.parser.parse(sys.argv[2])
logger.info('Start %s' % t0)
logger.info('Stop %s' % t1)
# t0 = datetime.datetime(2009,1,25,13,50,0)
# t1 = datetime.datetime(2009,1,25,14,10,0)
# t0 = datetime.datetime(2009,1,25,14,00,0)
# t1 = datetime.datetime(2009,1,25,14,00,10)
# t0 = datetime.datetime(2009,1,25,14,50,0)
# t1 = datetime.datetime(2009,1,25,15,10,0)
import irfpy.swim.swim_start
obstime=irfpy.swim.swim_start.getobstime(timerange=[t0,t1])
# for t in obstime:
# swim_observation_table(t)
swim_mominteg(t0, t1)