""" Misc tools related to SPICE
"""
import datetime
import numpy as np
from irfpy.util import utc
from irfpy.spice import spice as _spice
[docs]def set_error_action(strategy):
"""
:param strategy: Strategy, 'abort', 'report', 'return', 'ignore', or 'default'.
See SPICE document (``erract_c()``) for details.
"""
_spice.erract('set', 10, strategy)
[docs]def get_state(t, target, origin, frame, correction='NONE'):
""" Get the state of ``target`` from ``origin`` in ``frame``.
:param t: datetime object.
:param target: SPICE name of the target
:param origin: SPICE name of the origin
:param frame: SPICE name of the frame
:keyword correction: SPICE aberration correction
>>> get_state(datetime.datetime(2014, 10, 5, 3, 20, 15), 'MARS', "SUN", 'J2000', correction='LT+S') # doctest: +SKIP
array([ 8.39495688e+07, -1.75490667e+08, -8.27591341e+07,
2.31543409e+01, 1.08699922e+01, 4.36065716e+00])
"""
et0 = et(t)
posvel, lt = _spice.spkezr(target, et0, frame, correction, origin)
return posvel
[docs]def state(t0, t1, dt, target, origin, frame, correction='NONE'):
""" Generator of the state vector.
:param t0: Start time
:param t1: End time (exclusive)
:param dt: Time delta
:param target: SPICE name of target
:param origin: SPICE name of origin
:param frame: SPICE name of the frame
:param correction: SPICE aberration correction
:return: Generator, yielding the tuple (t, pos, vel).
>>> t0 = datetime.datetime(2014, 10, 1, 12, 0)
>>> t1 = datetime.datetime(2014, 10, 1, 12, 10)
>>> dt = datetime.timedelta(minutes=2)
>>> for t, xv in state(t0, t1, dt, 'SUN', 'MARS', 'J2000', correction='LT+S'): # doctest: +SKIP
... print(t, xv[:3], xv[3:])
2014-10-01 12:00:00 [ -7.66126240e+07 1.78785069e+08 8.40721063e+07] [-23.50762679 -10.09017858 -3.99343982]
2014-10-01 12:02:00 [ -7.66154449e+07 1.78783858e+08 8.40716271e+07] [-23.50749856 -10.09047781 -3.99358053]
2014-10-01 12:04:00 [ -7.66182658e+07 1.78782647e+08 8.40711478e+07] [-23.50737033 -10.09077705 -3.99372125]
2014-10-01 12:06:00 [ -7.66210866e+07 1.78781436e+08 8.40706686e+07] [-23.50724209 -10.09107629 -3.99386196]
2014-10-01 12:08:00 [ -7.66239075e+07 1.78780225e+08 8.40701893e+07] [-23.50711385 -10.09137552 -3.99400267]
"""
from irfpy.util import utc
for t in utc.trange(t0, t1, dt):
s = get_state(t, target, origin, frame, correction=correction)
yield (t, s)
[docs]def get_position(t, target, origin, frame, correction="NONE"):
""" Get the position of ``target`` from ``origin`` in ``frame``.
:param t: datetime object.
:param target: SPICE name of the target
:param origin: SPICE name of the origin
:param frame: SPICE name of the frame
:keyword correction: Correction allowed by SPICE
"""
et0 = et(t)
posvel, lt = _spice.spkezr(target, et0, frame, correction, origin)
return posvel[:3]
[docs]def position(t0, t1, dt, target, origin, frame, correction='NONE'):
""" Generator of the position vectors.
:param t0: Start time
:param t1: End time (exclusive)
:param dt: Time delta
:param target: SPICE name of target
:param origin: SPICE name of origin
:param frame: SPICE name of the frame
:param correction: SPICE aberration correction
:return: Generator, yielding the tuple (t, pos).
>>> t0 = datetime.datetime(2014, 10, 1, 12, 0)
>>> t1 = datetime.datetime(2014, 10, 1, 12, 10)
>>> dt = datetime.timedelta(minutes=2)
>>> for t, x in position(t0, t1, dt, 'SUN', 'MARS', 'J2000', correction='LT+S'): # doctest: +SKIP
... print(t, x)
2014-10-01 12:00:00 [ -7.66126240e+07 1.78785069e+08 8.40721063e+07]
2014-10-01 12:02:00 [ -7.66154449e+07 1.78783858e+08 8.40716271e+07]
2014-10-01 12:04:00 [ -7.66182658e+07 1.78782647e+08 8.40711478e+07]
2014-10-01 12:06:00 [ -7.66210866e+07 1.78781436e+08 8.40706686e+07]
2014-10-01 12:08:00 [ -7.66239075e+07 1.78780225e+08 8.40701893e+07]
"""
from irfpy.util import utc
for t in utc.trange(t0, t1, dt):
s = get_position(t, target, origin, frame, correction=correction)
x = s[:3]
yield (t, x)
[docs]def get_velocity(t, target, origin, frame, correction="NONE"):
""" Get the velocity of ``target`` relative to ``origin`` in ``frame``.
:param t: datetime object.
:param target: SPICE name of the target
:param origin: SPICE name of the origin
:param frame: SPICE name of the frame
:keyword correction: Correction allowed by SPICE
"""
et0 = et(t)
posvel, lt = _spice.spkezr(target, et0, frame, correction, origin)
return posvel[3:]
[docs]def velocity(t0, t1, dt, target, origin, frame, correction='NONE'):
""" Generator of the velocity vectors.
:param t0: Start time
:param t1: End time (exclusive)
:param dt: Time delta
:param target: SPICE name of target
:param origin: SPICE name of origin
:param frame: SPICE name of the frame
:param correction: SPICE aberration correction
:return: Generator, yielding the tuple (t, vel).
>>> t0 = datetime.datetime(2014, 10, 1, 12, 0)
>>> t1 = datetime.datetime(2014, 10, 1, 12, 10)
>>> dt = datetime.timedelta(minutes=2)
>>> for t, v in velocity(t0, t1, dt, 'SUN', 'MARS', 'J2000', correction='LT+S'): # doctest: +SKIP
... print(t, v)
2014-10-01 12:00:00 [-23.50762679 -10.09017858 -3.99343982]
2014-10-01 12:02:00 [-23.50749856 -10.09047781 -3.99358053]
2014-10-01 12:04:00 [-23.50737033 -10.09077705 -3.99372125]
2014-10-01 12:06:00 [-23.50724209 -10.09107629 -3.99386196]
2014-10-01 12:08:00 [-23.50711385 -10.09137552 -3.99400267]
"""
for t in utc.trange(t0, t1, dt):
s = get_velocity(t, target, origin, frame, correction=correction)
yield (t, s)
[docs]def get_distance(t, target, origin, frame='J2000', correction='NONE'):
""" Return the distance.
"""
pos = get_position(t, target, origin, frame, correction=correction)
return np.linalg.norm(pos)
[docs]def distance(t0, t1, dt, target, origin, frame='J2000', correction='NONE'):
""" Generator for distance.
>>> for t, h in distance(datetime.datetime(2014, 10, 1, 15),
... datetime.datetime(2014, 10, 1, 16), datetime.timedelta(minutes=1),
... 'MEX', 'MARS', frame='MSO'): # doctest: +SKIP
... print(t, h)
"""
for t in utc.trange(t0, t1, dt):
d = get_distance(t, target, origin, frame=frame, correction=correction)
yield (t, d)
[docs]def get_sza(t, target, origin, frame='J2000', correction='NONE'):
""" Return the solar zenith angle in degrees.
"""
post = get_position(t, target, origin, frame, correction=correction) # (3,) array
poss = get_position(t, 'SUN', origin, frame, correction=correction)
vpost = post / np.linalg.norm(post)
vposs = poss / np.linalg.norm(poss)
inner = np.clip((vpost * vposs).sum(), -1, 1) # Rare case, slightly above 1 or less than 1 is possible due to the numerical error.
angle = np.arccos(inner)
return np.rad2deg(angle)
[docs]def sza(t0, t1, dt, target, origin, frame='J2000', correction='NONE'):
for t in utc.trange(t0, t1, dt):
s = get_sza(t, target, origin, frame=frame, correction=correction)
yield (t, s)
[docs]def get_conversion_matrix(t, frame0, frame1):
""" Return the frame conversion matrix (3x3)
If you want to convert from the MEX_SPACECRAFT frame to MSO frame,
you may get a conversion matrix as follows.
>>> m = get_conversion_matrix(datetime.datetime(2014, 1, 23, 15, 30), 'MEX_SPACECRAFT', 'MSO') # doctest: +SKIP
>>> print(m.shape) # doctest: +SKIP
(3, 3)
"""
et0 = et(t)
m = _spice.pxform(frame0, frame1, et0)
return np.array(m)
[docs]def conversion_matrix(t0, t1, dt, frame0, frame1):
for t in utc.trange(t0, t1, dt):
yield (t, get_conversion_matrix(t, frame0, frame1))
[docs]def get_conversion_matrix_6d(t, frame0, frame1):
""" Return the frame conversion matrix (6x6)
If you want to convert from the J2000 frame to IAU_MARS frame,
you may get a conversion matrix as follows.
>>> m = get_conversion_matrix(datetime.datetime(2014, 1, 23, 15, 30), 'J2000', 'IAU_MARS') # doctest: +SKIP
>>> print(m.shape) # doctest: +SKIP
(6, 6)
"""
et0 = et(t)
m = _spice.sxform(frame0, frame1, et0)
return np.array(m)
[docs]def conversion_matrix_6d(t0, t1, dt, frame0, frame1):
for t in utc.trange(t0, t1, dt):
yield (t, get_conversion_matrix_6d(t, frame0, frame1))
[docs]def et(t):
""" Get the ``et`` from the datetime object.
Preparation.
>>> import datetime
>>> from irfpy.spice import localize
>>> localize.furnsh_tls()
Execution.
>>> print(et(datetime.datetime(2020, 1, 1)))
631108869.1839073
>>> print(et([datetime.datetime(2020, 1, 1), datetime.datetime(2020, 1, 2)]))
(631108869.1839073, 631195269.1839362)
"""
try:
et = _spice.str2et(t.strftime('%FT%T'))
except AttributeError:
et = tuple([_spice.str2et(_t.strftime('%FT%T')) for _t in t])
return et