'''
Created on 16 mai 2018
@author: Colley Jean-Marc, APC/IN2P3/CNRS
'''
import logging
import numpy as np
import transforms3d.quaternions as tq
import transforms3d.euler as te
import astropy.coordinates as coord
logger = logging.getLogger(__name__)
[docs]def xyz_to_direlev(p_xyz, deg=False):
"""Convert unit cartesian vector in spherical coordinate with CNES convention
CNES convention :
* elevation angle is the angle between vector and plane (y,z)
* direction is angle between y axis and projection vector in plane (y,z)
:param p_xyz: array of unit cartesian vector
:type p_xyz: numpy array with shape (n,3)
:return : array of spherical coordinate direction and elevation in degree
:rtype :numpy array with shape (n,2)
"""
dir_elev = np.empty((p_xyz.shape[0], 2), dtype=np.float64)
dir_elev[:, 0] = np.arctan2(p_xyz[:, 2], p_xyz[:, 1])
assert np.all(p_xyz[:, 0] <= np.float64(1))
dir_elev[:, 1] = np.arcsin(p_xyz[:, 0])
if deg:
return np.rad2deg(dir_elev)
else:
return dir_elev
[docs]def xyz_to_thetaphi(p_xyz):
"""Convert unit cartesian vector to physical spherical coordinate convention
Spherical coordinate convention:
* theta angle is the angle between the pointing vector X and the earth center
* phi angle is the angle between y axis and projection vector in plane (y,z)
:param p_xyz: array of unit cartesian vector
:type p_xyz: numpy array with shape (3)
:return : theta, phi in degree
:rtype :numpy array with shape (2)
"""
theta = np.arccos(p_xyz[0])
phi = np.arctan2(p_xyz[2], p_xyz[1])
if phi < 0:
phi += 2 * np.pi
assert 0 <= theta <= np.pi
assert 0 <= phi <= 2 * np.pi
return np.rad2deg(theta), np.rad2deg(phi)
[docs]def thetaphi_to_direlev(p_thetaphi):
"""Assume theta and phi are given in deg
"""
theta_, phi_ = np.deg2rad(p_thetaphi)
dir_ = phi_
elev_ = (np.pi / 2) - theta_
assert 0 <= elev_ <= 90
assert -180 <= dir_ <= 180
return np.rad2deg(dir_), np.rad2deg(elev_)
[docs]def direlev_to_xyz(p_direlev, deg=False):
"""Convert spherical coordinate with CNES convention to unit cartesian vector
CNES convention :
* elevation angle is the angle between vector and plane (y,z)
* direction is angle between y axis and projection vector in plane (y,z)
:param p_direlev: array of spherical coordinate direction and elevation
:type p_direlev: numpy array with shape (n,2)
:return : array of unit cartesian vector
:rtype : numpy array with shape (n,3)
"""
if deg:
direlev = np.deg2rad(p_direlev)
else:
direlev = p_direlev
if len(p_direlev.shape) == 1:
c_d = np.cos(direlev[0])
s_d = np.sin(direlev[0])
c_e = np.cos(direlev[1])
s_e = np.sin(direlev[1])
return np.array([s_e, c_e * c_d, c_e * s_d])
else:
xyz = np.empty((p_direlev.shape[0], 3), dtype=np.float64)
c_d = np.cos(direlev[:, 0])
s_d = np.sin(direlev[:, 0])
c_e = np.cos(direlev[:, 1])
s_e = np.sin(direlev[:, 1])
xyz[:, 0] = s_e
xyz[:, 1] = c_e * c_d
xyz[:, 2] = c_e * s_d
# xyz[:,0] = c_e
# xyz[:,1] = s_e*s_d
# xyz[:,2] = s_e*c_d
logger.debug(f'XYZ shape: {xyz.shape}')
return xyz
[docs]def compute_local_frame_on_sphere(p_ptg_radec, deg=False):
"""Return base of tangent local frame of unit sphere for
a given direction with astronomer orientation convention
:param p_ptg_radec: array ra, dec
:type p_ptg_radec: numpy array float
:return : 3 3d array : e1, e2, e3
:rtype : 3 3d numpy array
"""
if deg:
radec = np.deg2rad(p_ptg_radec)
else:
radec = p_ptg_radec
c_ra = np.cos(radec[0])
s_ra = np.sin(radec[0])
s_dec = np.sin(radec[1])
c_dec = np.cos(radec[1])
# # local frame on celestial sphere to compute orientation
# e1 to North
e1 = np.array([-s_dec * c_ra, -s_dec * s_ra, c_dec], dtype=np.float64)
logger.debug(f'e1: {e1}')
# e3 inverse pointing
# astronomer orientation convention
e3 = np.array([-c_dec * c_ra, -c_dec * s_ra, -s_dec], dtype=np.float64)
logger.debug(f'e3: {e3}')
# e2 to East
e2 = np.cross(e3, e1)
return e1, e2, e3
[docs]def ptgori_to_quaternion(p_ptg_radec, p_ori, deg=False):
"""Compute quaternion associated to rotation defined by a pointing direction and orientation
X instru axis is the pointing axis , ** CNES choice **
Z instru axis determine orientation in sky, not clear in CNES document, so it's my choice
Orientation of Z instru axis follow IAU recommandation to measure polarisation angle
example:
1) orientation is 0 so Z instru axis pointed to celestial North pole
2) orientation is 90 deg so Z instru axis pointed
following a parallel of celestial sphere with East direction
return : quaternion attitude
:param p_ptg_radec: ra, dec
:type p_ptg_radec: list
:param p_ori: orientation angle with convention defined above
:type p_ori: number
:return : unit quaternion
:rtype : list/array of float
"""
if deg:
radec = np.deg2rad(p_ptg_radec)
else:
radec = p_ptg_radec
e1, e2, e3 = compute_local_frame_on_sphere(radec)
r1 = -e3
r3 = np.cos(p_ori) * e1 + np.sin(p_ori) * e2
r2 = np.cross(r3, r1)
rot_attitude = np.empty((3, 3), dtype=np.float64)
rot_attitude[:, 0] = r1
rot_attitude[:, 1] = r2
rot_attitude[:, 2] = r3
logger.debug(f'rot_attitude: {rot_attitude}')
q_attitude = tq.mat2quat(rot_attitude)
return q_attitude
[docs]class SVOMPointingError(Exception):
""" Unused class
"""
pass
[docs]class AttitudeVT(object):
"""Attitude for VT instrument
"""
def __init__(self, p_frame_cat='gcrs', no_biais=False):
"""
doc CNES 6.1.1 inertial reference frame
For SVOM applications, J2000 and GCRF can be considered as identical,
J2000 is considered the baseline coordination frame.
:param p_frame_cat: astropy keyword to defined a frame
:type p_frame_cat: string
"""
self.catsmf = 'gcrs' # astronomical star mapper frame
self.catuf = p_frame_cat # astronomical user frame
# todo : check p_frame_cat astropy key word and raise SVOMPointingError is problem
self.instru_name = "VT"
if no_biais:
self._q_grf_instf = np.array([1.0, 0, 0, 0], dtype=np.float64)
else:
self._set_rot_instf_to_grf()
def _set_rot_instf_to_grf(self):
"""
set quaternion instrument frame to svom frame
VTlos is svom frame, see CNES documentation "8.2 AAV bulletin"
GRF = Guidance Reference Frame (which is VT sighting reference frame R VTlos as a
baseline, see CSC-R-4.2-0150)
"""
# so q_grf_instf is identity quaternion
self._q_grf_instf = np.array([1.0, 0, 0, 0], dtype=np.float64)
[docs] def set_attitude_svom_quater(self, p_quater):
"""
Set attitude svom satellite with quatenion in the same convention as the message AAV
Quaternion given by AAV bulletin: J2000->GRF
1) with symbol '->' convention CNES documentation "3.3 quaternion"
* x [B] = M B->A x [A]
where M is rotation matrix
* Q R1->R3 = Q R1->R2 . Q R2->R3
2) GRF = Guidance Reference Frame (which is VT sighting reference frame R VTlos as a
baseline, see CSC-R-4.2-0150)
3) For SVOM applications, J2000 and GCRF can be considered as identical, J2000 is considered
the baseline coordination frame.
M out->in
:param p_quater: unit quaternion
:type p_quater: list/array
"""
assert len(p_quater) == 4
self._q_j2000_grf = p_quater.copy()
self._q_j2000_instf = tq.qmult(self._q_j2000_grf, self._q_grf_instf)
self._r_j2000_instf = tq.quat2mat(self._q_j2000_instf)
self._compute_ptg_instru()
self._compute_orientation_instru()
[docs] def set_attitude_instru_quater(self, p_quater):
"""Set attitude instru (no biais) satellite with quaternion in the same convention
as the message AAV
Quaternion given by AAV bulletin: J2000->GRF
1) with symbol '->' convention CNES documentation "3.3 quaternion"
* x [B] = M B->A x [A]
where M is rotation matrix
* Q R1->R3 = Q R1->R2 . Q R2->R3
2) GRF = Guidance Reference Frame (which is VT sighting reference frame R VTlos as a
baseline, see CSC-R-4.2-0150)
3) For SVOM applications, J2000 and GCRF can be considered as identical, J2000 is considered
the baseline coordination frame.
M out->in
:param p_quater: unit quaternion
:type p_quater: list/array
"""
q_j2000_grf = tq.qmult(p_quater, tq.qinverse(self._q_grf_instf))
self.set_attitude_svom_quater(q_j2000_grf)
[docs] def set_attitude_svom_ptgori(self, p_ptg_ori):
"""Set attitude SVOM with pointing vector and orientation
:param p_ptg_ori: ra, dec, ori in degree
:type p_ptg_ori: list
"""
ptg_ori_rad = np.deg2rad(p_ptg_ori)
q_j2000_grf = ptgori_to_quaternion(ptg_ori_rad[:2], ptg_ori_rad[2])
self.set_attitude_svom_quater(q_j2000_grf)
[docs] def set_attitude_instru_ptgori(self, p_ptg_ori):
"""Set attitude instrument (no biais) with pointing vector and orientation
:param p_ptg_ori: ra, dec, ori in degree
:type p_ptg_ori: list
"""
ptg_ori_rad = np.deg2rad(p_ptg_ori)
q_j2000_instf = ptgori_to_quaternion(ptg_ori_rad[:2], ptg_ori_rad[2])
q_j2000_grf = tq.qmult(q_j2000_instf, tq.qinverse(self._q_grf_instf))
self.set_attitude_svom_quater(q_j2000_grf)
[docs] def instru_to_cat_user(self, p_ptg_instruf, deg=False):
"""
in : xyz (n,3)
out : ra dec (n,2)
:param p_ptg_instruf: pointing in unit cartesian
:type p_ptg_instruf: numpy array with shape (n,3)
:return : ra ,dec
:rtype : numpy array with shape (n,2)
"""
assert p_ptg_instruf.shape[1] == 3
nb_elt = p_ptg_instruf.shape[0]
logger.debug(f'nb_elt: {nb_elt}')
ptg_j2000_c = np.dot(self._r_j2000_instf, np.transpose(p_ptg_instruf))
ptg = coord.SkyCoord(x=ptg_j2000_c[0, :],
y=ptg_j2000_c[1, :],
z=ptg_j2000_c[2, :],
frame=self.catsmf,
representation_type='cartesian')
if self.catuf == self.catsmf:
ptg_out = ptg
else:
ptg_out = ptg.transform_to(self.catuf)
ptg_out.representation_type = 'unitspherical'
ra_dec_out = np.empty((nb_elt, 2), dtype=np.float64)
if deg:
ra_dec_out[:, 0] = ptg_out.ra.to('deg')
ra_dec_out[:, 1] = ptg_out.dec.to('deg')
else:
ra_dec_out[:, 0] = ptg_out.ra.to('rad')
ra_dec_out[:, 1] = ptg_out.dec.to('rad')
return ra_dec_out
[docs] def cat_user_to_instru(self, p_ptg_catuf, deg=False):
"""
in : ra dec (n,2) [deg]
out : xyz (n,3)
:param p_ptg_catuf: ra ,dec
:type p_ptg_catuf: numpy array with shape (n,2)
:return : pointing in unit cartesian
:rtype : numpy array with shape (n,3)
"""
assert p_ptg_catuf.shape[1] == 2
if deg:
a_unit = 'deg'
else:
a_unit = 'rad'
ptg_u = coord.SkyCoord(ra=p_ptg_catuf[:, 0], dec=p_ptg_catuf[:, 1], \
unit=a_unit, frame=self.catuf)
# foo.cartesian to convert in cartesian
# foo.xyz to return numpy array (3,)
if self.catuf == self.catsmf:
ptg_u.representation_type = 'cartesian'
ptg_j2000_c = ptg_u.cartesian.xyz
else:
ptg_j2000_c = ptg_u.transform_to(self.catsmf).cartesian.xyz
ptg_inst_c = np.transpose(np.dot(np.transpose(self._r_j2000_instf), ptg_j2000_c))
assert ptg_inst_c.shape[1] == 3
sh_0 = ptg_inst_c.shape
# print(sh_0)
# return with unit normalization
ptg_inst_c /= np.linalg.norm(ptg_inst_c, axis=1)[:, np.newaxis]
assert ptg_inst_c.shape == sh_0
return ptg_inst_c
def _compute_ptg_instru(self):
"""
X instru axis is the pointing axis , ** CNES choice **
"""
self.ptg_radec = self.instru_to_cat_user(np.array([[1, 0, 0]], dtype=np.float64))
[docs] def ptg_instru(self, deg=False):
"""Return pointing of instrument
:return : pointing of instrument ra, dec
:rtype : numpy array with shape (1,2)
"""
if deg:
return np.rad2deg(self.ptg_radec)
else:
return self.ptg_radec
def _compute_orientation_instru(self):
"""Compute orientation of instrument
Orientation convention is IAU/polarisation in rad.
fill attribut : ori_rad
"""
self.e1, self.e2, self.e3 = compute_local_frame_on_sphere(self.ptg_radec[0])
# oriention is angle of z in e1 e2
z_radec = self.instru_to_cat_user(np.array([[0, 0, 1]], dtype=np.float64))
c_ra = np.cos(z_radec[0, 0])
s_ra = np.sin(z_radec[0, 0])
s_dec = np.sin(z_radec[0, 1])
c_dec = np.cos(z_radec[0, 1])
z_c = np.array([c_dec * c_ra, c_dec * s_ra, s_dec], dtype=np.float64)
c_ori = np.dot(z_c, self.e1)
s_ori = np.dot(z_c, self.e2)
self.ori_rad = np.arctan2(s_ori, c_ori)
[docs] def ori_instru(self, deg=False):
""" return orientation of instrument in IAU/pol convention
:return : angle orientation
:rtype : float
"""
if deg:
return np.rad2deg(self.ori_rad)
else:
return self.ori_rad
[docs]class AttitudeECLAIRs(AttitudeVT):
"""Attitude for ECLAIRs instrument
"""
def __init__(self, p_frame_cat='gcrs', no_biais=False):
'''**Constructor**
'''
super().__init__(p_frame_cat, no_biais)
self.instru_name = "ECLAIRs"
def _set_rot_instf_to_grf(self):
"""set quaternion instrument frame to svom frame
VTlos is svom frame, see CNES documentation
"8.2 AAV bulletin"
GRF = Guidance Reference Frame (which is VT sighting reference frame R VTlos as a
baseline, see CSC-R-4.2-0150)
SV-FD-ECL-002 page 51
"""
# todo biais matrix between frame VTLos and ECLAIRsLos
# q_grf_sl = te.euler2quat(np.deg2rad(0.001), np.deg2rad(-0.002), np.deg2rad(0.0005), 'szyx')
# q_sl_ecl_los = te.euler2quat(np.deg2rad(-0.001), np.deg2rad(0.0005), np.deg2rad(0.001), 'szyx')
# self._q_grf_instf = tq.qmult(q_grf_sl, q_sl_ecl_los)
self._q_grf_instf = np.array([1.0, 0, 0, 0], dtype=np.float64)
[docs]class AttitudeMXT(AttitudeVT):
"""Attitude for MXT instrument
"""
def __init__(self, p_frame_cat='gcrs', no_biais=False):
'''**Constructor**
'''
super().__init__(p_frame_cat, no_biais)
self.instru_name = "MXT"
def _set_rot_instf_to_grf(self):
"""set quaternion instrument frame to svom frame
VTlos is svom frame, see CNES documentation :
"8.2 AAV bulletin"
GRF = Guidance Reference Frame (which is VT sighting reference frame R VTlos as a
baseline, see CSC-R-4.2-0150)
SV-FD-MXT-003 page 54
"""
# todo biais matrix between frame VTLos and MXTlos
# for example q_grf_mxt_los = te.euler2quat(np.deg2rad(0.0005), np.deg2rad(-0.002), np.deg2rad(0.001), 'szyx')
q_grf_mxt_los = np.array([1.0, 0, 0, 0], dtype=np.float64)
self._q_grf_instf = q_grf_mxt_los