# -*-coding:utf-8 -*
'''
Description:
A wrapper class handling Rosetta SPICE ancillary data. RPC-ICA specific.
To create one object::
SPICESet = ICASPICEData(dates, additionnalPath='/<bli>/<bla>/ICAPy/') ## (Absolute) path to the ICAPy directory.
SPICESet.initTimeDt('20141128' ,'20141204' , N=1000) ## Initialize time from date 1 to date 2 with 1000 points.
OR::
SPICESet.timeDt = ICADataSet.timeDtScan ## Synchronizes the SPICE data with ICA data (here FOVs).
Arguments:
additionnalPath (default ’’) Where SPICE related files can be found (metakernel, kernels themselves, 3D nucleus shape optionally, etc)
metakernelFile (default ‘ros_kernels_3dtool.txt’) The name of the metakernel file (not the path).
Dependencies:
irfpy.ica.modules
irfpy.ica.tools
SpiceyPy
'''
from irfpy.ica.modules import *
import irfpy.ica.tools as icaT
#import spiceypy as spice
try:
import SpiceyPy as spice
except ImportError:
import spiceypy as spice
[docs]class SPICEClass:
def __init__(self, additionnalPath=rcSPICEPath, metakernelFile='ros_kernels_3dtool.txt'):
print('\nSPICEClass________________________________________________________________\n')
if additionnalPath[-1] == '/':
self.additionnalPath = additionnalPath
else:
self.additionnalPath = additionnalPath + '/'
self.metakernelFile = metakernelFile
spice.furnsh(self.additionnalPath + self.metakernelFile)
if spice.failed():
sys.exit('\nUnable to load metakernel. Game over.\n')
# Date format. This is the master vector: everything will be initialized
# according to this one.
self.timeDt = None
# Can be either directly specified by the user (to sync SPICE informations with actual data)
# or using the self.initTimeDt(t1,t2,N) (see bellow).
self.timeET = None
self.sunElev = None # RADIAN.
self.sunAzim = None # RADIAN.
self.cometElev = None # RADIAN.
self.cometAzim = None # RADIAN.
# Expressed in the ICA reference frame (Xsc = Xica, Ysc = Zica, Zsc = -Yica).
self.sunDirection_ICA = None
self.cometDirection_ICA = None
# Expressed in CSEQ reference frame. For the sun, only for checking purposes..
self.sunDirection_CSEQ = None
self.cometDirection_CSEQ = None
self.phySA = None # Solar array angle. RADIAN.
self.distanceSun = None # Helio centrique distance. Initialised in ICASPICEData.initEphemSunJ2000()
self.AU = 149597871. # Astronomical unit in km.
self.ephem_Ros_Chu_Chu = None # Ros seen by Chu in Chu-FIXED.
self.ephem_Ros_Chu_J2000 = None # Ros seen by Chu in J2000.
self.ephem_Sun_Chu_J2000 = None # Sun seen by Chu in J2000.
self.ephem_Ros_Sun_J2000 = None # Ros seen by Sun in J2000.
self.ephem_Chu_Sun_J2000 = None # Chu seen by Sun in J2000.
self.ephem_Ear_Sun_J2000 = None # Earth seen by Sun in J2000.
self.ephem_Mar_Sun_J2000 = None # Mars seen by Sun in J2000.
self.ephem_Ros_Chu_CSO = None # Ros seen by Chu in CSO, initialized in initEphemCSO.
self.ephem_Sun_Chu_CSO = None # Check purpose.
# Radius, longitude and lattitude in CSO, initialized in initEphemCSO.
self.radLonLat_CSO = None
self.ephem_Ros_Chu_CSEQ = None # Ros seen by Chu in CSEq, initialized in initEphemCSEQ.
self.ephem_Ros_Sun_CSEQ = None # Ros seen by Sun in CSEq, initialized in initEphemCSEQ.
self.ephem_Ros_Chu_CK = None # Ros seen by Sun in C-G_FIXED, initialized in initEphemFIXED.
# Radius, longitude and lattitude in CSEQ, initialized in initEphemCSEQ.
self.radLonLat_CSEQ = None
self.radLonLat_CK = None
self.rotate_CSEQ_ICA = None # Rotation matrix --FROM-- CSEQ --TO-- ICA .
self.rotate_ICA_CSEQ = None # Rotation matrix --FROM-- ICA --TO-- CSEQ .
self.rotate_SC_CSEQ = None # Rotation matrix --FROM-- SC --TO-- CSEQ .
self.rotate_CSEQ_SC = None # Rotation matrix --FROM-- CSEQ --TO-- SC .
self.rotate_ICA_CK = None # Rotation matrix --FROM-- ICA --TO-- FIXED .
self.rotate_CK_ICA = None # Rotation matrix --FROM-- FIXED --TO-- ICA .
##########################################################################
##########################################################################
##
# From 'http://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/pxform_c.html'
##
# rotate is the matrix that transforms position vectors from the
# reference frame `from' to the frame `to' at epoch `et'.
# If (x, y, z) is a position relative to the frame `from'
# then the vector ( x', y', z') is the same position
# relative to the frame `to' at epoch `et'. Here the
# vector ( x', y', z' ) is defined by the equation:
##
# - - - - - -
# | x' | | | | x |
# | y' | = | rotate | | y |
# | z' | | | | z |
# - - - - - -
##
##########################################################################
##########################################################################
[docs] def unloadAll(self):
r''' Unload all kernels.'''
spice.unload(self.additionnalPath + self.metakernelFile)
spice.unload('/Users/etiennebehar/Documents/Rosetta/Data/SPICE/kernels/fk/ROS_V26.TF')
[docs] def initTimeDt(self, t1, t2, N=100):
r''' Initialize ICASPICEData.timeDt with
t1: 'yyyymmdd' ,
t2: 'yyyymmdd' ,
N: number of points. '''
t1 = datetime.date(int(t1[:4]), int(t1[4:6]), int(t1[6:8])).toordinal()
t2 = datetime.date(int(t2[:4]), int(t2[4:6]), int(t2[6:8])).toordinal()
tt = np.linspace(t1, t2, N)
self.timeDt = np.array(num2date(tt))
# def initEphemChuCK(self):
# r''' Initialize ephemeris in the comet-FIXED reference frame, according to ICASPICEData.timeDt . '''
# self.timeET = self.dtToET(self.timeDt)
# self.ephem_Ros_Chu_Chu = np.array([spice.spkezr('-226', t, '67P/C-G_CK', 'None', '1000012')[0][:3] for t in self.timeET]).T # SC seen by Chury in Chury FIXED frame
[docs] def initEphemChuJ2000(self):
r''' Initialize ephemeris in J2000 centered on 67P, according to ICASPICEData.timeDt . '''
self.timeET = self.dtToET(self.timeDt)
self.ephem_Ros_Chu_J2000 = np.array([spice.spkezr(
'-226', t, 'J2000', 'None', '1000012')[0][:3] for t in self.timeET]).T # SC seen by Chury in J2000
self.ephem_Sun_Chu_J2000 = np.array([spice.spkezr('SUN', t, 'J2000', 'None', '1000012')[0][:3] for t in self.timeET]).T # Sun seen by Chury in J2000
[docs] def initEphemSunJ2000(self):
r''' Initialize ephemeris in J2000 centered on the Sun, according to ICASPICEData.timeDt . '''
self.timeET = self.dtToET(self.timeDt)
# self.ephem_Ros_Sun_J2000 = np.array([spice.spkezr(
# '-226', t, 'J2000', 'None', 'SUN')[0][:3] for t in self.timeET]).T # SC seen by the Sun in J2000
# self.distanceSun = np.array([spice.reclat(v)[0]
# for v in self.ephem_Ros_Sun_J2000.T]).T / self.AU
self.ephem_Chu_Sun_J2000 = np.array([spice.spkezr('1000012', t, 'J2000', 'None', 'SUN')[0][:3] for t in self.timeET]).T # 67P seen by the Sun in J2000
self.ephem_Ear_Sun_J2000 = np.array([spice.spkezr('EARTH', t, 'J2000', 'None', 'SUN')[0][:3] for t in self.timeET]).T # 67P seen by the Sun in J2000
self.ephem_Mar_Sun_J2000 = np.array([spice.spkezr('MARS', t, 'J2000', 'None', 'SUN')[0][:3] for t in self.timeET]).T # 67P seen by the Sun in J2000
[docs] def initEphemCSO(self):
r''' Initialize ephemeris in CSO centered on 67P, according to ICASPICEData.timeDt . '''
self.timeET = self.dtToET(self.timeDt)
self.ephem_Ros_Chu_CSO = np.array([spice.spkezr('-226', t, '67P/C-G_CSO', 'None', '1000012')[0][:3] for t in self.timeET]).T # SC seen by the nucleus in CSO
self.ephem_Sun_Chu_CSO = np.array([spice.spkezr('SUN', t, '67P/C-G_CSO', 'None', '1000012')[0][:3] for t in self.timeET]).T # Sun seen by the nucleus in CSO (check)
self.radLonLat_CSO = np.array([spice.reclat(v) for v in self.ephem_Ros_Chu_CSO.T]).T
[docs] def initEphemCSEQ(self):
r''' Initialize ephemeris in CSEQ centered on 67P, according to ICASPICEData.timeDt. '''
self.timeET = self.dtToET(self.timeDt)
self.ephem_Ros_Chu_CSEQ = np.array([spice.spkezr('-226', t, '67P/C-G_CSEQ', 'None', '1000012')[0][:3] for t in self.timeET]).T # SC seen by the nucleus in CSEQ
self.radLonLat_CSEQ = np.array([spice.reclat(v) for v in self.ephem_Ros_Chu_CSEQ.T]).T
[docs] def initEphemCK(self):
r''' Initialize ephemeris in 67P/C-G_CK centered on 67P, according to ICASPICEData.timeDt. '''
self.timeET = self.dtToET(self.timeDt)
self.ephem_Ros_Chu_CK = np.array([spice.spkezr('-226', t, '67P/C-G_CK', 'None', '1000012')[0][:3] for t in self.timeET]).T # SC seen by the nucleus in CK
self.radLonLat_CK = np.array([spice.reclat(v) for v in self.ephem_Ros_Chu_CK.T]).T
[docs] def initRotateCSEQ(self):
r''' Initialize rotation matrices, according to ICASPICEData.timeDt:
'ROS_RPC_ICA' <-> '67P/C-G_CSEQ'
'ROS_SPACECRAFT' <-> '67P/C-G_CSEQ'
'ROS_SPACECRAFT' -> 'ROS_SA-Y'
'''
self.timeET = self.dtToET(self.timeDt)# From CSEQ to ICA, WARNING, see Bernhart document.
self.rotate_CSEQ_ICA = np.array([spice.pxform('67P/C-G_CSEQ', 'ROS_RPC_ICA', t) for t in self.timeET]) ## From ICA to CSEQ, WARNING, see Bernhart document.
self.rotate_ICA_CSEQ = np.array([spice.pxform('ROS_RPC_ICA', '67P/C-G_CSEQ', t) for t in self.timeET]) ## From CSEq to ICA.
# self.rotate_CSEQ_SC = np.array( [spice.pxform('67P/C-G_CSEQ', 'ROS_SPACECRAFT', self.dtToET(t)) for t in self.timeDt] ) ## From CSEQ to SC.
# self.rotate_SC_CSEQ = np.array( [spice.pxform('ROS_SPACECRAFT', '67P/C-G_CSEQ', self.dtToET(t)) for t in self.timeDt] ) ## From SC to CSEQ.
self.phySA = np.array([np.arccos(np.array(spice.pxform('ROS_SPACECRAFT', 'ROS_SA-Y', self.dtToET(t)))[0, 0] * -1) for t in self.timeDt])
# self.phySA = np.zeros(self.timeDt.size)
[docs] def initRotateCK(self):
r''' Initialize rotation matrices, according to ICASPICEData.timeDt:
'ROS_RPC_ICA' <-> '67P/C-G_CK'
'ROS_SPACECRAFT' -> 'ROS_SA-Y'
'''
self.timeET = self.dtToET(self.timeDt)# From CSEQ to ICA, WARNING, see Bernhart document.
self.rotate_CK_ICA = np.array([spice.pxform('67P/C-G_CK', 'ROS_RPC_ICA', t) for t in self.timeET]) ## From ICA to CK, WARNING, see Bernhart document.
self.rotate_ICA_FIXED = np.array([spice.pxform('ROS_RPC_ICA', '67P/C-G_CK', t) for t in self.timeET]) ## From CK to ICA.
self.phySA = np.array([np.arccos(np.array(spice.pxform('ROS_SPACECRAFT', 'ROS_SA-Y', self.dtToET(t)))[0, 0] * -1) for t in self.timeDt])
# self.phySA = np.zeros(self.timeDt.size)
[docs] def initSunCometDirection(self):
r''' Initialize direction of the Sun and the nucleus in ICA reference frame, according to ICASPICEData.timeDt. '''
self.timeET = self.dtToET(self.timeDt)
# Sun in ICA reference frame seen by SC.
self.sunDirection_ICA = np.array([spice.spkezr('SUN', t, 'ROS_RPC_ICA', 'None', '-226')[0][0:3] for t in self.timeET]).T
# self.sunDirection_CSEQ = np.array([spice.spkezr( 'SUN', t,
# '67P/C-G_CSEQ', 'None', '-226')[0][0:3] for t in self.timeET]).T ##
# Check..
self.sunAzim, self.sunElev = self.vec2ang(self.sunDirection_ICA.copy())
# Comet in ICA reference frame seen by SC.
self.cometDirection_ICA = np.array([spice.spkezr('1000012', t, 'ROS_RPC_ICA', 'None', '-226')[0][0:3] for t in self.timeET]).T
# Comet in CSEq seen by spacecraft.
self.cometDirection_CSEQ = np.array([spice.spkezr('1000012', t, '67P/C-G_CSEQ', 'None', '-226')[0][0:3] for t in self.timeET]).T
self.cometAzim, self.cometElev = self.vec2ang(self.cometDirection_ICA.copy())
[docs] def changeFrame_ICA_to_CSEQ(self, v):
r''' Rotate vector v from ICA frame to CSEQ frame. v.shape = 3 x n '''
return np.array([np.dot(M, v[:, i]) for i, M in enumerate(self.rotate_ICA_CSEQ)]).T
[docs] def vec2ang(self, v):
''' From a vector or an array of vector to (azim,elev) angles.
v.shape = 3xN'''
v /= icaT.norm(v)
elev = np.arcsin(v[2])
azim = np.arccos(v[0] / np.cos(elev))
azim[v[1] < 0] = 2 * np.pi - azim[v[1] < 0]
return azim, elev
[docs] def dtToET(self, t):
if type(t) is datetime.datetime:
return spice.str2et(t.strftime('%FT%T'))
elif type(t) is np.ndarray:
return np.fromiter(map(lambda x: spice.str2et(x.strftime('%FT%T')), t), dtype=float)
[docs] def ordinalToET(self, t):
return spice.str2et(num2date(t).strftime('%FT%T'))