'''
Created on August 26 2019
author: Philippe Bacon, 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
import transforms3d.utils as tu
from ecpi.common.mission.attitude import AttitudeVT, AttitudeECLAIRs
from ecpi.common.mission.attitude import xyz_to_thetaphi, direlev_to_xyz
from ecpi.common.instru.ref_geom import ECLAIRsGeom
from ecpi.common.instru.model_geom import InstruECLAIRs
[docs]class EclairsOrbit(object):
def __init__(self, orbit_filename, attitude_filename):
"""
**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
"""
assert isinstance(orbit_filename, str)
assert isinstance(attitude_filename, str)
self.orbit_filename = orbit_filename
self.attitude_filename = attitude_filename
self._load()
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
"""
self.logger = logging.getLogger(__name__)
self.logger.info('creating an instance of EclairsOrbit')
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 J2000',
'velocity': '(VX,VY,VZ) velocity vector',
'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, quat, pos, n):
"""Set quaternion and position of n-th row in astorpy table.
For tests purposes: likely to be removed in the end.
.. warning: to be removed in the end.
:param quat: quaternion (q0, q1, q2, q3)
:type quat: array
:param pos: instrument position (X, Y, Z)
:type pos: array
:param n: row index in astropy table
:type n: int
"""
self.data['qparam'][n] = quat
self.data['position'][n] = pos
[docs] def get_earth_position(self, n):
"""Return Earth angular position (theta_E, phi_E)
in field of view.
Read positions in fits file are in the J2000 ref frame.
Need to take conjugate as matrices in algorithm are transposed.
See section 9.1 (p.50/65) in SV-SY-NT-242-JPO document (v4).
See mail/code from Sujay 17 Sept 2019
:param n: row index in astropy table
:type n: int
:return: theta and phi of the Earth center in ECLAIRs frame
:rtype: array(float) of shape (2,)
"""
assert n<=self.nrows
sat_gcrs = self.data['position'][n]
Q = self.data['qparam'][n]
u_ecllos_earth = from_j2000_to_ecllos(sat_gcrs, Q, mode='Earth')
return xyz_to_thetaphi(u_ecllos_earth)
[docs] def get_limb_angle(self, n):
"""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)
.. note: Assumes the position vector is NOT normalized
:param n: row index in astropy table
:type n: int
:return: limb angle in degrees
:rtype: float
:return: Earth center to satellite distance in m
:rtype: float
"""
assert n<=self.nrows
sat_gcrs = self.data['position'][n]
Q = self.data['qparam'][n]
u_ecllos_earth = from_j2000_to_ecllos(sat_gcrs, Q)
d_earth_sat = np.linalg.norm(u_ecllos_earth)
arg = (1e-3 * R_earth.value) / (d_earth_sat + (1e-3 * R_earth.value))
limb_angle = np.rad2deg(np.arcsin(arg))
return limb_angle, d_earth_sat
[docs] def is_earth_in_fov(self, n):
"""Check whether Earth appears in the instrument FOV.
See AlG slides for scheme and maths.
:param n: row index in astropy table
:type n: int
:return: if Earth is in FOV
:rtype: bool
"""
# compute ECLAIRs instrument half solid angle
ecl_geom = ECLAIRsGeom()
arg = (ecl_geom.massize + ecl_geom.distdetmask) / (np.sqrt(2)* ecl_geom.distdetmask)
omega_half = np.rad2deg(np.arctan(arg))
limb_angle, _ = self.get_limb_angle(n)
theta_earth, _ = self.get_earth_position(n)
return np.abs(theta_earth) <= limb_angle + omega_half
[docs] def create_open_fov_map(self, n):
"""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 and 1 otherwise.
:param n: row index in astropy table
:type n: int
:return: open fov 199x199 map
:rtype: 2D array(bool)
"""
instru_ecl = InstruECLAIRs()
size_fov_map = 199
open_fov_map = np.empty((size_fov_map, size_fov_map),
dtype=np.int8)
# compute Earth vector.
theta_earth, phi_earth = self.get_earth_position(n)
u_ecllos_earth = direlev_to_xyz(np.array([phi_earth, theta_earth]))
# compute limb angle
limb_angle, _ = self.get_limb_angle(n)
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(idx_y, idx_z)
ptg_ij = direlev_to_xyz(np.array([dir_ij, elev_ij]))
# check condition and fill open fov map
bool_cond = np.cos(np.dot(ptg_ij, u_ecllos_earth)) < np.cos(np.rad2deg(limb_angle))
open_fov_map[idx_y,idx_z] = bool_cond
return open_fov_map
[docs] def is_in_saa(self):
"""Check whether the satellite is in the South Atlantic Anomaly (SAA).
Likely to be provided by a .fits file. No need to be recomputed.
"""
raise NotImplementedError
[docs]def from_j2000_to_ecllos(v, q_j2000_sat, mode='limb'):
"""Succesive transformations of a vector v in J2000 frame
to ECLAIRs frame with appropriate quaternion.
:param v: array of unit cartesian vector in J2000 frame
:type v: numpy array with shape (3)
:param q_j2000_sat: quaternion in J2000 frame
:type q_j2000_sat: numpy array with shape (4)
:param mode: Earth or limb. Default is limb
:type mode: string
:return: array of unit cartesian vector in ECLAIRs frame
:rtype: numpy array with shape (3)
"""
# transform earth direction from J2000 to sat/VT frame.
att_vt = AttitudeVT()
att_vt.set_attitude_svom_quater(q_j2000_sat)
u_vt_earth = tq.rotate_vector(-v, tq.qconjugate(q_j2000_sat))
if mode is 'Earth':
u_vt_earth = tu.normalized_vector(u_vt_earth)
# transform from sat/VT to ECLLoS frame.
# For the moment this transformation is identity
att_ecl = AttitudeECLAIRs()
q_sat_ecllos = att_ecl._q_grf_instf
u_ecllos_earth = tq.rotate_vector(u_vt_earth, tq.qconjugate(q_sat_ecllos))
if mode is 'Earth':
u_ecllos_earth = tu.normalized_vector(u_ecllos_earth)
return u_ecllos_earth