'''SwimSciCount module is an implementation for bridging Sadppac.
SwimSciCount module is an implementation for bridging Sadppac.
Sadppac is python implementation of libsadppac for packet analysis library for SARA.
libsadppac is based on packet structure, therefor, all the available information from the sensor can be obtained.
However, for scientific analysis matrix form is more easier than the packet class.
Here, I would provide a sequence of routines to handle scientific count data.
The source code can be used as well for a sample of using Sadppac module.
>>> print(getBmuUri(939)) # doctest: +SKIP
file///path/to/orb_0939/orb_0939.swim
>>> loader=SwimLoader()
>>> t,x=loader.getStartMatrix(939)
>>> print(x.shape)
(885, 16, 16)
'''
import logging
from Sadppac import *
import datetime
from numpy import *
from irfpy.swim import SwimEnergy
import urllib.request
import urllib.parse
import urllib.error
from irfpy.util.irfpyrc import Rc
import irfpy.util.utc as utc
import irfpy.util.julday as julday
[docs]class SwimLoader:
''' A factory class to load SWIM data.
The data base is open using ${HOME}/.irfpyrc setting.
USAGE:
loader = SwimLoader()
data3000=loader.getStartMatrix(3000)
'''
def __init__(self, rcfile=None):
''' Initialize the loader.
'''
self.__rc=Rc()
if rcfile:
self.__rc.append(rcfile)
self.__baseuri=self.__rc.get('swim', 'bmuobase')
[docs] def getStartMatrix(self, orbnr):
''' Obtain a start counter data matrix of the specified orbit number.
The returne format is same as bmufile2startmatrix. Tuple of (time, startmatx).
@retval time Time is N-element array of datetime instances.
@retval startmatx Matrix in the start counter, which has the NxE(=16)xD(=16) elements.
Here N is the number of data samples.
E and D is the energy and directions.
'''
uri=getBmuUri(orbnr, baseuri=self.__baseuri)
fnam, inf=urllib.request.urlretrieve(uri)
logging.debug('Downloaded: %s %s'%(fnam,inf))
return bmufile2startmatrix(fnam, start=datetime.datetime.min, stop=datetime.datetime.max)
[docs] def getMassMatrix(self, orbnr):
''' Obtain a mass matrix of the specified orbit number.
'''
uri = getBmuUri(orbnr, baseuri=self.__baseuri)
fnam, inf = urllib.request.urlretrieve(uri)
logging.debug('Downloaded: %s %s' % (fnam, inf))
return bmufile2massmatrix(fnam, start=datetime.datetime.min, stop=datetime.datetime.max)
[docs] def getObservationTime(self, orbnr):
''' Get a tuple of observations time in the format of Julday
'''
matx = self.getStartMatrix(orbnr)
t = matx[0]
jd=[]
for t_dt in t:
juld=utc.convert(t_dt, julday.Julday)
jd.append(juld)
return jd
[docs]def getBmuUri(orbnr, baseuri=None):
''' Obtain the uriname of the BMU file of SWIM based on the orbit number.
The base URI is obtained from ${HOME}/.irfpyrc file.
You must have the file and in [swim] section, bmuobase entry have to
specify the correct location.
If the irfpyrc has a problem, RuntimeError will be returned.
'''
if baseuri is None:
rc=Rc()
baseuri=rc.get('swim', 'bmuobase')
if baseuri is None:
raise RuntimeError('The program cannot detect any BMU file name.\nPlease make sure that the file ${HOME}/.irfpyrc exists and readable,\nand bmuobase entry in [swim] section specifies the correct URI.')
uri='%s/orb_%04d/orb_%04d.swim'%(baseuri, orbnr, orbnr)
return uri
[docs]def bmufile2massmatrix(filename, start=datetime.datetime.min, stop=datetime.datetime.max):
logger = logging.getLogger('bmufile2massmatrix')
# logger.setLevel(logging.DEBUG)
reader = SwimBmuReader()
reader.readFile(filename)
sarray = reader.getArray()
ndat = sarray.size()
time_list = []
matrix_list = []
for idx in range(ndat):
try:
logger.debug('Processing %d' % idx)
pac = sarray.getPacket(idx)
utc = pac.getUtc()
dtm = datetime.datetime.utcfromtimestamp(utc)
logger.debug('T = %s' % dtm.strftime('%FT%T'))
if dtm < start or dtm > stop:
continue
try:
pac = pac.decompress()
except:
logger.error('ERROR OF DECOMPRESSION.')
raise
massmatx = pac.getSwimSciAcc()
massmatx = massmatx.reform16x16()
cntr = zeros([16, 16, 32]) # E, D, M
for iE in range(16):
for iD in range(16):
for iM in range(32):
cntr[iE, iD, iM] = massmatx.get(iE, iD, iM)
time_list.append(dtm)
matrix_list.append(cntr)
except:
logger.warn("Skipped by any reasons")
return (time_list, array(matrix_list))
[docs]def bmufile2startmatrix(filename, start=datetime.datetime.min, stop=datetime.datetime.max):
'''The simplest way of obtaining start matrix.
@param[in] filename Filename
@param[in] start Start time of the data to be read.
@param[in] stop Stop time of the data to be read.
@retval time Time is N-element array of datetime instances.
@retval startmatx Matrix in the start counter, which has the NxE(=16)xD(=16) elements.
Here N is the number of data samples.
E and D is the energy and directions.
'''
### First, open the BMU file
reader=SwimBmuReader()
reader.readFile(filename)
sarray=reader.getArray()
ndat=sarray.size()
time_list=[]
matrix_list=[]
for idx in range(ndat):
try:
logging.info('processing %d'%idx)
pac=sarray.getPacket(idx)
utc=pac.getUtc()
dtm = datetime.datetime.utcfromtimestamp(utc)
logging.info(dtm)
if dtm < start or dtm > stop:
logging.info('Skipped data because the time is out of interest.')
continue;
atot=pac.getAccumTotMatx()
atot_16=atot.reform16x16()
time_list.append(dtm)
cntr=zeros([16,16], double) ### 2-dim array. ExD
for ie in range(16):
for id in range(16):
cntr[ie,id]=atot_16.get(ie,id,SwimAccumTotMatx.START_COUNTER)
matrix_list.append(cntr) ### 3-dim array. NxExD
except:
logging.warning('Skipped data because of any reasons.')
matrix_array=array(matrix_list)
return (time_list, matrix_array)
[docs]def simple_vtk_dumper_trial(time_list, matrix_array, filename='swim_simple_vtk_dump.vtk', start=datetime.datetime.min, end=datetime.datetime.max):
'''A very primitive dump into VTK format.
'''
### Velocity list
ene = SwimEnergy.EnergyTable.table().getTable()
vel = SwimEnergy.eV2kms(ene)
### Assumed that data start at x-z plane. No attitude control is assumed.
### Assumed that orbit is exact circle with 119min period.
### So the angular speed is 360[deg]/(119*60[sec]) ~ 0.05 [deg/sec] ~ 8.8e-4 [rad/sec]
###
# The exact information should be obtained by attitude data,
# but the work is not simple. So here, as a simple approach,
# I used 'fixed' parameter.
ang_speed = 360./(119.*60.) * pi/180.
logging.debug('Assumed angular speed of %e'%ang_speed)
### SWIM FoV assumed in X-Z plane.
### Angles assumed 0 to 180 in equally separated.
fovAng = array([0.01, 12, 24, 36, 48, 60, 72, 84, 96, 108, 120, 132, 144, 156, 168, 179.99,], double)
fovRad = fovAng * pi /180.
xsc = -cos(fovRad)
ysc = zeros(16)
zsc = -sin(fovRad)
### First data = no rotation. Second data = rotation compared to nadir (+x).
lvvec=[]
for i in range(len(time_list)):
# for i in range(3):
dt=time_list[i]-time_list[0]
rotAng=(dt.seconds+dt.microseconds/1e6)*ang_speed
### xsc ysc zsc is the original. Rotate rotAng wrt zsc. xyznow is unit vectors for each direction in 'real' space.
xnow = zeros_like(xsc)
ynow = zeros_like(ysc)
znow = zsc
for ix in range(16):
xnow[ix] = xsc[ix] * cos(rotAng) - ysc[ix] * sin(rotAng)
ynow[ix] = xsc[ix] * sin(rotAng) + ysc[ix] * cos(rotAng)
vvec=zeros([16, 16, 3]) # E D xyz
for ie in range(16):
for id in range(16):
vvec[ie,id,0] = xnow[id] * vel[ie]
vvec[ie,id,1] = ynow[id] * vel[ie]
vvec[ie,id,2] = znow[id] * vel[ie]
lvvec.append(vvec)
lvvec=array(lvvec)
### Now lvvec is the matrix of the observation velocities.
### Nelem = T=Nobs * E=16 * D=16 * xyz=3
fp=open(filename, 'w')
print("""# vtk DataFile Version 2.0
SWIM data display - very basic conceptial test by spaghetti code
ASCII
""", file=fp)
### Dump the coordinate system of the dump.
nTime =lvvec.shape[0]
nPoint=lvvec.size/3
print('''DATASET UNSTRUCTURED_GRID
POINTS %d float'''%(nPoint), file=fp)
ted2idx = lambda t,e,d: t*256+e*16+d
for it in range(nTime):
for ie in range(16):
for id in range(16):
print('%f %f %f'%(lvvec[it,ie,id,0],lvvec[it,ie,id,1],lvvec[it,ie,id,2],), file=fp)
# ### First step (experimental) is just using POINT.
# print >> fp, 'CELLS %d %d'%(nPoint, nPoint*2)
# for i in range(nPoint):
# print >> fp, '1 %d'%(i)
# print >> fp, "CELL_TYPES %d"%nPoint
# for i in range(nPoint):
# print >> fp, '2'
### Next challenge is using hexagon.
boxvortex=[]
for it in range(nTime-1):
for ie in range(15):
for id in range(15):
#connection order important???
box=[ted2idx(it,ie,id), ted2idx(it,ie,id+1), ted2idx(it,ie+1,id+1), ted2idx(it,ie+1,id), ted2idx(it+1,ie,id), ted2idx(it+1,ie,id+1), ted2idx(it+1,ie+1,id+1), ted2idx(it+1,ie+1,id)]
boxvortex.append(box)
ncell=len(boxvortex)
print('CELLS %d %d'%(ncell, ncell*9), file=fp)
for i in range(ncell):
print('8 ', end=' ', file=fp)
for j in range(8):
print('%d '%boxvortex[i][j], end=' ', file=fp)
print(file=fp)
print("CELL_TYPES %d"%nPoint, file=fp)
for i in range(nPoint):
print('12', file=fp)
### Define data
print("""POINT_DATA %d
SCALARS log_count_rate float
LOOKUP_TABLE default"""%(nPoint), file=fp)
# for i in range(nPoint):
# print >> fp, i
cnts=log(matrix_array+0.1)
for it in range(nTime):
for ie in range(16):
for id in range(16):
print('%f'%cnts[it,ie,id], file=fp)
fp.close()