"""
@author: BACON Philippe, APC/IN2P3/CNRS
"""
import logging
import numpy as np
from astropy.io import fits
import astropy.units as u
import astropy.table as at
from astropy.constants import R_earth
import transforms3d.quaternions as tq
from scipy.linalg import norm
import ecpi.common as ec
from ecpi.common.instru.model_geom import InstruECLAIRs
from ecpi.common.mission.attitude import xyz_to_direlev, direlev_to_xyz
from ecpi.common.mission.attitude import AttitudeECLAIRs
logger = logging.getLogger(__name__)
FILES = {
'orbit': ec.add_path_data_ref_eclairs('instru/SVO-ORB-CNV-20190121T125959_MXT.0.fits'),
'attitude': ec.add_path_data_ref_eclairs('instru/SVO-ATT-CNV-20190121T125959_MXT.0.fits')
}
[docs]def earth_limb(position):
"""
Compute alpha limb angle
Limb angle is defined to be the angle between the earth direction and
a tangent line (see AlG slides for scheme and maths)
documentation:
* simulation_correction_cxb_earth_15_05_2020_v2_alg.pdf (version 2) p.5
* simulation_correction_pipeline_cxb_earth_15_05_2020_v2_pb.pdf (version 2) p.3
:param position: instrument position (X, Y, Z) in geocentric frame.
Assume coordinates are given in km. Default is None.
:type position: array
:return: limb angle in degrees
:rtype: float
"""
d_earth_sat = np.linalg.norm(position)
logger.info(f'distance Earth center - sat: {d_earth_sat} km')
arg = (1e-3 * R_earth.value) / (d_earth_sat)
return np.rad2deg(np.arcsin(arg))
[docs]class EclairsOrbit(object):
"""Manage orbital aspects for the SVOM mission.
"""
def __init__(self,
orbit_filename=None,
attitude_filename=None):
"""
**Constructor**
:param orbit_filename: name of the orbit .fits file
:type orbit_filename: str
:param attitude_filename: name of the attitude .fits file
:type attitude_filename: str
"""
init_orb = False
init_att = False
if orbit_filename is not None:
assert isinstance(orbit_filename, str)
self.orbit_filename = orbit_filename
init_orb = True
if attitude_filename is not None:
assert isinstance(attitude_filename, str)
self.attitude_filename = attitude_filename
init_att = True
self._init_logger()
if init_orb and init_att:
self._load()
def _init_logger(self):
self.logger = logging.getLogger(__name__)
self.logger.info('creating an instance of EclairsOrbit')
def _load(self):
"""Read SVOM orbital SVO-ORB-CNV and SVO-ATT-CNV attitude files.
.. warning: some attributes defined below are likely to be unused in the end. TBC
"""
with fits.open(self.orbit_filename) as orf_:
orbital_data = orf_[1].data
self.logger.info(f'reading orbit file {self.orbit_filename} ...')
with fits.open(self.attitude_filename) as atf_:
attitude_data = atf_[1].data
self.logger.info(f'reading attitude file {self.attitude_filename} ...')
assert len(orbital_data)==len(attitude_data)
self.nrows = len(orbital_data)
self.metadata = {
'time_pvt': 'time from PVT packet',
'orbit_id': 'orbital identifier',
'position': '(X,Y,Z) position vector in m in J2000',
'velocity': '(VX,VY,VZ) velocity vector in deg/s',
'time': 'mean time AAV PVT used to compute last angles',
'elv': 'angle between pointing and earth limb',
'qparam': 'attitude quaternion'
}
self.names = self.metadata.keys()
self.dtypes = ('f8', 'U20', np.ndarray, np.ndarray, 'f8', 'f8', np.ndarray)
self.data = at.Table(
[
orbital_data['time_pvt'],
orbital_data['orbit_id'],
orbital_data['position'],
orbital_data['velocity'],
orbital_data['time'],
orbital_data['elv'],
attitude_data['qparam']
],
names=self.names,
meta=self.metadata,
dtype=self.dtypes
)
self.data['time_pvt'].unit = u.s
self.data['orbit_id'].unit = u.dimensionless_unscaled
self.data['position'].unit = u.km
self.data['velocity'].unit = u.km/u.s
self.data['time'].unit = u.s
self.data['elv'].unit = u.deg
self.data['qparam'].unit = u.dimensionless_unscaled
[docs] def get_field(self, field):
"""
Get field column from already loaded data.
:param field: name of the required field
:type field: str
"""
if not field in self.names:
self.logger.error(f'{field} field not understood.')
return None
return self.data[field]
def _set_field(self, position, quater, idx_row):
"""Set quaternion and position of n-th row in astropy table.
For tests purposes: likely to be removed in the end.
.. warning: to be removed in the end.
:param position: instrument position (X, Y, Z) in geocentric frame.
Assume coordinates are given in km. Default is None.
:type position: array
:param quater: quaternion (q0, q1, q2, q3). Default is None.
:type quater: array
:param idx_row: row index in astropy table
:type idx_row: int
"""
self.data['qparam'][idx_row] = quater
self.data['position'][idx_row] = position
# TODO: all methods below have duplicates in class EarthInFOV
[docs] def earth_position(self, position=None, quater=None, idx_row=None):
"""Return Earth angular position (theta_E, phi_E)
in field of view.
See section 9.1 (p.50/65) in SV-SY-NT-242-JPO document (v4).
See mail/code from Sujay 17 Sept 2019
.. note :: transformation (dir, elev) <-> (theta, phi) is
theta = 90 - elev and phi = dir.
:param position: instrument position (X, Y, Z) in geocentric frame.
Assume coordinates are given in km. Default is None.
:type position: array
:param quater: quaternion (q0, q1, q2, q3). Default is None.
:type quater: array
:param n: row index in astropy table. Default is None.
:type n: int
:return: unit vector pointing to the center of the Earth in FOV frame
:rtype: array(float) of shape (3,)
"""
if idx_row is not None:
assert idx_row<=self.nrows
position = self.data['position'][idx_row]
quater = self.data['qparam'][idx_row]
else:
assert position is not None
assert quater is not None
self.logger.info(f'[Xsat={position[0]}, Ysat={position[1]}, Zsat={position[2]}]')
if not isinstance(position, np.ndarray):
sat_gcrs = np.array(position)
else:
sat_gcrs = position
mag = norm(sat_gcrs)
if not mag > R_earth.value * 1e-3:
msg = f"read position is inside Earth ({sat_gcrs})"
self.logger.error(msg)
raise ValueError(msg)
else:
self.logger.info(f'sat_gcrs: {sat_gcrs}')
# account for ECLAIRS/VT misalignment
att = AttitudeECLAIRs()
att.set_attitude_instru_quater(quater)
quater = att._q_j2000_instf
# Earth centre dirn in ECL frame is reverse of sat pos vec
earth_gcrs = -sat_gcrs/mag
earth_sat = tq.rotate_vector(earth_gcrs, tq.qconjugate(quater))
# convert from sat gcrs position to (theta, phi)
direlev = xyz_to_direlev(earth_sat[np.newaxis, :], deg=True)[0]
self.logger.info(f'dir Earth={direlev[0]} deg | elev earth={direlev[1]} deg')
# add reflection
self.logger.warning('added reflection in Y (cf. AlG/UGTS)')
print('added reflection in Y (cf. AlG/UGTS)')
#TODO: array must have ecpi convention, see ticket https://forge.in2p3.fr/issues/41460
earth_sat[1] = -earth_sat[1]
return earth_sat
[docs] def limb_angle(self, position=None, idx_row=None):
"""Compute alpha limb angle at different times so as to
determine whether a given pixel falls on the Earth surface or not.
Limb angle is defined to be the angle between the earth direction and
a tangent line (see AlG slides for scheme and maths)
documentation:
* simulation_correction_cxb_earth_15_05_2020_v2_alg.pdf (version 2) p.5
* simulation_correction_pipeline_cxb_earth_15_05_2020_v2_pb.pdf (version 2) p.3
:param position: instrument position (X, Y, Z) in geocentric frame.
Assume coordinates are given in km. Default is None.
:type position: array
:param idx_row: row index in astropy table. Default is None.
:type idx_row: int
:return: limb angle in degrees
:rtype: float
"""
if idx_row is not None:
assert idx_row<=self.nrows
position = self.data['position']#[n]
else:
assert position is not None
d_earth_sat = norm(position)
if not d_earth_sat > R_earth.value * 1e-3:
msg = f"read position is inside Earth ({d_earth_sat})"
self.logger.error(msg)
raise ValueError(msg)
else:
self.logger.info(f'sat_gcrs: {d_earth_sat}')
return earth_limb(position)
[docs] def is_earth_in_fov(self, position=None, quater=None, idx_row=None):
"""Check whether Earth appears in the instrument FOV.
documentation:
* simulation_correction_cxb_earth_15_05_2020_v2_alg.pdf (version 2) p.6
* simulation_correction_pipeline_cxb_earth_15_05_2020_v2_pb.pdf (version 2) p.3
:param position: instrument position (X, Y, Z) in geocentric frame.
Assume coordinates are given in km. Default is None.
:type position: array
:param quater: quaternion (q0, q1, q2, q3). Default is None.
:type quater: array
:param idx_row: row index in astropy table. Default is None.
:type idx_row: int
:return: if Earth is in FOV
:rtype: bool
"""
instru_ecl = InstruECLAIRs()
if idx_row is not None:
assert idx_row<=self.nrows
position = self.data['position'][idx_row]
quater = self.data['qparam'][idx_row]
else:
assert position is not None
assert quater is not None
mag = norm(position)
if not mag > R_earth.value * 1e-3:
msg = f"read position is inside Earth ({mag})"
self.logger.error(msg)
raise ValueError(msg)
else:
self.logger.info(f'sat_gcrs: {mag}')
size_mask_dpix = instru_ecl.mask_total_size + instru_ecl.p_nb_pix*instru_ecl.pitch_pix_size
arg = np.sqrt(2) * size_mask_dpix / (2 * instru_ecl.dist_mask_detect)
omega_half = np.rad2deg(np.arctan(arg))
self.logger.info(f'omega half: {omega_half} deg')
limb_angle = self.limb_angle(position)
u_earth = self.earth_position(position, quater)
theta_earth = np.rad2deg(np.arccos(u_earth[0]))
print('theta:', theta_earth)
cond = not (np.abs(theta_earth) > limb_angle + omega_half)
self.logger.info(f'Earth in FOV: {cond}')
return cond
[docs] def open_fov_map(self, position=None, quater=None, idx_row=None):
"""Create the open field of view map, ie. a detector size array
(199x199) filled with 0 and 1. Each pixel value is 0 if Earth is
present in the LoS and 1 otherwise.
docstring TBC
:param idx_row: row index in astropy table
:type idx_row: int
:return: open fov 199x199 map
:rtype: 2D array(bool)
"""
size_fov_map = 199
instru_ecl = InstruECLAIRs()
open_fov_map = np.empty(
shape=(size_fov_map, size_fov_map),
dtype=np.int16
)
if idx_row is not None:
assert idx_row<=self.nrows
limb = self.limb_angle(idx_row=idx_row)
vec_earth = self.earth_position(idx_row=idx_row)
else:
assert position is not None
assert quater is not None
limb = self.limb_angle(position)
vec_earth = self.earth_position(position, quater)
direlev_earth = xyz_to_direlev(vec_earth, deg=True)
dir_earth = direlev_earth[0,0]
elev_earth = direlev_earth[0,1]
print(f'limb angle: {limb} deg ({np.deg2rad(limb)} rad)')
print(f'u_earth: {vec_earth}')
print(f'dir_earth={dir_earth} deg | elev_earth={elev_earth} deg')
mag = norm(position)
if not mag > R_earth.value * 1e-3:
msg = f"read position is inside Earth ({mag})"
self.logger.error(msg)
raise ValueError(msg)
else:
self.logger.info(f'sat_gcrs: {mag}')
limb = np.deg2rad(limb)
for idx_y in range(size_fov_map):
for idx_z in range(size_fov_map):
# compute corresponding direction of a detector pixel
elev_ij, dir_ij = instru_ecl.skypix_to_elevdir(
(size_fov_map//2)-idx_y, (size_fov_map//2)-idx_z
)
#TODO: +180 ?? voir ticket https://forge.in2p3.fr/issues/41460
ptg_ij = direlev_to_xyz(
np.array([dir_ij+180, elev_ij]),
deg=True)
# check condition and fill open fov map
# no need for normalization as direlev_to_xyz returns unit vector.
theta_s = np.arccos(np.dot(ptg_ij, vec_earth))
bool_cond = theta_s > limb
open_fov_map[idx_y,idx_z] = bool_cond
#TODO: np.flipud ?? voir ticket https://forge.in2p3.fr/issues/41460
return np.flipud(open_fov_map)
[docs] def compute_earth_frac(self, open_fov_map):
"""Compute earth fraction in FOV;
:param open_fov_map: sky map with directions occulted by Earth
:type open_fov_map: array(bool)
:return: earth fraction
:rtype: float
"""
size_fov_map = 199
frac = np.round(np.sum(open_fov_map) / size_fov_map**2, 3)
print(f'fraction of fov covered by Earth: {1-frac}')
return 1-frac
[docs] def is_in_saa(self):
"""Check whether the satellite is in the South Atlantic Anomaly (SAA).
.. note :: likely to be provided by a .fits file. No need to be recomputed.
"""
raise NotImplementedError