''' CENA mass matrix with time-axis information.
:mod:`cena_mass2` module is the packet base library.
Now the :mod:`cena_mass_time` module is prepared to
implement the time series data of the mass mode matrix.
'''
import datetime
import numpy as np
import logging
from irfpy.cena import cena_mass2 as cm
[docs]def getfullHdeflux(t0, t1, validation=None, mainhv=2750):
''' Get the time series of the flux between ``t0`` and ``t1``.
:param t0: Time start (inclusive)
:type t0: ``datetime.datetime``
:param t1: Time end (inclusive)
:type t1: ``datetime.datetime``.
:keyword mainhv: Main HV value should be more than this value. Otherwise, ignored.
:returns: A set of (*obstime*, *matrix of the differential flux of H in the unit of #/cm2 eV sr s*).
The differential flux matrix has the shape of (16, 7, N) where N is the number of observation time.
See also :meth:`getfulldata` and :meth:`irfpy.cena.cena_mass2.getHdefluxE16`.
>>> t0 = datetime.datetime(2009, 7, 13, 11, 17)
>>> t1 = datetime.datetime(2009, 7, 13, 11, 17, 10)
>>> timelist, fulldata = getfullHdeflux(t0, t1)
>>> fulldata.shape
(16, 7, 3)
Simple way of averaging.
>>> aveflux = fulldata.mean(2)
>>> print(aveflux[0, 3])
--
>>> print('%.2f' % aveflux[8, 3])
59527.07
>>> print('%.2f %.2f %.2f' % (fulldata[8, 3, 0], fulldata[8, 3, 1], fulldata[8, 3, 2]))
48703.97 64938.62 64938.62
'''
logger = logging.getLogger('getfullHdefluxE16')
tlist = cm.getobstime(timerange=(t0, t1))
tlist_valid = []
flux_valid = []
for t in tlist:
try:
mhv = cm.getHvMain(t).getData()
if mhv < mainhv:
logger.warn('HV not enough at %s. Ignored.' % t)
continue
flx = cm.getHdefluxE16(t, validation=validation)
tlist_valid.append(t)
flux_valid.append(flx.getData())
except RuntimeError as e:
logger.warn('No or invalid data at %s. Skipped' % t)
continue
if len(tlist_valid) == 0:
raise RuntimeError('No data between the time range [%s, %s]' %
(t0, t1))
flux_valid = np.ma.array(flux_valid)
flux_valid = flux_valid.transpose((1, 2, 0))
return tlist_valid, flux_valid
[docs]def getfulldata(t0, t1, validation=None, mainhv=2750):
''' Get the data obtained between ``t0`` and ``t1``.
:param t0: Time start (inclusive)
:type t0: ``datetime.datetime``
:param t1: Time end (inclusive)
:type t1: ``datetime.datetime``.
:keyword mainhv: Main HV value should be more than this value. Otherwise, ignored.
:returns: A set of (*obstime*, *matrix of the counts*). Counts are just raw value.
But if validation is specified, invalid data is masked.
:rtype: A set of (list of ``datetime`` with shape of (N, ), ``numpy.ma.ndarray`` with shape of (16, 7, 128, N)),
where N is the number of observation.
If the data is not found during the time period,
``RuntimeError`` will be raised.
>>> t0 = datetime.datetime(2009, 7, 13, 11, 17)
>>> t1 = datetime.datetime(2009, 7, 13, 11, 17, 10)
>>> timelist, fulldata = getfulldata(t0, t1)
>>> fulldata.shape
(16, 7, 128, 3)
>>> print(len(fulldata.compressed())) # Number of valid data
21504
>>> print(fulldata.sum()) # Total counts. May change after DB changes.
374.0
The corresponding data was obtained by SV table 2, energy range 0-3 and 12-16 is
masked. Check this fact.
>>> print(fulldata[0, 4, 77, 2])
--
>>> len(fulldata[:4, :, :, :].compressed())
0
>>> print(fulldata[14, 6, 27, 0])
--
>>> len(fulldata[12:, :, :, :].compressed())
0
Hydrogen count can be get via summing mass channel 1 to 63.
>>> hfulldata = fulldata[:, :, 1:64, :]
>>> print(hfulldata.sum()) # H total counts. May change after DB cahnges.
346.0
Examine the validation.
>>> timelist2, fulldata2 = getfulldata(t0, t1, validation=[cm.invalid_mask.high_heavy_mask,
... cm.invalid_mask.high_count_mask])
>>> timelist == timelist2
True
>>> fulldata2.shape
(16, 7, 128, 3)
>>> len(fulldata2.compressed()) # Valid bins, smaller than fulldata.
15360
>>> fulldata2.sum() # Total count valid. Different from fulldata.sum()
268.0
Examine the HV mask.
>>> t0 = datetime.datetime(2009, 6, 13, 8, 0, 0)
>>> t1 = datetime.datetime(2009, 6, 13, 8, 15, 0)
>>> tlist3, fdata3 = getfulldata(t0, t1, mainhv=0) # All the data retrieval
>>> print(len(tlist3))
209
>>> print(fdata3.shape)
(16, 7, 128, 209)
>>> tlist4, fdata4 = getfulldata(t0, t1, mainhv=2750) # Only HV up data.
>>> print(len(tlist4))
183
>>> print(fdata4.shape)
(16, 7, 128, 183)
'''
logger = logging.getLogger('getfulldata')
tlist = cm.getobstime(timerange=(t0, t1))
tlist_valid = []
count_valid = []
for t in tlist:
try:
mhv = cm.getHvMain(t).getData()
if mhv < mainhv:
logger.warn('HV not enough at %s. Ignored.' % t)
continue
data = cm.getdataE16(t, validation=validation)
count_valid.append(data.getData())
tlist_valid.append(t)
except RuntimeError as e:
logger.warn('No or invalid data on %s. Ignored' % t)
continue
if len(tlist_valid) == 0:
raise RuntimeError('No data between the time range [%s, %s]' %
(t0, t1))
count_valid = np.ma.array(count_valid)
count_valid = count_valid.transpose((1, 2, 3, 0))
return tlist_valid, count_valid