r''' G-factor of Vex/IMA
- :func:`gf_af`: G-factor of VEX/IMA in Andrei's definition.
- :func:`gf_mf`: G-factor of VEX/IMA in Markus's definition.
- :func:`gf`: G-factor of ``irfpy`` definition.
**Geometric factor and conversion to differential flux**
A practical use of geometric factor is as follows.
.. math::
J = \frac{C}{G E_c \Delta t}
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 :math:`\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 (:mod:`irfpy.vima.rawdata`).
The simplest way to get the IMA data is using :class:`IMA data center <irfpy.vima.rawdata.DataCenterCount3d>`.
.. 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
.. The count rate can also be obtained from the IMA extra dataset (:mod:`irfpy.vima.imaextra`).
.. warning::
Contents of IMA extra have defined by IRAP, and IRF has no responsibility on that.
Any concerns should be notified to the PI (IRF) before any publication!!
For the case of IMA extra, you can get the (mass separated) count rate as follows.
>>> from irfpy.vima.imaextra import iterextra
>>> for t, c in iterextra(t0, t0 + datetime.timedelta(minutes=3)):
... pass # A little bit tricky here, but just get the data at the time of your preference. Need update.
>>> print(t)
2011-03-26 14:28:10.930000
>>> cH = c['HpSpec']
>>> cO = c['HeavySpec'] # These are the count data. The shape is (P16,E16,A16), a reversed order from the raw data... To be updated.
>>> print(cH.shape)
(16, 96, 16)
>>> print(cH.max())
3283.5754
>>> print(cO.shape)
(16, 96, 16)
>>> print(cO.max())
2.2132316
Note that the order is (P, E, A), the reversed order of the RAW dataset
.. todo::
IMA extra module and functions should follow (MAEPT) order approach.
A new module / functions will be introduced for that purpose.
**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 :func:`gf_irf_H` or :func:`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
'''
import logging
import numpy as np
from irfpy.vima import energy
import irfpy.vima.fov
import irfpy.vima.massring
def __gf_irf_elevation_table():
""" A helper function that is used for :func:`gf_irf_H`.
.. code-author:: Moa Persson
"""
# ===================================================
# Version 3 - corresponds to energy table 3!
# NOTE: [P16, E96]
def _elevation_table(version='3', direction='velocity'):
poltbl = np.zeros((16, 96))
for enidx in range(96):
poltbl[:, enidx] = [irfpy.vima.fov.get_elevation_table(polidx, enidx, version, direction=direction)
for polidx in range(16)]
poltbl[10:, 18:22][poltbl[10:, 18:22] < 0] = np.nan # Change the strange negative values in the table to nan, as I don't know what they mean
return poltbl
elevTable = _elevation_table()
# Elevation angle separation between two neighbouring elevation steps, add the last value twice, approx. correct
delta_elev = np.zeros(elevTable.shape)
delta_elev[:-1, :] = elevTable[:-1, :] - elevTable[1:, :]
delta_elev[-1, :] = delta_elev[-2, :]
# Remove values > 45 deg
idx_above45 = (elevTable > 45.) + (elevTable < -45.)
elevTable[idx_above45] = np.nan
delta_elev[idx_above45] = np.nan
# Change to radians
elevTable_v3 = np.deg2rad(elevTable)
delta_elev_v3 = np.deg2rad(delta_elev)
# ===================================================
# Version 2 - corresponds to energy table 1!
# NOTE: [P16, E96]
def _elevation_table(version='2', direction='velocity'):
poltbl = np.zeros((16, 96))
for enidx in range(96):
poltbl[:, enidx] = [irfpy.vima.fov.get_elevation_table(polidx, enidx, version, direction=direction)
for polidx in range(16)]
poltbl[10:, 18:22][poltbl[10:, 18:22] < 0] = np.nan # Change the strange negative values in the table to nan, as I don't know what they mean
return poltbl
elevTable = _elevation_table()
# Elevation angle separation between two neighbouring elevation steps, add the last value twice, approx. correct
delta_elev = np.zeros(elevTable.shape)
delta_elev[:-1, :] = elevTable[:-1, :] - elevTable[1:, :]
delta_elev[-1, :] = delta_elev[-2, :]
# Remove values > 45 deg
idx_above45 = (elevTable > 45.) + (elevTable < -45.)
elevTable[idx_above45] = np.nan
delta_elev[idx_above45] = np.nan
# Change to radians
elevTable_v2 = np.deg2rad(elevTable)
delta_elev_v2 = np.deg2rad(delta_elev)
return elevTable_v3, delta_elev_v3, elevTable_v2, delta_elev_v2
[docs]def gf_irf_H(prom_index, pacc=6):
""" Calculate the G-factor for proton as a table.
:param prom_index: PROM/EEPROM index. Either 0, 15, or 16. See explanation below.
:return: (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 :class:`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
"""
if pacc != 6:
raise ValueError('Only PACC=6 is supported for VEX/IMA')
gfac_H = np.zeros((16, 96, 16)) # A16, E96, P16
m_per_q = 1
elevTable_v3, delta_elev_v3, elevTable_v2, delta_elev_v2 = __gf_irf_elevation_table()
if prom_index == 0:
etbl = energy.get_default_table_v1()
elevTable = elevTable_v2
elif prom_index in (15, 16):
etbl = energy.get_default_table_v3()
elevTable = elevTable_v3
else:
raise ValueError('Only PROM/EEPROM 0, 15, or 16 is supported')
for azimidx in range(16):
for Eidx, e_per_q in enumerate(etbl):
for elevidx, elevation in enumerate(elevTable[:, Eidx]):
Rm = int(np.floor(irfpy.vima.massring.rm(e_per_q=e_per_q, m_per_q=m_per_q, pi=pacc)))
gfac_H[azimidx, Eidx, elevidx] = gf(e_per_q, pacc, Rm, elevation, azimidx)
return gfac_H
[docs]def gf_irf_O(prom_index, pacc=6):
""" Calculate the G-factor for O+ as a table.
:param prom_index: PROM/EEPROM index. Either 0, 15, or 16. See explanation in :func:`gf_irf_H`.
:return: (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
"""
if pacc != 6:
raise ValueError('Only PACC=6 is supported for VEX/IMA')
gfac_O = np.zeros((16, 96, 16)) # A16, E96, P16
m_per_q = 16
elevTable_v3, delta_elev_v3, elevTable_v2, delta_elev_v2 = __gf_irf_elevation_table()
if prom_index == 0:
etbl = energy.get_default_table_v1()
elevTable = elevTable_v2
elif prom_index in (15, 16):
etbl = energy.get_default_table_v3()
elevTable = elevTable_v3
else:
raise ValueError('Only PROM/EEPROM 0, 15, or 16 is supported')
for azimidx in range(16):
for Eidx, e_per_q in enumerate(etbl):
for elevidx, elevation in enumerate(elevTable[:, Eidx]):
Rm = int(np.floor(irfpy.vima.massring.rm(e_per_q=e_per_q, m_per_q=m_per_q, pi=pacc)))
gfac_O[azimidx, Eidx, elevidx] = gf(e_per_q, pacc, Rm, elevation, azimidx)
return gfac_O
class _matlab_array:
""" Emulate matlab array, starting from index 1
Only index access functions are overrided.
>>> python_list = ['one', 'two', 'three']
>>> matlab_list = _matlab_array(python_list)
>>> print(len(matlab_list))
3
>>> print(matlab_list[1])
one
>>> matlab_list[2] = 'tvaa'
>>> print(matlab_list[2])
tvaa
>>> print(python_list[1])
two
"""
def __init__(self, python_list):
self._list = list(python_list)
def __getitem__(self, index):
if index <= 0:
raise IndexError('Matlab support index from 1')
return self._list[index - 1]
def __setitem__(self, index, value):
if index <= 0:
raise IndexError('Matlab support index from 1')
self._list[index - 1] = value
def __len__(self):
return len(self._list)
class _Constants:
dGf2Gfe = 0.292 # ; % GFe (return value) = dGf(center)/dPhi * dGf2Gfe, see VexImaCalLog
Esatur = [0.0, 300.0, 600.0] # ;
Agfe = [-0.58600, -0.58600, -0.58600] # ;
Bgfe = [-5.19938, -5.19938, -5.19938] # ;
CorrGFRM = [-1.0,-1.0,-1.0,7.0,5.0,3.6,2.0,1.2,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0, 1.0,1.3,1.6,1.9,2.2,2.8,3.5,4.0,5.0] # ;
ElevP = [ -46.7, -39.3, -32.3, -25.2, -17.8, -7.6, 2.2, 11.7, 21.4, 28.5, 34.9, 41.3, 46.9] # ;
ElevCorr = [ 0.038, 0.468, 1.182, 0.902, 0.924, 0.983, 0.999, 0.960, 0.963, 1.070, 0.976, 0.904, 0.371] # ;
SectCorr = [ 1.027, 1.198, 1.339, 1.428, 1.452, 1.407, 1.300, 1.148, 0.973, 0.802, 0.661, 0.572, 0.548, 0.593, 0.700, 0.852] # ;
[docs]def gf_af(Energy, pac, Rm, Elevation, Sector, dGf2Gfe=None, Esatur=None, Agfe=None, Bgfe=None, CorrGFRM=None, ElevP=None, ElevCorr=None, SectCorr=None):
""" Andrei's gfactor.
:param energy: Energy in the unit of ``eV``
:param pac: PAC index. Usually, it should be 6, but depending on the PAC setting. (0, 3, 6) is used for VEX.
:param rm: Mass ring number. 0 to 31.
:param elevation: Elevation angle.
:param sector: Sector index. 0 to 15.
:keyword dgf2gfe: Not sure what this parameter is, but use default.
:keyword esatur: Not sure what this parameter is, but use default.
:keyword agfe: Not sure what this parameter is, but use default.
:keyword bgfe: Not sure what this parameter is, but use default.
:keyword corrGFRm: Not sure what this parameter is, but use default.
:keyword elevP: Not sure what this parameter is, but use default.
:keyword elevCorr: Not sure what this parameter is, but use default.
:keyword 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 :file:`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
"""
# ;=============================================================
# ; Fill the common constants at the first time
# ;=============================================================
# if n_elements(Esatur) LT 1 then begin
# print,'VexGf.pro V.2.0, 10 Feb 2008'
# dGf2Gfe = 0.292 ; GFe (return value) = dGf(center)/dPhi * dGf2Gfe, see VexImaCalLog
# Esatur = [0.0, 300.0, 600.0]
# Agfe = [-0.58600, -0.58600, -0.58600]
# Bgfe = [-5.19938, -5.19938, -5.19938]
# CorrGFRM = [-1.0,-1.0,-1.0,7.0,5.0,3.6,2.0,1.2,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0, 1.0,1.3,1.6,1.9,2.2,2.8,3.5,4.0,5.0]
# ElevP = [ -46.7, -39.3, -32.3, -25.2, -17.8, -7.6, 2.2, 11.7, 21.4, 28.5, 34.9, 41.3, 46.9]
# ElevCorr = [ 0.038, 0.468, 1.182, 0.902, 0.924, 0.983, 0.999, 0.960, 0.963, 1.070, 0.976, 0.904, 0.371]
# SectCorr = [ 1.027, 1.198, 1.339, 1.428, 1.452, 1.407, 1.300, 1.148, 0.973, 0.802, 0.661, 0.572, 0.548, 0.593, 0.700, 0.852]
# endif
#
if dGf2Gfe is None: dGf2Gfe = _Constants.dGf2Gfe
if Esatur is None: Esatur = _Constants.Esatur
if Agfe is None: Agfe = _Constants.Agfe
if Bgfe is None: Bgfe = _Constants.Bgfe
if CorrGFRM is None: CorrGFRM = _Constants.CorrGFRM
if ElevP is None: ElevP = _Constants.ElevP
if ElevCorr is None: ElevCorr = _Constants.ElevCorr
if SectCorr is None: SectCorr = _Constants.SectCorr
logger = logging.getLogger('gf_af')
# logger.setLevel(logging.DEBUG)
## logger.debug('Rm={}'.format(Rm))
# ;============================================================
# ; Check position first of all
# ;============================================================
# if Rm LT 3.0 then return, -1.0
if Rm < 3.0:
return -1.0
#
# ;=============================================================
# ; Calculation of original GF value
# ;=============================================================
# case PAC of
# 0: ipac = 0
# 3: ipac = 1
# 6: ipac = 2
# else: ipac = 2
# endcase
if pac == 0:
ipac = 0
elif pac == 3:
ipac = 1
elif pac == 6:
ipac = 2
else:
ipac = 2
# if Energy LT Esatur(ipac) then Gfl = exp(Agfe(ipac) * alog(Esatur(ipac)) + Bgfe(ipac)) else $
# Gfl = exp(Agfe(ipac) * alog(Energy) + Bgfe(ipac))
if Energy < Esatur[ipac]:
Gfl = np.exp(Agfe[ipac] * np.log(Esatur[ipac]) + Bgfe[ipac])
## logger.debug('Gfl(if0)={}'.format(Gfl))
else:
Gfl = np.exp(Agfe[ipac] * np.log(Energy) + Bgfe[ipac])
## logger.debug('Gfl(if1)={}'.format(Gfl))
#
# ;=============================================================
# ; Position correction
# ;=============================================================
# firm = 2.99
firm = 2.99
# while (firm LT 31.0) AND (firm LT Rm) do firm += 1.0
while firm < 31.0 and firm < Rm:
firm += 1.0
## logger.debug('firm={}'.format(firm))
# irm = fix(firm+0.2) < 31
irm = int(np.floor(firm + 0.2))
if irm > 31:
irm = 31
## logger.debug('irm={}'.format(irm))
# RmFactor = (CorrGFRM(irm) - CorrGFRM(irm-1))*(Rm - float(irm-1)) + CorrGFRM(irm-1)
RmFactor = (CorrGFRM[irm] - CorrGFRM[irm - 1]) * (Rm - (irm - 1.0)) + CorrGFRM[irm - 1]
## logger.debug('RmFactor={}'.format(RmFactor))
# Gfl = GFl/RmFactor
Gfl = Gfl / RmFactor
#
# ;=============================================================
# ; Elevation correction
# ;=============================================================
# Nel = n_elements(ElevP)
Nel = len(ElevP)
# if Elevation LT ElevP(1) then ie = 1 else if Elevation GT ElevP(Nel-2) then ie = Nel-1 else begin
if Elevation < ElevP[1]:
ie = 1
elif Elevation > ElevP[Nel - 2]:
ie = Nel - 1
else:
# ie = 1
# while (ie LT Nel-1) AND (Elevation GT ElevP(ie)) do ie++
ie = 1
while ie < Nel - 1 and Elevation > ElevP[ie]:
ie += 1
# endelse
## logger.debug('ie={}'.format(ie))
# ElevFactor = (ElevCorr(ie) - ElevCorr(ie-1))/(ElevP(ie) - ElevP(ie-1))*(Elevation - ElevP(ie-1)) + ElevCorr(ie-1)
ElevFactor = (ElevCorr[ie] - ElevCorr[ie - 1]) / (ElevP[ie] - ElevP[ie - 1]) * (Elevation - ElevP[ie - 1]) + ElevCorr[ie - 1]
## logger.debug('ElevFactor={}'.format(ElevFactor))
# Gfl = GFl*ElevFactor
Gfl = Gfl * ElevFactor
#
# ;=============================================================
# ; Azimuth correction
# ;=============================================================
# Gfl = GFl*SectCorr(Sector)
Gfl = Gfl * SectCorr[Sector]
#
# return, GFl * dGf2Gfe
return Gfl * dGf2Gfe
[docs]def gf_af_matlab(Energy, pac, Rm0, Elevation, Sector0, dGf2Gfe=None, Esatur=None, Agfe=None, Bgfe=None, CorrGFRM=None, ElevP=None, ElevCorr=None, SectCorr=None):
""" Andrei's gfactor
.. note::
This function is out of date.
Keep for debug purpose.
Use :func:`gf_af`.
:param energy: Energy in the unit of ``eV``
:param pac: PAC index. Usually, it should be 6, but depending on the PAC setting. (0, 3, 6) is used for VEX.
:param rm0: Mass ring number. 0 to 31.
:param elevation: Elevation angle.
:param sector0: Sector index. 0 to 15.
:keyword dgf2gfe: Not sure what this parameter is, but use default.
:keyword esatur: Not sure what this parameter is, but use default.
:keyword agfe: Not sure what this parameter is, but use default.
:keyword bgfe: Not sure what this parameter is, but use default.
:keyword corrGFRm: Not sure what this parameter is, but use default.
:keyword elevP: Not sure what this parameter is, but use default.
:keyword elevCorr: Not sure what this parameter is, but use default.
:keyword 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 :func:`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 :file:`gfacVEX_AF.m`.
.. todo::
Vectorize, or make it as an easy class (to make the calculation faster)
"""
Rm = Rm0 + 1 # Rm is 1 based
# Elevation = Elevation0 + 1 # Elevation input looks angle in degrees.
Sector = Sector0 + 1
if dGf2Gfe is None: dGf2Gfe = _Constants.dGf2Gfe
if Esatur is None: Esatur = _Constants.Esatur
Esatur = _matlab_array(Esatur)
if Agfe is None: Agfe = _Constants.Agfe
Agfe = _matlab_array(Agfe)
if Bgfe is None: Bgfe = _Constants.Bgfe
Bgfe = _matlab_array(Bgfe)
if CorrGFRM is None: CorrGFRM = _Constants.CorrGFRM
CorrGFRM = _matlab_array(CorrGFRM)
if ElevP is None: ElevP = _Constants.ElevP
ElevP = _matlab_array(ElevP)
if ElevCorr is None: ElevCorr = _Constants.ElevCorr
ElevCorr = _matlab_array(ElevCorr)
if SectCorr is None: SectCorr = _Constants.SectCorr
SectCorr = _matlab_array(SectCorr)
#### % Check position first of all
#### if (Rm <4.0)
#### result= -1.0;
#### return
#### end
if Rm < 4.0:
return -1.0
####
#### % Calculation of original GF value
####
#### switch PAC
#### case 0
#### ipac = 1;
#### case 3
#### ipac = 2;
#### case 6
#### ipac = 3;
#### otherwise
#### ipac = 3;
#### end
if pac == 0:
ipac = 1
elif pac == 3:
ipac = 2
elif pac == 6:
ipac = 3
else:
ipac = 3
####
#### if Energy <Esatur(ipac)
#### Gfl = exp(Agfe(ipac) * log(Esatur(ipac)) + Bgfe(ipac));
#### else
#### Gfl = exp(Agfe(ipac) * log(Energy) + Bgfe(ipac));
#### end
if Energy < Esatur[ipac]:
Gfl = np.exp(Agfe[ipac] * np.log(Esatur[ipac]) + Bgfe[ipac])
else:
Gfl = np.exp(Agfe[ipac] * np.log(Energy) + Bgfe[ipac])
####
#### %disp(strcat('Gfl=',mat2str(Gfl)));
####
#### % Position correction
#### firm = 3.99 ;
firm = 3.99
#### while (firm <32.0) && (firm <Rm)
#### firm = firm+1.0;
#### end;
while firm < 32 and firm < Rm:
firm = firm + 1
#### irm = floor(firm+0.2); % Should be maximum 32
irm = np.floor(firm + 0.2)
#### if (irm>32)
#### irm=32;
#### end
if irm > 32:
irm = 32
irm = int(irm)
#### RmFactor = double((CorrGFRM(irm) - CorrGFRM(irm-1))*(Rm - double(irm-1)) + CorrGFRM(irm-1));
RmFactor = float((CorrGFRM[irm] - CorrGFRM[irm - 1]) * (Rm - float(irm - 1)) + CorrGFRM[irm - 1] )
#### Gfl = Gfl/RmFactor;
Gfl = Gfl / RmFactor
#### % Elevation correction
#### Nel = length(ElevP);
#### if Elevation <ElevP(2)
#### ie = 1; # It must be ie=2, (YF, 161106)
#### elseif Elevation> ElevP(Nel-1)
#### ie = Nel;
#### else
#### ie = 1; # It must be ie=2, (YF, 161106)
#### while (ie <Nel) && (Elevation> ElevP(ie))
#### ie=ie+1;
#### end
#### end
#### ElevFactor = double(ElevCorr(ie) - ElevCorr(ie-1))/(ElevP(ie) - ElevP(ie-1))*(Elevation - ElevP(ie-1)) + ElevCorr(ie-1);
# ElevFactor = _elevation_factor_yf(Elevation) # Intention of AF code is just an interpolation.
ElevFactor = _elevation_factor_af(Elevation) # Comparison confirmed between AF impl and YF impl, assuming the
# translation error (?) in ie=1...
# _yf uses simpler way (using scipy interpolation) but takes more exec time
# As long as _af version is "verified", I use _af version here.
#### Gfl = Gfl*ElevFactor;
Gfl = Gfl * ElevFactor
#### % Azimuth correction
#### Gfl = Gfl*double(SectCorr(Sector));
Gfl = Gfl * float(SectCorr[Sector])
#### result=Gfl * dGf2Gfe;
result = Gfl * dGf2Gfe
#### end
return result
def _elevation_factor_af(Elevation):
""" (Examination code. Do not use!!)
This is an examination code to understand what is elevation factor
and how it is implemented.
Comparison with _elevation_factor_yf was done (_script_elevcheck function).
Conclusion: The implementation by _elevation_factor_yf is more readable;
short, clear, and intuitive.
On the other hand, af version is faster.
"""
ElevP = None
if ElevP is None: ElevP = _Constants.ElevP
ElevP = _matlab_array(ElevP)
ElevCorr = None
if ElevCorr is None: ElevCorr = _Constants.ElevCorr
ElevCorr = _matlab_array(ElevCorr)
#### % Elevation correction
#### Nel = length(ElevP);
Nel = len(ElevP)
#### if Elevation <ElevP(2)
#### ie = 1;
if Elevation < ElevP[2]:
#ie = 1 ###YF: THIS WILL CAUSE PROBLEM. COULD BE ie=2??
ie = 2
#### elseif Elevation> ElevP(Nel-1)
#### ie = Nel;
elif Elevation > ElevP[Nel - 1]:
ie = Nel
#### else
#### ie = 1;
#### while (ie <Nel) && (Elevation> ElevP(ie))
#### ie=ie+1;
#### end
#### end
else:
ie = 1
while ie < Nel and Elevation > ElevP[ie]:
ie = ie + 1
#### ElevFactor = double(ElevCorr(ie) - ElevCorr(ie-1))/(ElevP(ie) - ElevP(ie-1))*(Elevation - ElevP(ie-1)) + ElevCorr(ie-1);
ElevFactor = float(ElevCorr[ie] - ElevCorr[ie -1]) / (ElevP[ie] - ElevP[ie -1]) * (Elevation - ElevP[ie -1]) + ElevCorr[ie - 1]
return ElevFactor
def _elevation_factor_yf(Elevation):
""" (Examination code. Do not use)
What is the elevation factor?
Conclusion: It is just a linear interpolation.
"""
ElevP = _Constants.ElevP
ElevCorr = _Constants.ElevCorr
from scipy import interpolate
interpfunc = interpolate.interp1d(ElevP, ElevCorr, fill_value='extrapolate')
return interpfunc(Elevation)
def _script_elevcheck():
""" (Examination code. Do not use)
Plot the af version and yf version of elevation factor.
"""
factdic = {}
factdic2 = {}
for elev in range(-60, 61, 2): # Every 2 deg
try:
efact = _elevation_factor_af(elev)
except IndexError:
efact = 1e-3
try:
efact2 = _elevation_factor_yf(elev)
except IndexError:
efact2 = 1e-3
# print(elev, efact)
factdic[elev] = efact
factdic2[elev] = efact2
for p, f in zip(_Constants.ElevP, _Constants.ElevCorr):
factdic[p] = f
factdic2[p] = f
# print(factdic)
keys = [p for p in factdic]
keys.sort()
x = []
y = []
for p in keys:
print(p, factdic[p])
x.append(p)
y.append(factdic[p])
import matplotlib.pyplot as plt
plt.plot(x, y, 'x')
keys = [p for p in factdic2]
keys.sort()
x = []
y = []
for p in keys:
print(p, factdic2[p])
x.append(p)
y.append(factdic2[p])
import matplotlib.pyplot as plt
plt.plot(x, y, 'o')
plt.yscale('log')
plt.yscale('log')
plt.savefig('elevaf.png')
[docs]def gfactor_matrix(verbose=False):
""" 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 :meth:`gf_af` 3x16x96x16x32 = 2359296 times.
.. note::
(Developing memo) the code is translated from Matlab ``gfacVEX.m``.
"""
### %%% MATLABized version of Andrei's vexgf.pro, V.2.0
### %%% Creates G-factor matrix gfac(paccindex,elevation,energy,sector,Rm)
### %%% Output of gfacVEX_AF is 'Gfe'=energy-geometrical factor for one sector
### %%% Note: Gfe=0.292*dGF
### %%% See also:
### %%% http://goezog.cesr.fr:8080/DD_VEX/BIBLES/IMA/CALIBRREP/VexImaCalLog.html
###
###
### % These constants should be read only once.
### disp('Runs gfacVEX...')
### disp(' To apply to Rest matrix: MATLABized version of VexGf.pro V.2.0, 10 Feb 2008');
### dGf2Gfe = 0.292; % GFe (return value) = dGf(center)/dPhi * dGf2Gfe, see VexImaCalLog
### Esatur = [0.0, 300.0, 600.0];
### Agfe = [-0.58600, -0.58600, -0.58600];
### Bgfe = [-5.19938, -5.19938, -5.19938];
### CorrGFRM = [-1.0,-1.0,-1.0,7.0,5.0,3.6,2.0,1.2,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0, 1.0,1.3,1.6,1.9,2.2,2.8,3.5,4.0,5.0];
### ElevP = [ -46.7, -39.3, -32.3, -25.2, -17.8, -7.6, 2.2, 11.7, 21.4, 28.5, 34.9, 41.3, 46.9];
### ElevCorr = [ 0.038, 0.468, 1.182, 0.902, 0.924, 0.983, 0.999, 0.960, 0.963, 1.070, 0.976, 0.904, 0.371];
### SectCorr = [ 1.027, 1.198, 1.339, 1.428, 1.452, 1.407, 1.300, 1.148, 0.973, 0.802, 0.661, 0.572, 0.548, 0.593, 0.700, 0.852];
###
###
### pacc=[0 3 6];
pacc = [6, 3, 0]
###
### % Allocating:
### ALLgfac=zeros(3,16,96,16,32);
ALLgfac = np.zeros([3, 16, 96, 16, 32])
###
### for pci=1:3 % Pacc index
### for elev=0:15
### for ee=1:96
### for sec=0:15
### for mch=1:32
ene_table = energy.get_default_table_v1()
logger = logging.getLogger('gfactor_matrix')
if verbose:
logger.setLevel(logging.INFO)
for pci in range(3):
for elev in range(16): # 0-based index
for ee in range(96): # 0-based index
for sec in range(16): # 0-based index. Argument should be really the angle though.
for mch in range(32): # 0-based index.
### ALLgfac(pci,elev+1,ee,sec+1,mch)=gfacVEX_AF(E(ee), ...
### pacc(pci),mch,elev+1,sec+1,dGf2Gfe,Esatur,Agfe, ...
### Bgfe,CorrGFRM,ElevP,ElevCorr,SectCorr);
ALLgfac[pci, elev, ee, sec, mch] = gf_af(ene_table[ee], pacc[pci], mch, (elev - 7.5) * 5.625, sec)
logger.info('{} {} {} {} {} {}'.format(ALLgfac[pci, elev, ee, sec, mch], pci, elev, ee, sec, mch))
### end
### end
### end
### end
### end
return ALLgfac
[docs]def gf_mf(E, m):
""" Markus version of g-factor
.. warning::
It is not recommended to use this g-factor.
Use :func:`gf_af`.
:param E: Energy in eV.
:param 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.
"""
### %%% This is a MATLABized version of Markus
### %%% Fränz' function 'asp_ima_gfac.pro' (only VEX-part, reduced for work
### %%% with IMAEXTRA data which have only two different products (HpSpec and
### %%% HeavySpec).
###
### disp(' To apply to HPSPEC and HEAVYSPEC: scheme adopted from M Fränz asp_ima_gfac.pro');
### gfac_arr=zeros(96,3);
### species=['HPSPEC'; 'HEAVYS'];
###
### for spi=1:2
### for pci=1:3
G0 = [5.5e-5, 7.63944e-5][m]
mq = [1.0, 16.0][m]
G1 = 6.09935e-5 -1.43219e-6*mq;
B = np.log(G1 / G0) / np.log(10.0 / 3.0)
A = np.log(G0) - B * np.log(300.0)
if E < 100.0: # Orignal code says 100 eV, but 300 eV makes more sense, indeed.
gfac = E * 0.0 + G0
else:
gfac = np.exp(A + B * np.log(E))
return gfac
### % Formulae for VEX from Fedorov email to M Fränz 20071009
### if strcmp(species(spi,1:6),'HPSPEC')
### G0 = 5.5e-05; mq=1;
### else
### G0 = 7.63944e-05; mq=16.0;
### end
### G1 = 6.09935e-05 -1.43219e-06*mq;
### B = log(G1/G0)/log(10.0/3.0);
### A = log(G0) - B * log(300.0);
### highin=find(E > 100.0);
### gfac=E*0.0+G0;
### gfac(highin)= exp(A + B * log(E(highin)));
###
### gfac_arr(:,pci)=gfac;
### end
### if spi==1
### HpSpecGfac=gfac_arr;
### else
### HeavySpecGfac=gfac_arr;
### end
### end]]
###
###
###
gf = gf_af
""" G-factor to be used in general.
This is an alias to :func:`gf_af`.
"""