"""Solid angle (with cos theta) of ECLAIRs
"""
from itertools import product
import logging
import numpy as np
from numba import njit
from numba.typed import List
from ecpi.common.instru.model_geom import InstruECLAIRs
from ecpi.common.num.integration_tools import solidangle_rectangle_ana_numba
from ecpi.common.sky.extented_body import EarthInFov
EPS_LIMIT = 2.3e-8
# warnings.simplefilter('ignore', category=NumbaPendingDeprecationWarning)
[docs]def plot(arr_pos): # pragma: no cover
"""..warning: for diagnostic purposes
"""
import matplotlib.pyplot as plt
plt.figure()
for x, y in arr_pos:
plt.scatter(x, y)
plt.show()
[docs]class SolidAngleComputationWithoutEarth(object):
"""Solid angle computation for the ECLAIRs instrument.
"""
def __init__(self, instru_ecl: InstruECLAIRs, use_numba=False):
"""**constructor**
..note :: numba mode only concerns the calculator method
for the moment
:param instru_ecl: description of the mask of ECLAIRs
:type instru_ecl: InstruECLAIRs object
:param use_numba: whether or not to optimize the detector-mask loops
:type use_numba: bool. Default is False.
"""
assert isinstance(instru_ecl, InstruECLAIRs)
self.l_submasks = [instru_ecl._l_mask[m_idx] for m_idx in range(4)]
self.pixpitch = instru_ecl.pitch_pix_size
self.detpix = instru_ecl.p_nb_pix
self.detmas = instru_ecl.dist_mask_detect
self.massize = instru_ecl.mask_total_size
self.masnerv = instru_ecl.half_cross * 2
self.maspsize = instru_ecl.mask_pix_size
self.maspix = instru_ecl.mask_nb_pix
self.detsize = self.pixpitch * self.detpix
self.detpsize = instru_ecl.det_pix_size
self.instru_ecl = instru_ecl
# declare logger
self.logger = logging.getLogger(__name__)
# decide which function to call when numba is (in)active
if use_numba:
self.logger.info('running optimized version')
self.loops_computation_func = _run_over_detector_mask_pixels_numba
else:
self.logger.info('running default version: this may take a while...')
self.loops_computation_func = _run_over_detector_mask_pixels
[docs] def set_mask_configuration(self, sub_mask_arr, fully_opened_mask):
"""Sets the mask configuration: whether or not to compute
solid angle with the real mask geometry or with open submasks.
:param sub_mask_arr: array specifying submask geometry
:type sub_mask_arr: 2D array(int). Filled with 0s and 1s.
:param fully_opened_mask: whether or not to consider the opened mask.
The ECLAIRs maks cross remains !
:type fully_opened_mask: bool
"""
if fully_opened_mask:
self.logger.info(f'Using mask configuration: fully opened submasks')
sub_mask_holes_row, sub_mask_holes_col = np.array(np.where(sub_mask_arr != 2))
else:
self.logger.info(f'Using mask configuration: real mask')
sub_mask_holes_row, sub_mask_holes_col = np.array(np.where(sub_mask_arr == 1))
return sub_mask_holes_row, sub_mask_holes_col
[docs] def run_calculator(self, fully_opened_mask=False, debug=False):
"""Compute solid angle map without the Earth in the FOV.
documentation:
* simulation_correction_cxb_earth_15_05_2020_v2_alg.pdf (version 2) p.10->11
* simulation_correction_pipeline_cxb_earth_15_05_2020_v2_pb.pdf (version 2) p.4
:param fully_opened_mask: whether or not to consider the opened mask.
The ECLAIRs mask cross remains ! Default is False.
:type fully_opened_mask: bool
:return: solid angle map
:rtype: 2D array(float) with shape 80x80
..note: the computation includes the cos theta effect. It follows the sum of the
solid angle is 2199.0922 instead of 2452.0531 (without cos theta case).
"""
solid_angle_with_costheta = np.zeros((self.detpix, self.detpix))
self.logger.info('starting computation of solid angle with Earth')
for sub_mask in self.l_submasks:
sub_mask_arr = sub_mask.ar_val
sub_mask_holes_row, sub_mask_holes_col = self.set_mask_configuration(
sub_mask_arr, fully_opened_mask
)
sub_mask_pos = np.array([sub_mask.idx2pos(np.array([k, l])) \
for (k, l) in zip(sub_mask_holes_row, sub_mask_holes_col)])
self.logger.info(f'processing submask {sub_mask.name}...')
l_detpix = list(product(range(self.detpix), repeat=2))
sa_array = self.loops_computation_func(l_detpix, sub_mask_pos, self.detpix,
self.pixpitch, self.detmas, self.maspsize, self.detpsize)
solid_angle_with_costheta += sa_array
if debug:
self.logger.debug(sub_mask_pos)
plot(sub_mask_pos)
return np.flipud(solid_angle_with_costheta)
def _run_over_detector_mask_pixels(
l_detpix, l_submaskpos, detpix,
pixpitch, detmas, maspsize, detpsize
): # pragma: no cover
"""Loop over detector-mask pixels so as to effectively compute the solid angle.
Pure python implementation.
"""
ox_unit = np.array([1, 0, 0]) # define (Ox) axis normal to detector plane.
sa_array = np.zeros((detpix, detpix))
for (i, j) in l_detpix:
# determine physical center of detector pixel d(i,j)
y_det_c = (i + 0.5) * pixpitch - 0.5 * (detpix * pixpitch)
z_det_c = (j + 0.5) * pixpitch - 0.5 * (detpix * pixpitch)
for (y_msk_c, z_msk_c) in l_submaskpos: # loop over mask holes
# compute solid angle subtended by a pixel m(k,l) of mask as seen by
# the center of detector pixel d(i,j)
yp = y_msk_c - y_det_c
zp = z_msk_c - z_det_c
xp = detmas # + 0.03
# compute solid angle without Earth.
sa_mask_pix = solidangle_rectangle_ana_numba(maspsize, maspsize, yp, zp, xp)
# compute direction from detector pixel to mask pixel.
v_unit = np.array([xp, yp, zp]) / np.linalg.norm([xp, yp, zp])
# compute angles.
cos_theta = np.abs(np.dot(v_unit, ox_unit))
sa_array[i, j] += sa_mask_pix * cos_theta
return sa_array
@njit
def _run_over_detector_mask_pixels_numba(
l_detpix, l_submaskpos, detpix,
pixpitch, detmas, maspsize, detpsize
): # pragma: no cover
"""Loop over detector-mask pixels so as to effectively compute the solid angle.
Python/C implementation (numba).
..warning: seems upgrading to numba==0.48.XXX could cure use of linalg.norm and cos. TBC
"""
sa_array = np.zeros((detpix, detpix))
for (i, j) in l_detpix:
y_det_c = (i + 0.5) * pixpitch - 0.5 * (detpix * pixpitch)
z_det_c = (j + 0.5) * pixpitch - 0.5 * (detpix * pixpitch)
for k in range(len(l_submaskpos)):
y_msk_c = l_submaskpos[k][0]
z_msk_c = l_submaskpos[k][1]
yp = y_msk_c - y_det_c
zp = z_msk_c - z_det_c
xp = detmas
sa_mask_pix = solidangle_rectangle_ana_numba(maspsize, maspsize, yp, zp, xp)
v_unit = np.array([xp, yp, zp]) / (xp ** 2 + yp ** 2 + zp ** 2) ** 0.5
cos_theta = v_unit[0]
sa_array[i, j] += sa_mask_pix * cos_theta
return sa_array
[docs]class SolidAngleComputationWithEarth(SolidAngleComputationWithoutEarth):
"""Solid angle computation for the ECLAIRs instrument
in the presence of Earth in the FOV.
"""
def __init__(self, instru_ecl: InstruECLAIRs, use_numba=False):
"""**constructor**
..note :: numba mode only concerns the calculator method
for the moment
:param instru_ecl: description of the mask of ECLAIRs
:type instru_ecl: InstruECLAIRs object
:param use_numba: whether or not to optimize the detector-mask loops
:type use_numba: bool. Default is False.
"""
super().__init__(instru_ecl)
# decide which function to call when numba is (in)active
if use_numba:
self.logger.info('running optimized version')
self.loops_computation_func = _run_over_detector_mask_pixels_with_earth_numba
else:
self.logger.info('running default version: this may take a while...')
self.loops_computation_func = _run_over_detector_mask_pixels_with_earth
self.use_numba = use_numba
[docs] def run_calculator(self, limb_angle, u_earth, fully_opened_mask=False, debug=False):
"""Compute solid angle map with the Earth in the FOV using
the mean CXB shape calculation.
..note: when angle(vec direction, vec Earth) very close to limb angle
may lead to numerical errors in comparsion. introduced a bias in
solid angle computation. Check has been added that raises a warning.
:param limb_angle: limb angle (cf. doc) in degree
:type limb_angle: float
:param u_earth: unit vector pointing towards Earth center
:type u_earth: [float, float, float]
:param fully_opened_mask: whether or not to consider the opened mask.
The ECLAIRs maks cross remains !
:type fully_opened_mask: bool
:return: solid angle map
:rtype: 2D array(float) with shape 80x80
"""
limb = np.deg2rad(limb_angle)
solid_angle_with_costheta = np.zeros((self.detpix, self.detpix))
self.logger.info('starting computation of solid angle with Earth')
for sub_mask in self.l_submasks:
sub_mask_arr = sub_mask.ar_val
sub_mask_holes_row, sub_mask_holes_col = self.set_mask_configuration(
sub_mask_arr, fully_opened_mask
)
sub_mask_pos = np.array([sub_mask.idx2pos(np.array([k, l]))
for (k, l) in zip(sub_mask_holes_row, sub_mask_holes_col)])
logger = self.logger
self.logger.info(f'processing submask {sub_mask.name}...')
l_detpix = list(product(range(self.detpix), repeat=2))
if self.use_numba:
sa_array = self.loops_computation_func(
l_detpix, sub_mask_pos, self.detpix, self.pixpitch,
self.detmas, self.maspsize, self.detpsize, limb, u_earth
)
else:
sa_array = self.loops_computation_func(
l_detpix, sub_mask_pos, self.detpix, self.pixpitch,
self.detmas, self.maspsize, self.detpsize, limb, u_earth, logger
)
solid_angle_with_costheta += sa_array
if debug:
logger.debug(sub_mask_pos)
plot(sub_mask_pos)
return np.flipud(solid_angle_with_costheta)
def _run_over_detector_mask_pixels_with_earth(
l_detpix, l_submaskpos, detpix,
pixpitch, detmas, maspsize, detpsize,
limb, u_earth, logger
): # pragma: no cover
"""Loop over detector-mask pixels so as to effectively compute the solid angle
in presence of Earth in FOV.
Pure python implementation.
"""
ox_unit = np.array([1, 0, 0]) # define (Ox) axis normal to detector plane.
sa_array = np.zeros((detpix, detpix))
for (i, j) in l_detpix:
# determine physical center of detector pixel d(i,j)
y_det_c = (i + 0.5) * pixpitch - 0.5 * (detpix * pixpitch)
z_det_c = (j + 0.5) * pixpitch - 0.5 * (detpix * pixpitch)
for (y_msk_c, z_msk_c) in l_submaskpos: # loop over mask holes
# compute solid angle subtended by a pixel m(k,l) of mask as seen by
# the center of detector pixel d(i,j)
yp = y_msk_c - y_det_c
zp = z_msk_c - z_det_c
xp = detmas
# compute solid angle without Earth.
sa_mask_pix = solidangle_rectangle_ana_numba(maspsize, maspsize, yp, zp, xp)
sa_array[i, j] += sa_mask_pix
# compute direction from detector pixel to mask pixel.
v_unit = np.array([xp, yp, zp]) / np.linalg.norm([xp, yp, zp])
# compute angles.
cos_alpha = np.dot(v_unit, u_earth)
cos_theta = np.abs(np.dot(v_unit, ox_unit))
alpha = np.arccos(cos_alpha)
if np.abs(limb - alpha) < EPS_LIMIT:
logger.warning(f'Detected pixel with alpha close to limb at pixel ({i}, {j})')
logger.warning(f'limb angle={limb} rad | alpha angle={alpha} rad')
if not alpha < limb:
sa_array[i, j] += sa_mask_pix * cos_theta
return sa_array
@njit
def _run_over_detector_mask_pixels_with_earth_numba(
l_detpix, l_submaskpos, detpix,
pixpitch, detmas, maspsize, detpsize,
limb, u_earth
): # pragma: no cover
"""Loop over detector-mask pixels so as to effectively compute the solid angle
in presence of Earth in FOV.
Python/C implementation (numba).
..warning: seems upgrading to numba==0.48.XXX could cure use of linalg.norm and cos. TBC
"""
# TODO PB: cannot communicate with logger (numba).
sa_array = np.zeros((detpix, detpix))
for (i, j) in l_detpix:
y_det_c = (i + 0.5) * pixpitch - 0.5 * (detpix * pixpitch)
z_det_c = (j + 0.5) * pixpitch - 0.5 * (detpix * pixpitch)
for k in range(len(l_submaskpos)):
y_msk_c = l_submaskpos[k][0]
z_msk_c = l_submaskpos[k][1]
yp = y_msk_c - y_det_c
zp = z_msk_c - z_det_c
xp = detmas
sa_mask_pix = solidangle_rectangle_ana_numba(maspsize, maspsize, yp, zp, xp)
v_unit = np.array([xp, yp, zp]) / (xp ** 2 + yp ** 2 + zp ** 2) ** 0.5
cos_alpha = v_unit[0] * u_earth[0] + v_unit[1] * u_earth[1] + v_unit[2] * u_earth[2] # np.dot(v_unit, u_earth)
cos_theta = v_unit[0]
alpha = np.arccos(cos_alpha)
if not alpha < limb:
sa_array[i, j] += sa_mask_pix * cos_theta
return sa_array
[docs]class SolidAngle(SolidAngleComputationWithoutEarth):
"""Solid angle computation for the ECLAIRs instrument
in the presence of Earth in the FOV.
"""
def __init__(self, instru_ecl: InstruECLAIRs, fully_opened_mask=False):
"""**constructor**
..note :: numba mode only concerns the calculator method
for the moment
:param instru_ecl: description of the mask of ECLAIRs
:type instru_ecl: InstruECLAIRs object
:param use_numba: whether or not to optimize the detector-mask loops
:type use_numba: bool. Default is False.
"""
super().__init__(instru_ecl)
# liste de couple d'indice des pixels detecteur
l_detpix = list(product(range(self.detpix), repeat=2))
self.l_detpix = List()
[self.l_detpix.append(x) for x in l_detpix]
# list des solid angle by hole mask
self.l_sabhm = List()
self.l_sub_mask_pos = List()
self.do_solid_angle_by_hole(fully_opened_mask=fully_opened_mask)
self.no_earth_sa = self.run_calculator(0.0, np.array([-1.0, 0.0, 0.0]))
[docs] def do_solid_angle_by_hole(self, fully_opened_mask=False):
self.logger.info('starting computation of solid angle for each pixel det-mask')
l_sabhm = []
l_sub_mask_pos = []
for sub_mask in self.l_submasks:
sub_mask_arr = sub_mask.ar_val
sub_mask_holes_row, sub_mask_holes_col = self.set_mask_configuration(
sub_mask_arr, fully_opened_mask
)
sub_mask_pos = np.array([sub_mask.idx2pos(np.array([k, l])) \
for (k, l) in zip(sub_mask_holes_row, sub_mask_holes_col)])
if len(sub_mask_pos) == 0:
self.logger.info(f'empty mask {sub_mask.name}')
continue
l_sub_mask_pos.append(sub_mask_pos)
sabhm = solid_angle_for_mask_hole_numba(self.l_detpix,
sub_mask_pos,
self.detpix,
self.pixpitch,
self.detmas,
self.maspsize,
self.detpsize)
l_sabhm.append(sabhm)
[self.l_sabhm.append(x) for x in l_sabhm]
[self.l_sub_mask_pos.append(x) for x in l_sub_mask_pos]
[docs] def run(self, earth_fov):
"""Run and select with or without earth model
:param earth_fov:
"""
assert isinstance(earth_fov, EarthInFov)
if earth_fov.is_in_fov():
return self.run_calculator(earth_fov.limb_angle(), earth_fov.body_xyz())
else:
self.logger.debug('No Earth : return pre-calculated solid angle')
return self.no_earth_sa
[docs] def run_calculator(self, limb_angle, u_earth, fully_opened_mask=False, debug=False):
"""Compute solid angle map with the Earth in the FOV using
the mean CXB shape calculation.
..note: when angle(vec direction, vec Earth) very close to limb angle
may lead to numerical errors in comparsion. introduced a bias in
solid angle computation. Check has been added that raises a warning.
:param limb_angle: limb angle (cf. doc) in degree
:type limb_angle: float
:param u_earth: unit vector pointing towards Earth center
:type u_earth: [float, float, float]
:param fully_opened_mask: whether or not to consider the opened mask.
The ECLAIRs maks cross remains !
:type fully_opened_mask: bool
:return: solid angle map
:rtype: 2D array(float) with shape 80x80
"""
limb = np.deg2rad(limb_angle)
self.logger.info('starting computation of solid angle with Earth')
sa_array = solid_angle_with_earth_numba(self.l_detpix,
self.l_sub_mask_pos,
self.detpix,
self.pixpitch,
self.detmas,
self.maspsize,
self.detpsize,
limb, u_earth,
self.l_sabhm)
# return np.flipud(sa_array)
return sa_array
[docs]class SolidAngleSingleton(SolidAngle):
"""Singleton SolidAngle class, call solid angle by hole one time
"""
instance = None
def __new__(cls, *args, **kargs):
"""
used singleton design pattern
"""
if cls.instance is None:
return object.__new__(cls)
else:
return cls.instance
def __init__(self, instru_ecl: InstruECLAIRs, fully_opened_mask=False):
"""**constructor**
"""
if SolidAngleSingleton.instance is None:
super().__init__(instru_ecl, fully_opened_mask)
SolidAngleSingleton.instance = self
[docs]@njit
def solid_angle_for_mask_hole_numba(
l_detpix, l_submaskpos, detpix,
pixpitch, detmas, maspsize, detpsize
): # pragma: no cover
"""Loop over detector-mask pixels so as to effectively compute the solid angle
in presence of Earth in FOV.
Python/C implementation (numba).
..warning: seems upgrading to numba==0.48.XXX could cure use of linalg.norm and cos. TBC
"""
sa_array = np.zeros((detpix, detpix, len(l_submaskpos)), dtype=np.float32)
for (i, j) in l_detpix:
y_det_c = (i + 0.5) * pixpitch - 0.5 * (detpix * pixpitch)
z_det_c = (j + 0.5) * pixpitch - 0.5 * (detpix * pixpitch)
for k in range(len(l_submaskpos)):
y_msk_c = l_submaskpos[k][0]
z_msk_c = l_submaskpos[k][1]
yp = y_msk_c - y_det_c
zp = z_msk_c - z_det_c
xp = detmas
sa_mask_pix = solidangle_rectangle_ana_numba(maspsize, maspsize, yp, zp, xp)
v_unit = np.array([xp, yp, zp]) / np.sqrt(xp ** 2 + yp ** 2 + zp ** 2)
cos_theta = v_unit[0]
sa_array[i, j, k] = sa_mask_pix * cos_theta
return sa_array
[docs]@njit(cache=True)
def solid_angle_with_earth_numba(
l_detpix, l_submaskpos, detpix,
pixpitch, detmas, maspsize, detpsize,
limb, u_earth, l_sabhm
): # pragma: no cover
"""Loop over detector and all submask.
..note: by Colley JM, optimization of _run_over_detector_mask_pixels_with_earth_numba,
* used precomputing of solid angle of each pixel
det-mask, it s a constant of the system
* include all submask
* earth presence: test on cos angle not angle directly (remove arccos)
* => ratio time ~ 10
"""
sa_array = np.zeros((detpix, detpix))
cos_limb = np.cos(limb)
xp = detmas
xp2 = xp ** 2
for sabhm, submaskpos in zip(l_sabhm, l_submaskpos):
for (i, j) in l_detpix:
y_det_c = (i + 0.5) * pixpitch - 0.5 * (detpix * pixpitch)
z_det_c = (j + 0.5) * pixpitch - 0.5 * (detpix * pixpitch)
for k in range(len(submaskpos)):
y_msk_c = submaskpos[k][0]
z_msk_c = submaskpos[k][1]
yp = y_msk_c - y_det_c
zp = z_msk_c - z_det_c
norm_p = np.sqrt(xp2 + yp ** 2 + zp ** 2)
# scalar product : cos(angle)*norm_p
cos_alpha_m = xp * u_earth[0] + yp * u_earth[1] + zp * u_earth[2]
if cos_alpha_m < cos_limb * norm_p:
sa_array[i, j] += sabhm[i, j, k]
return sa_array