'''swim_start package is a collection for analyzing SWIM start counter.
SWIM start counter is the most outer counter for SWIM.
But since start counter is placed after the energy separation,
the data is really useful for scientific analysis, while
no mass analysis is available.
In this module, softwares providing SWIM start counter are implemented.
Importing swim_start package, and use getdata(), getdeflux(), getvdf() functions.
SwimEnergy, SwimFov, and SwimTime modules are also useful.
Problems:
swim_start package has yet a problem.
Main problem is memory leak.
Multiple swim_start.getdata() causes a terrible memory leak.
Most probably the problem occurs in Sadppac package or the interface with Python.
I haven't specified suspective point in the code...
The work is progressed in the branch of memory_leak_swimstart.
'''
from Sadppac import *
import datetime
import numpy
from irfpy.util.utc import *
from irfpy.util.julday import *
from irfpy.util.fivenumsum import *
from irfpy.cy1orb.Cy1OrbitNr import Cy1OrbitNr
from irfpy.swim.SwimSciCount import *
from irfpy.swim.flux import *
import logging
from irfpy.swim import SwimTime
from irfpy.swim import SwimGfact
from irfpy.swim import SwimEnergy
from irfpy.swim.swim_cache import SwimDataCache
[docs]def getobstime(orbit=None, timerange=None):
''' Obtain the observation time of SWIM for specified orbit (primary) or specified time range (secondary).
@se SwimDataCache##getObstime() method.
'''
obst=__swimcache.cache().getObstime(orbit=orbit, timerange=timerange)
return obst
[docs]def getdata(t, filter=[]):
''' Obtain SWIM start count rate data of the specified time.
If no data is available, RuntimeError is returned.
@param t Time
@param filter A list of callable to filter the data matrix.
@retval JdObject instance with E=16xD=16 numpy.array object.
Filter:
Some useful filters are implemented in filter class.
You can use it like
getdata(datetime(2009,1,25,15,0,0), filter=[filter.remove_one_count, ])
Customized filter functions can be made. The callable have to have an argument of (16,16) shaped numpy-array, and returns (16,16) shaped numpy-array.
'''
# orbnr=Cy1OrbitNr()
# orbnr.setFromDefaultUri()
# print __swimcache.cache()
time_dt = convert(t, datetime.datetime)
matx_jdo=__swimcache.cache().getStartMatrix(time_dt)
matx_jd=matx_jdo.getJd() # Julian day
matx=matx_jdo.getData() # numpy.array 16x16
for filt in filter:
matx=filt(matx)
return JdObject(matx_jd, matx)
[docs]class filter:
''' Collection of filters to be set to the count rate.
'''
info='' # Saved information.
[docs] @classmethod
def remove_one_count(self, raw):
''' All the bin with count '1' is set to be '0'. Dark count can provide occasionally count=1.
'''
return where(raw <= 1, 0, raw)
[docs] @classmethod
def remove_isolated_one_count(self, raw):
''' A bin with count '1' is set to be '0' if the four-neighbors of the bin are 0.
This is more preferable way than self.remove_one_count, but takes computer time.
'''
shape=raw.shape
if len(shape)!=2 or shape[0]!=16 or shape[1]!=16:
logging.warn('Filter remove_ioslated_one_count only support 16x16 matrix')
return numpy.array(raw).copy()
m = zeros([18,18]) # m is a 18x18 matrix in order to treat 'edge'.
m[1:17,1:17]=raw.copy() # 16x16 raw matrix is inserted to m[1..16,1..16]. m[0 or 17, 0 or 17] is 0.
for ie in range(1, 17):
for id in range(1, 17):
if m[ie,id]<=1 and m[ie+1,id]==0 and m[ie-1,id]==0 and m[ie,id+1]==0 and m[ie,id-1]==0:
m[ie,id]=0
val=m[1:17,1:17].copy()
logging.debug(val)
return val
[docs] @classmethod
def remove_background_5numsum(self, raw, factor=1.):
''' Remove a background using 5 value summary.
Median is considered as a background level.
Maximum inner fence is considered as a threshold of the background level.
'''
# print 'Removing background!'
shape=raw.shape
if len(shape)!=2 or shape[0]!=16 or shape[1]!=16:
logging.warn('Filter remove_ibackground_5numsum only support 16x16 matrix')
return numpy.array(raw).copy()
rsum=numpy.array(raw).sum()
m=numpy.array(raw).copy()
m=m.reshape(m.size)
five = fivenumsum(m)
bg=five[0] # Background is median
# bg=five[1] # Background is lower fence
bgfence = five[4] # Background threshold is maximum inside.
bgambig=bgfence-bg
bgambig = bgambig*float(factor)
bgthr=bg+bgambig # Maximum of inner fence is the threshold.
logging.debug('Background level = %f'%bg)
logging.debug('Background threshold = %f'%bgthr)
m=m.reshape(16,16)
for ix in range(16):
for iy in range(16):
if m[ix][iy] <= bgthr:
m[ix][iy]=0;
else:
m[ix][iy]=m[ix][iy]-bg
rsum2=numpy.array(raw).sum()
assert rsum==rsum2
self.info='BG=%f BGth=%f'%(bg,bgthr)
return m
[docs]def getdeflux1cnt(strategy="draft_simple_gfactor", etblVersion=2):
''' Returns the one count level for SWIM start counter.
@retval One count level in E=16xD=16 numpy array.
'''
if strategy == 'draft_simple_gfactor':
#Count rate assuming 1.
cnt=ones([16,16])
geav = SwimGfact.getSimpleAbsG0()
grel = SwimGfact.gfactor().getRelatG()
ge = [ geav * grel[i] for i in range(16)]
dt=SwimTime.duty_time()
e = SwimEnergy.EnergyTable.table(version=etblVersion).getTable()
flx=[]
for ie in range(16):
f=[]
for id in range(16):
f = [ simple_c2j(cnt[ie][id], dt, ge[id], 1, e[ie]) for id in range(16) ]
flx.append(f)
return numpy.array(flx)
[docs]def getdeflux(time_dt, strategy="draft_simple_gfactor", etblVersion=2, filter=[]):
''' Obtain differential flux from SWIM start counter.
@param strategy You can specify the strategy to be used. Presently only "draft_simple_gfactor" is available and is the default.
@param filter You can specify a list of 'callable' in this argument. This is called after getting the raw data to pre-process. Refer to self.getdata() method for details.
@retval Differentian energy flux in the unit of [# / cm^2 sr s eV] for E=16xD=16 matrix.
@exception If no data is available, RuntimeError is returned as self.getdata() returns.
@exception If the obtained data is not in maximum mode(16x16), RuntimeError is returned.
Strategy:
"draft_simple_gfactor":
Using draft values of g-factor.
Relative efficiencies between looking directions are supported.
But relative difference of g-factor in energy bin directions are not supported.
It uses G0(avg) = 5e-5 [cm^2 sr eV/eV] for absolute value (SwimGfact.getSimpleAbsG0()).
Anycase, G0(avg) is in the errors of G0 in the draft cal rep issued in July 2009.
'''
if strategy == 'draft_simple_gfactor':
# Get the count rate date.
jdobj = getdata(time_dt, filter=filter) # 16x16 tuple, otherwise, RuntimeError!
jd=jdobj.getJd()
cnt=jdobj.getData()
if len(cnt) != 16 or len(cnt[0]) != 16:
raise RuntimeError('The obtained data is not in the maximum mode (E=16/D=16). E=%d D=%d'%(len(cnt), len(cnt[0])))
# Parameter
geav = SwimGfact.getSimpleAbsG0()
grel = SwimGfact.gfactor().getRelatG()
ge = [ geav * grel[i] for i in range(16)]
# print ge
dt = SwimTime.duty_time()
e = SwimEnergy.EnergyTable.table(version=etblVersion).getTable()
flx=[]
for ie in range(16):
f=[]
for id in range(16):
f = [ simple_c2j(cnt[ie][id], dt, ge[id], 1, e[ie]) for id in range(16) ]
flx.append(f)
return JdObject(jd, numpy.array(flx))
else:
logging.error('No stragety %s supported'%strategy)
raise RuntimeError('No strategy %s supported.'%strategy)
[docs]def getvdf1cnt(strategy="draft_simple_gfactor", etblVersion=2):
''' Returns the one count level for SWIM start counter.
@param strategy You can specify the strategy to be used. Presently only "draft_simple_gfactor" is available and is the default. Refer to self.getdeflux() method for details.
'''
deflux = getdeflux1cnt(strategy=strategy, etblVersion=etblVersion)
e=SwimEnergy.EnergyTable.table(version=etblVersion).getTable()
vdf=[]
for ie in range(16):
vdfe = [simple_j2f(deflux[ie][id], e[ie]) for id in range(16)]
vdf.append(vdfe)
return numpy.array(vdf)
[docs]def getvdf(time_dt, strategy="draft_simple_gfactor", etblVersion=2, filter=[]):
''' Returns velocity distribution function from SWIM start counter.
@param strategy You can specify the strategy to be used. Presently only "draft_simple_gfactor" is available and is the default. Refer to self.getdeflux() method for details.
@param filter You can specify a list of 'callable' in this argument. This is called after getting the raw data to pre-process. Refer to self.getdata() method for details.
@retval Differentian energy flux in the unit of [ / m^6 s^3] for E=16xD=16 matrix.
@exception If no data is available, RuntimeError is returned as self.getdata() returns.
@exception If the obtained data is not in maximum mode(16x16), RuntimeError is returned.
'''
# Get the differential flux
deflux = getdeflux(time_dt, strategy=strategy, etblVersion=etblVersion, filter=filter)
jd = deflux.getJd()
J = numpy.array(deflux.getData()) # J is the differential flux in / cm^2 sr eV s.
# Energy table
e = SwimEnergy.EnergyTable.table(version=etblVersion).getTable()
### Calculation
vdf=[]
# One can assume that the result is E=16xD=16 matrix.
for ie in range(16):
vdfe = [ simple_j2f( J[ie][id], e[ie] ) for id in range(16)] # VDF in MKSA.
vdf.append(vdfe)
return JdObject(jd, numpy.array(vdf))
class __swimcache:
instance = None
@classmethod
def cache(cls):
if cls.instance is None:
cls.instance = SwimDataCache()
return cls.instance
[docs]def clearCache():
logging.debug('Clearing cache')
__swimcache.cache().clearCache()
[docs]def statusCache():
return __swimcache.cache().statusCache()