''' MEX spice related module
This module provide the MEX-related spice access simpler.
``.irfpyrc`` file should have following entry.
..
[mexpvat]
spiceroot = /Volumes/scidata/data/mars/spice/kernels/
The above path is looked through as root path of SPICE kernel.
There are three levels of the kernel set.
- 'full' kernel set is to load all the kernels.
- 'noatt' kernel set does not load attitude file. One may use this
option if one only wants the orbit data.
- 'min' kernel set provide minimum set, namely no attitude/no orbit
but for example, Mars position or spacecraft frames are defined.
Obviously if you choose 'min' set, it is much faster than 'full' set
for calculation.
>>> init(level='min')
You may overwrite the level. In this case, only the kernel not to
be furnshed is loaded.
>>> init(level='noatt')
Each file can also be furnshed.
.. code-block:: python
mexspice.furnsh('sp/EARTH.BSP')
A sample code is seen in :mod:`snippet.mexspicesample`.
.. todo::
SpiceFactory should not check the kernels if they have been loaded.
In SPICE, the last kernel loaded is in priority,
user have to add specific kernel in the end, if necessarily.
.. note::
Since March 2017, ESA looks to provide MK file, which contains the
meta data of SPICE. The current version just ignores this file.
For the future development, the MK file will be used.
'''
import os
import logging
import datetime
import dateutil.parser
import re
import operator
import numpy as np
import spiceypy as spice
import irfpy.util.irfpyrc
import irfpy.util.exception as ex
logger = logging.getLogger(__name__)
[docs]def isdb():
rc = irfpy.util.irfpyrc.Rc()
spiceroot = rc.get('mexpvat', 'spiceroot')
if spiceroot is None:
return False
if not os.path.isdir(spiceroot):
return False
return True
[docs]class Attfiles:
r''' Attitude file handling class.
Intention is to implement the automatic furnsh routine for
attitude file from time, but it may need a moderate job.
Now, the implementaion of that functionality is on hold.
.. todo::
Implement the auto-load function of attitude file
(and orbit file) from the given time.
'''
[docs] class Attfile:
def __init__(self, fn):
sea = re.search(r'ATNM_PTR(\d\d\d\d\d)_(\d\d\d\d\d\d)_(\d\d\d)\.BC', fn)
if sea is None:
raise RuntimeError('Not a good file name convention.')
ptr, t0, ver = sea.groups()
self.ptr = int(ptr, 10)
self.t0 = dateutil.parser.parse(t0)
self.ver = int(ver, 10)
self.filename = fn
def __str__(self):
return '{} {} {} {}'.format(self.filename, self.ptr, self.t0, self.ver)
def __repr__(self):
return self.__str__()
def __init__(self):
self.rc = irfpy.util.irfpyrc.Rc()
self.list = self.makelist()
[docs] def makelist(self):
spiceroot = self.rc.get('mexpvat', 'spiceroot')
if spiceroot is None:
logger.warn('Specify the root of spice in .irfpyrc'
'by [mexpvat] spiceroot entry')
raise IOError('No entry of spiceroot.')
fs = []
for p, d, f in os.walk(spiceroot):
if p.find('former_versions') >= 0:
continue
for ff in f:
try:
attf = self.Attfile(os.path.join(p, ff))
fs.append(attf)
except RuntimeError:
continue
fs = sorted(fs, key=operator.attrgetter('t0'))
return fs
[docs]class MexSpiceFactory:
''' Use :func:`create` function for creating the object.
'''
def __init__(self):
'''
'''
self.rc = irfpy.util.irfpyrc.Rc()
self.loaded = [] # A fullpath list of the loaded kernels
self.exclude_kernel = ['MAR063.BSP', 'MAR080S.BSP', 'MAR085.BSP', ]
''' List of filename to remove from the get_xxxx method
'''
self.spice = spice
[docs] def showstatus(self):
for k in self.loaded:
print(' -- %s' % k)
print('Loaded kernels = %d' % len(self.loaded))
[docs] def get(self, level='full', dataonly=True):
''' Return the filename list of kernel.
:param level: "full", "noatt", or "min".
:meth:`get_full`, :meth:`get_noatt`, or :meth:`get_min`.
'''
corres = {'full': self.get_full, 'noatt': self.get_noatt, 'min': self.get_min}
if not level in corres:
raise KeyError('No %s in the kernel set.' % level)
return corres[level](dataonly=dataonly)
[docs] def get_full(self, dataonly=True):
''' Get the kernel file name list.
'''
spiceroot = self.rc.get('mexpvat', 'spiceroot')
if spiceroot is None:
logger.warning('Specify the root of spice in .irfpyrc'
'by [mexpvat] spiceroot entry')
raise IOError('No entry of spiceroot.')
logger.debug('spiceroot= %s' % spiceroot)
kfiles = []
for p, d, f in os.walk(spiceroot):
if p.find('former_versions') >= 0:
continue
for ff in f:
if ff.lower().endswith('.txt'):
continue
if ff.lower().endswith('.tm'):
continue
if ff.endswith('~'):
continue
kfiles.append(os.path.join(p, ff))
if dataonly:
logger.debug('Dataonly ON')
removes = []
for fabsname in kfiles:
fname = os.path.basename(fabsname)
if re.search(r'ORMC_FD_.*\.ORB', fname) is not None:
removes.append(fabsname)
elif re.search(r'ORMF_.*\.ORB', fname) is not None:
removes.append(fabsname)
elif re.search(r'MAR033_HRSC_.*\.BSP', fname) is not None:
removes.append(fabsname)
elif re.search(r'MARSAT_ESA_.*\.BSP', fname) is not None:
removes.append(fabsname)
elif re.search(r'MEX_ROB_.*\.BSP', fname) is not None:
removes.append(fabsname)
elif re.search(r'ORMF_.*\.BSP', fname) is not None:
removes.append(fabsname)
elif re.search(r'ORMC_FD_.*\.BSP', fname) is not None:
removes.append(fabsname)
elif fname in self.exclude_kernel:
removes.append(fabsname)
for fname in removes:
logger.debug('Remove {}'.format(fname))
kfiles.remove(fname)
# Addition of the right kernel for the comet SS encounter to the end
logger.debug('SS enounter kernel in the end')
for fabsname in kfiles:
fname = os.path.basename(fabsname)
if re.search('ATNM_PTR00667_141013_001.BC', fname) is not None:
kfiles.remove(fabsname)
kfiles.append(fabsname)
return kfiles
[docs] def get_noatt(self, dataonly=True):
''' Update the nonatt MK file.
Update MK file that does not include attitude file.
'''
kfiles = self.get_full(dataonly=dataonly)
removes = []
for fabsname in kfiles:
fname = os.path.basename(fabsname)
if re.search(r'ATNM_.*\.BC', fname) is not None:
removes.append(fabsname)
for fname in removes:
logger.debug('Remove {}'.format(fname))
kfiles.remove(fname)
return kfiles
[docs] def get_min(self, dataonly=True):
''' Update the minimum MK file.
Minimum MK file contains all the files except for orbit and attitude.
:keyword dataonly: If set to True (default), only the kernels
related to data analysis is returned. Otherwise, all the
data below the specific directory is returned.
'''
kfiles = self.get_noatt(dataonly=dataonly)
removes = []
for fabsname in kfiles:
fname = os.path.basename(fabsname)
if re.search(r'ORMM_.*\.BSP', fname):
removes.append(fabsname)
for fname in removes:
logger.debug('Remove {}'.format(fname))
kfiles.remove(fname)
return kfiles
[docs] def furnsheach(self, fname):
fabsname = os.path.abspath(fname)
if not fabsname in self.loaded:
if not os.path.exists(fabsname):
raise IOError("Kernel file {} not found".format(fabsname))
logger.info('Furnsh: {}'.format(fabsname))
spice.furnsh(fabsname)
self.loaded.append(fabsname)
[docs] def furnsh(self, fnames):
''' Furnsh a file, file list, or the name of the set
:param fnames: Single file, file list, or the name of the set.
Name of the set is eithr "full", "noatt", "min".
'''
try:
fnames = self.get(level=fnames)
except KeyError:
pass
if isinstance(fnames, str):
self.furnsheach(fnames)
else:
for fn in fnames:
try:
self.furnsheach(fn)
except spice.stypes.SpiceyError:
logger.warning("Failed to load {}. Usually no problem. Let's just continue.".format(fn))
__factory = MexSpiceFactory()
[docs]def init_full(update=True):
''' Furnsh default kernels.
:param update: No effect anymore.
'''
__factory.furnsh('full')
[docs]def init_noatt(update=True):
''' Furnsh default kernels.
:param update: No effect anymore.
'''
__factory.furnsh('noatt')
[docs]def init_min(update=True):
''' Furnsh minimum kernels.
:param update: No effect anymore.
'''
__factory.furnsh('min')
[docs]def init(level='full'):
corres = {'full': init_full, 'noatt': init_noatt, 'min': init_min}
corres[level]()
[docs]def furnsh(kernel_file):
''' Furnsh additional kernel file
'''
__factory.furnsh(kernel_file)
[docs]def updatemk(noatt=False, mkfile=None):
''' (DEPRECATED) Update the mk file.
Walk throught the path under ``spiceroot``,
create ``mk`` file.
Excluding ``*.txt``
>>> updatemk()
'''
import warnings
warnings.simplefilter('default')
warnings.warn('MK file is not used any more. Disregarded this call.',
category=DeprecationWarning)
[docs]def get_position(t, target='MEX', origin='MARS', frame='MSO', correction='LT+S'):
''' Get the position. By default, MEX in MSO.
>>> if not isdb():
... pass # doctest: +SKIP
... else:
... pos = get_position(datetime.datetime(2006, 12, 2, 15, 30))
... print('%.2f %.2f %.2f' % (pos[0], pos[1], pos[2])) # may change slightly
6265.30 -4041.08 -5409.83
'''
return get_posvel(t, target=target, origin=origin, frame=frame, correction=correction)[:3]
[docs]def height(t):
""" The height at the time of t.
:param t:
:return: The height above Mars (r=3396)
"""
p = get_position(t)
return np.sqrt((p ** 2).sum()) - 3396
[docs]def get_positions(tlist, target='MEX', origin='MARS', frame='MSO', correction='LT+S'):
''' Get the time series of the position.
:param tlist: Iterable of ``datetime.datetime`` object.
:returns: (N, 3) array of the position vector.
>>> tlist = [datetime.datetime(2012, 1, 1, 12), datetime.datetime(2012, 1, 1, 13)]
>>> positions = get_positions(tlist)
>>> print(positions.shape)
(2, 3)
>>> print("{p[0]:.2f} {p[1]:.2f} {p[2]:.2f}".format(p=positions[0, :]))
-1459.99 -4650.47 -7384.98
>>> print("{p[0]:.2f} {p[1]:.2f} {p[2]:.2f}".format(p=positions[1, :]))
-1461.35 3435.04 -1096.03
'''
poss = [get_position(t, target=target, origin=origin, frame=frame, correction=correction) for t in tlist]
return np.array(poss)
[docs]def state_generator(t0, t1=None, dt=datetime.timedelta(minutes=1), target='MEX', origin='MARS', frame='MSO', correction='LT+S'):
""" Return a generator to get state vector.
:keyword t1: Time to stop (Default is None, meaning 100*dt)
>>> init_noatt()
>>> t0 = datetime.datetime(2014, 10, 19, 18, 30)
>>> dt = datetime.timedelta(minutes=1)
>>> t1 = datetime.datetime(2014, 10, 19, 18, 33)
>>> for t, pv in state_generator(t0, dt=dt, t1=t1):
... print("{t} {pv[0]:.1f} {pv[1]:.1f} {pv[2]:.1f}".format(t=t, pv=pv))
2014-10-19 18:30:00 1833.9 2689.1 2343.6
2014-10-19 18:31:00 1927.6 2840.0 2178.1
2014-10-19 18:32:00 2016.9 2984.4 2007.5
2014-10-19 18:33:00 2101.7 3122.2 1832.6
"""
t = t0
if t1 is None:
t1 = t0 + 100 * dt
while t <= t1:
pv = get_posvel(t, target=target, origin=origin, frame=frame, correction=correction)
yield t, pv
t += dt
[docs]def pericenter(t0, dt=datetime.timedelta(seconds=1)):
""" Return the first pericenter around the time t0 from BC kernel.
>>> init_noatt()
>>> t0 = datetime.datetime(2014, 10, 19, 17)
>>> tp, hp = pericenter(t0) # The pericenter time after t0 from SPICE.
>>> print(tp, '{:.1f}'.format(hp))
2014-10-19 18:20:47 362.2
>>> tp, hp = pericenter(datetime.datetime(2014, 10, 18, 17)) # Given is ascending
>>> print(tp, '{:.1f}'.format(hp))
2014-10-18 14:22:52 362.2
See also :class:`irfpy.mexpvat.orbnum.OrbnumKernel`, but from Orbit file kernel.
"""
h0 = height(t0)
hmin = h0
tmin = t0
h1 = height(t0 + dt)
if h1 - h0 > 0: # Ascending
dt = -dt # Back in time
t1 = t0 + datetime.timedelta(hours=12)
for t, pv in state_generator(t0 + dt, t1=t1, dt=dt): # MEX MSO position
h1 = np.sqrt((pv[:3] ** 2).sum()) - 3396
if h1 < hmin:
hmin = h1
tmin = t
else:
return tmin, hmin
[docs]def when_at_height(t0, h, t1=None, dt=datetime.timedelta(seconds=1)):
""" Return the time after t0 when the altitude (from the surface) crosses the given altitude.
:param t0: Time to start calculation
:param h: Height in km.
:param end: Time to stop the investigate. Default is None, which gives 12 hours (to give more than 1 orbit).
:keywrod dt: Time resolution to calculate
:return: The time.
>>> init_noatt()
>>> t0 = datetime.datetime(2014, 10, 19, 16)
>>> t_at_1000 = when_at_height(t0, 1000) # Want to know the time of 1000 km altitude crossing
>>> print(t_at_1000)
2014-10-19 18:05:34
>>> p = get_position(datetime.datetime(2014, 10, 19, 18, 5, 34))
>>> print('{:.2f}'.format(np.sqrt((p ** 2).sum())-3396))
998.78
"""
if t1 is None:
t1 = t0 + datetime.timedelta(hours=12)
hc = h + 3396 # Add Martian radius
h0 = np.sqrt((get_posvel(t0)[:3] ** 2).sum()) # Height at the start
for t, pv in state_generator(t0 + dt, t1=t1, dt=dt):
h1 = np.sqrt((pv[:3] ** 2).sum())
if (hc - h0) * (hc - h1) <= 0:
return t
raise ex.IrfpyException('Height {} is did not reach between {} and {}'.format(h, t0, t1))
[docs]def get_posvel(t, target='MEX', origin='MARS', correction='LT+S', frame='MSO'):
et0 = spice.str2et(t.strftime('%FT%T.%f'))
posvel, lt = spice.spkezr(target, et0, frame, correction, origin)
return np.array(posvel)
[docs]def get_scannerangle(t):
''' Return the scanner angle in degrees.
'''
### The scanner angle is the angle between +Xsc and +Xsaf.
et0 = spice.str2et(t.strftime('%FT%T.%f'))
sc2saf = spice.pxform('MEX_SPACECRAFT', 'MEX_ASPERA_SAF', et0)
saf = spice.mxv(sc2saf, [1, 0, 0])
sdeg = np.rad2deg(np.arctan2(saf[1], saf[0]))
return sdeg
[docs]def et(t):
''' Return the et from the given datetime object.
'''
return spice.str2et(t.strftime('%FT%T.%f'))
[docs]def showstatus():
__factory.showstatus()
[docs]def loaded_kernels():
return __factory.loaded
[docs]def init_sidingspring():
""" Initialize the siding spring kernel. Also defines the name 'CSS'
>>> init() # doctest: +SKIP
>>> init_sidingspring() # doctest: +SKIP
>>> sspos = get_position(datetime.datetime(2014, 9, 18, 17, 30), target='CSS', origin='Mars', frame='MSO') # doctest: +SKIP
>>> print('{x[0]:.3e} {x[1]:.3e} {x[2]:.3e}'.format(x=sspos)) # doctest: +SKIP
4.079e+07 -1.221e+08 -7.222e+07
::
[mexpvat]
sidingspring_bsp = /path/to/siding_spring_s46.bsp
"""
rc = irfpy.util.irfpyrc.Rc(raises=True)
kernel = rc.get('mexpvat', 'sidingspring_bsp')
__factory.furnsh(kernel)
spice.boddef('CSS', 1003228)