irfpy.vima.gfactor
¶
G-factor of Vex/IMA
gf_af()
: G-factor of VEX/IMA in Andrei’s definition.gf_mf()
: G-factor of VEX/IMA in Markus’s definition.gf()
: G-factor ofirfpy
definition.
Geometric factor and conversion to differential flux
A practical use of geometric factor is as follows.
Here J is the differential flux (in /cm2 s sr eV), G is the geometric factor (in cm2 sr eV/eV), Ec is the central energy, and \(\Delta t\) is the duty time. C is the count of each corresponding bin. The geometric factor here includes the efficiency and per bin.
Count rate
The count rate can also be obtained from the raw data (irfpy.vima.rawdata
).
The simplest way to get the IMA data is using IMA data center
.
Note
However, one should be careful because the dataset is not calibrated. For example, one should pre-process the data, such as
background subtraction
mass composition separations
saturation correction
before getting the meaningful results. The processing depends on the purpose of analysis. These issues are not discussed here.
>>> import datetime
>>> t0 = datetime.datetime(2011, 3, 26, 14, 28)
>>> from irfpy.vima import rawdata
>>> dc = rawdata.DataCenterCount3d() # An instance of data center
>>> tdata, datavalue = dc.nextof(t0) # Get the data next of ``t0``. Data is stored to ``datavalue``.
>>> print(tdata)
2011-03-26 14:28:10.931640
>>> cnts = datavalue.matrix
>>> print(cnts.shape) # The data is sotred to an array with a shape of (M32, A16, E96, P16)
(32, 16, 96, 16)
>>> cnts = cnts.sum(0) # Mass collapsed.
>>> print(cnts.max()) # Check the maximum count in a bin
2227.0
>>> print(cnts.sum()) # Check the total count in a bin
60516.0
In this case, we assumed all the counts are from hproton. >>> cH = cnts
Central energy (Ec)
The energy table depends on EEPROM section, which can be retrieved as follows.
>>> promsection = datavalue.hk.promsection # PROM section, which defines the energy table
>>> print(promsection) # For this particular time, PROMSECTION=16(RAW) was used
16
Then, the energy table is known.
>>> etable = energy.get_default_table(16)
Duty time
Duty time for IMA is constantly 120.9 ms.
>>> dt = 0.1209
Geometric factor
For a simple use of the geometric factor is using functions gf_irf_H()
or gf_irf_O()
.
These functions return tables, corresponding to the EEPROM as shown above.
>>> gH0 = gf_irf_H(0)
>>> gH15 = gf_irf_H(15)
>>> gO0 = gf_irf_O(0)
>>> gO15 = gf_irf_O(15)
Any of the above table is (A16, E96, P16) shaped array. Note the order of the dimension.
The unit is in cm2 sr eV/eV
.
Flux
Thus, the differential flux is calculated for each bin.
>>> j = cH / (gH15 * etable[np.newaxis, :, np.newaxis] * dt) # Transpose for IMA extra dataset!!
>>> print(j.shape) # A16, E96, P16
(16, 96, 16)
>>> print('{:.3e}'.format(np.nanmax(j)))
4.070e+06
>>> print('{:.3e}'.format(np.nansum(j))) # NOTE: THIS OPERATION IS NOT MEANINGFUL FOR ANY PHYSICS!!
6.909e+07
- irfpy.vima.gfactor.gf_irf_H(prom_index, pacc=6)[source]¶
Calculate the G-factor for proton as a table.
- Parameters
prom_index – PROM/EEPROM index. Either 0, 15, or 16. See explanation below.
- Returns
(16, 96, 16) shaped array containing G-factor for protons. The unit is cm2 sr eV/eV, as a usual sense. Order of the array is (A=16, E=96, P=16) order. This is the same order as raw data matrix, while the reversed order as IMA EXTRA data.
>>> gf_v1 = gf_irf_H(0) # The first half of the mission. >>> gf_v3 = gf_irf_H(15) # The last half of the mission.
The argument is the PROM index. See the later part of this document.
>>> print(gf_v1.shape) (16, 96, 16) >>> print('{:.3e}'.format(gf_v1[10, 8, 5])) # For A=10, E=80, P=5, for the first table 3.713e-06 >>> print('{:.3e}'.format(gf_v3[10, 8, 5])) # For A=10, E=80, P=5, for the second table 4.843e-06
PROM index defines the energy and elevation tables in the instrument. It is embedded in HK stream, and the
rawdata <irfpy.imacommon.imascipac.CntMatrix
contain the information. An example as below.>>> from irfpy.vima import rawdata >>> import datetime >>> dc = rawdata.DataCenterCount3d() # An instance of data center >>> tdata, datavalue = dc.nextof(datetime.datetime(2011, 3, 26, 14, 28)) # Get the data next of ``t0`` >>> print(tdata) # Time slightly shifted from IMA extra, but should be ok... 2011-03-26 14:28:10.931640 >>> prom_index = datavalue.hk.promsection # PROM section, which defines the energy table >>> print(prom_index) # For this particular time, PROM/EEPROM index was 16(RAW). 16
- irfpy.vima.gfactor.gf_irf_O(prom_index, pacc=6)[source]¶
Calculate the G-factor for O+ as a table.
- Parameters
prom_index – PROM/EEPROM index. Either 0, 15, or 16. See explanation in
gf_irf_H()
.- Returns
(16, 96, 16) shaped array containing G-factor for oxygen ion. The unit is cm2 sr eV/eV, as a usual sense. Order of the array is (A=16, E=96, P=16) order. This is the same order as raw data matrix, while the reversed order as IMA EXTRA data.
>>> gfO_v1 = gf_irf_O(0) # The first half of the mission. >>> gfO_v3 = gf_irf_O(15) # The last half of the mission.
The argument is the PROM index. See the later part of this document.
>>> print(gfO_v1.shape) (16, 96, 16) >>> print('{:.3e}'.format(gfO_v1[10, 8, 5])) # For A=10, E=80, P=5, for the first table 7.426e-07 >>> print('{:.3e}'.format(gfO_v3[10, 8, 5])) # For A=10, E=80, P=5, for the second table 9.686e-07
- irfpy.vima.gfactor.gf_af(Energy, pac, Rm, Elevation, Sector, dGf2Gfe=None, Esatur=None, Agfe=None, Bgfe=None, CorrGFRM=None, ElevP=None, ElevCorr=None, SectCorr=None)[source]¶
Andrei’s gfactor.
- Parameters
energy – Energy in the unit of
eV
pac – PAC index. Usually, it should be 6, but depending on the PAC setting. (0, 3, 6) is used for VEX.
rm – Mass ring number. 0 to 31.
elevation – Elevation angle.
sector – Sector index. 0 to 15.
dgf2gfe – Not sure what this parameter is, but use default.
esatur – Not sure what this parameter is, but use default.
agfe – Not sure what this parameter is, but use default.
bgfe – Not sure what this parameter is, but use default.
corrGFRm – Not sure what this parameter is, but use default.
elevP – Not sure what this parameter is, but use default.
elevCorr – Not sure what this parameter is, but use default.
sectCorr – Not sure what this parameter is, but use default.
- Returns
G-factor in the unit of
cm^2 sr eV/eV
.
Note
The code is translated from
vexgf.pro
V3.0>>> print('{:.4e}'.format(gf_af(100, 6, 30.5, 0, 0))) 8.6235e-06 >>> print('{:.4e}'.format(gf_af(300, 6, 30.6, 0, 0))) 8.4360e-06
- irfpy.vima.gfactor.gf_af_matlab(Energy, pac, Rm0, Elevation, Sector0, dGf2Gfe=None, Esatur=None, Agfe=None, Bgfe=None, CorrGFRM=None, ElevP=None, ElevCorr=None, SectCorr=None)[source]¶
Andrei’s gfactor
Note
This function is out of date. Keep for debug purpose. Use
gf_af()
.- Parameters
energy – Energy in the unit of
eV
pac – PAC index. Usually, it should be 6, but depending on the PAC setting. (0, 3, 6) is used for VEX.
rm0 – Mass ring number. 0 to 31.
elevation – Elevation angle.
sector0 – Sector index. 0 to 15.
dgf2gfe – Not sure what this parameter is, but use default.
esatur – Not sure what this parameter is, but use default.
agfe – Not sure what this parameter is, but use default.
bgfe – Not sure what this parameter is, but use default.
corrGFRm – Not sure what this parameter is, but use default.
elevP – Not sure what this parameter is, but use default.
elevCorr – Not sure what this parameter is, but use default.
sectCorr – Not sure what this parameter is, but use default.
Note
The elevation looks angle in degree. Confirmed by AF.
Note
The inputs of this function is described in l1tol2 document. Or, see
gf_mf()
that calls this function.Note
To convert to differential flux from count rate will be given soon, but the source of information would be in l1tol2 document.
Note
The code is translated from
gfacVEX_AF.m
.Todo
Vectorize, or make it as an easy class (to make the calculation faster)
- irfpy.vima.gfactor.gfactor_matrix(verbose=False)[source]¶
Obtain G-factor matrix.
Obtain the g-factor matrix.
- Returns
(3, 16, 96, 16, 32) shaped array with g-factor information.
Todo
Write more. What is g-factor matrix. How to interpret etc.
Note
This function takes extremely long execution time. Well, it is natural. It calls
gf_af()
3x16x96x16x32 = 2359296 times.Note
(Developing memo) the code is translated from Matlab
gfacVEX.m
.
- irfpy.vima.gfactor.gf_mf(E, m)[source]¶
Markus version of g-factor
Warning
It is not recommended to use this g-factor. Use
gf_af()
.- Parameters
E – Energy in eV.
m – 0 for proton, 1 for oxygen
The g-factor has no dependency on directions and mass-bins (and even no dependency on PAC!).
Note
(Developing memo) the code is translated from Matlab
gfacVEX.m
.Note
The original code
gfacVEX.m
claims 100 eV is a “turn-around” energy, but it really looks like 300 eV, then continuous curves was found.Todo
Describe the findings. ;)
Todo
Vecotorize the function.