'''
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):
"""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])
dir_elev[:,1] = np.arcsin(p_xyz[:,0])
return np.rad2deg(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
return np.rad2deg(theta), np.rad2deg(phi)
[docs]def direlev_to_xyz(p_direlev):
"""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 in degree
:type p_direlev: numpy array with shape (n,2)
:return : array of unit cartesian vector
:rtype : numpy array with shape (n,3)
"""
if len(p_direlev.shape)==1:
direlev = np.deg2rad(p_direlev)
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)
direlev = np.deg2rad(p_direlev)
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
logger.debug(f'XYZ shape: {xyz.shape}')
return xyz
[docs]def compute_local_frame_on_sphere(p_ptg_radec):
"""Return base of tangent local frame of unit sphere for
a given direction with astronomer orientation convention
:param p_ptg_radec: array ra, dec in radian
:type p_ptg_radec: numpy array float
:return : 3 3d array : e1, e2, e3
:rtype : 3 3d numpy array
"""
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):
"""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 in rad
: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
"""
e1, e2, e3 = compute_local_frame_on_sphere(p_ptg_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):
pass
[docs]class AttitudeVT(object):
"""
"""
def __init__(self, p_frame_cat='gcrs'):
"""
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"
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
"""
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_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 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):
"""
in : xyz (n,3)
out : ra dec (n,2) [deg]
:param p_ptg_instruf: pointing in unit cartesian
:type p_ptg_instruf: numpy array with shape (n,3)
:return : ra ,dec in degree
: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)
ra_dec_out[:,0] = ptg_out.ra.to('deg')
ra_dec_out[:,1] = ptg_out.dec.to('deg')
return ra_dec_out
[docs] def cat_user_to_instru(self, p_ptg_catuf):
"""
in : ra dec (n,2) [deg]
out : xyz (n,3)
:param p_ptg_catuf: ra ,dec in degree
: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
ptg_u = coord.SkyCoord(ra=p_ptg_catuf[:,0], dec=p_ptg_catuf[:,1], unit='deg', 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
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):
"""Return pointing of instrument
:return : pointing of instrument ra, dec in degree
:rtype : numpy array with shape (1,2)
"""
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(np.deg2rad(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))
z_radec = np.deg2rad(z_radec)
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):
""" return orientation of instrument in IAU/pol convention
:return : angle orientation in rad
:rtype : float
"""
return self.ori_rad
[docs] def ori_instru_deg(self):
""" return orientation of instrument in IAU/pol convention
:return : angle orientation in degree
:rtype : float
"""
return np.rad2deg(self.ori_rad)
[docs]class AttitudeECLAIRs(AttitudeVT):
def __init__(self, p_frame_cat='gcrs'):
super().__init__(p_frame_cat)
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)
[docs]class AttitudeMXT(AttitudeVT):
def __init__(self, p_frame_cat='gcrs'):
super().__init__(p_frame_cat)
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
q_grf_mxt_los = te.euler2quat(np.deg2rad(0.0005), np.deg2rad(-0.002), np.deg2rad(0.001), 'szyx')
self._q_grf_instf = q_grf_mxt_los