197 lines
7.8 KiB
Python
197 lines
7.8 KiB
Python
"""Compute positions from an ephemeris installed as a Python package.
|
|
|
|
Note: This entire module is DEPRECATED. The idea of distributing JPL
|
|
ephemerides as Python packages proved to be impractical (they were much
|
|
too large for the Python Package Index to easily store and distribute),
|
|
and it forced `jplephem` users get their ephemerides from a different
|
|
source than mainline astronomers, who use SPICE files. This package's
|
|
documentation now recommends avoiding this old code, and using the
|
|
features now build directly into the `SPK` class instead.
|
|
|
|
"""
|
|
import os
|
|
import numpy as np
|
|
|
|
|
|
class DateError(ValueError):
|
|
"""Date input is outside the range covered by the ephemeris."""
|
|
|
|
|
|
class Ephemeris(object):
|
|
"""[DEPRECATED] JPL planetary ephemeris for computing positions on dates."""
|
|
|
|
def __init__(self, module):
|
|
self.name = module.__name__.upper()
|
|
self.dirpath = os.path.dirname(module.__file__)
|
|
self.names = tuple(sorted(
|
|
name.split('-')[-1].split('.')[0]
|
|
for name in os.listdir(self.dirpath)
|
|
if not name.startswith('constants') and name.endswith('.npy')
|
|
))
|
|
path = self.path('constants.npy')
|
|
self.__dict__.update((k.decode('ascii'), v) for k, v in np.load(path))
|
|
self.earth_share = 1.0 / (1.0 + self.EMRAT)
|
|
self.moon_share = self.EMRAT / (1.0 + self.EMRAT)
|
|
self.sets = {}
|
|
|
|
def path(self, filename):
|
|
"""[DEPRECATED] Compute the path to a particular file in the ephemeris."""
|
|
return os.path.join(self.dirpath, filename)
|
|
|
|
def load(self, name):
|
|
"""[DEPRECATED] Load the polynomial series for `name` and return it."""
|
|
s = self.sets.get(name)
|
|
if s is None:
|
|
self.sets[name] = s = np.load(self.path('jpl-%s.npy' % name))
|
|
return s
|
|
|
|
def position(self, name, tdb, tdb2=0.0):
|
|
"""[DEPRECATED] Compute the position of `name` at time ``tdb [+ tdb2]``.
|
|
|
|
The position is returned as a NumPy array ``[x y z]``.
|
|
|
|
The barycentric dynamical time `tdb` argument should be a float.
|
|
If there are many dates you want computed, then make `tdb` an
|
|
array, which is more efficient than calling this method multiple
|
|
times; the return value will be a two-dimensional array giving a
|
|
row of values for each coordinate.
|
|
|
|
For extra precision, the time can be split into two floats; a
|
|
popular choice is to use `tdb` for the integer or half-integer
|
|
date, and `tdb2` to hold the remaining fraction.
|
|
|
|
Consult the `names` attribute of this ephemeris for the values
|
|
of `name` it supports, such as ``'mars'`` or ``'earthmoon'``.
|
|
|
|
"""
|
|
bundle = self.compute_bundle(name, tdb, tdb2)
|
|
return self.position_from_bundle(bundle)
|
|
|
|
def position_and_velocity(self, name, tdb, tdb2=0.0):
|
|
"""[DEPRECATED] Compute the position and velocity of `name` at ``tdb [+ tdb2]``.
|
|
|
|
The position and velocity are returned in a 2-tuple::
|
|
|
|
([x y z], [xdot ydot zdot])
|
|
|
|
The barycentric dynamical time `tdb` argument should be a float.
|
|
If there are many dates you want computed, then make `tdb` an
|
|
array, which is more efficient than calling this method multiple
|
|
times; the return values will be two-dimensional arrays giving a
|
|
row of values for each coordinate.
|
|
|
|
For extra precision, the time can be split into two floats; a
|
|
popular choice is to use `tdb` for the integer or half-integer
|
|
date, and `tdb2` to hold the remaining fraction.
|
|
|
|
Consult the `names` attribute of this ephemeris for the values
|
|
of `name` it supports, such as ``'mars'`` or ``'earthmoon'``.
|
|
|
|
"""
|
|
bundle = self.compute_bundle(name, tdb, tdb2)
|
|
position = self.position_from_bundle(bundle)
|
|
velocity = self.velocity_from_bundle(bundle)
|
|
return position, velocity
|
|
|
|
def compute(self, name, tdb):
|
|
"""[DEPRECATED] Legacy routine that concatenates position and velocity vectors.
|
|
|
|
This routine is deprecated. Use the methods `position()` and
|
|
`position_and_velocity()` instead. This method follows the same
|
|
calling convention, but incurs extra copy operations in order to
|
|
return a single NumPy array::
|
|
|
|
[x y z xdot ydot zdot]
|
|
|
|
"""
|
|
bundle = self.compute_bundle(name, tdb, 0.0)
|
|
position = self.position_from_bundle(bundle)
|
|
velocity = self.velocity_from_bundle(bundle)
|
|
return np.concatenate((position, velocity))
|
|
|
|
def compute_bundle(self, name, tdb, tdb2=0.0):
|
|
"""[DEPRECATED] Return a tuple of coefficients and parameters for `tdb`.
|
|
|
|
The return value is a tuple that bundles together the
|
|
coefficients and other Chebyshev intermediate values that are
|
|
needed for the computation of either the position or velocity.
|
|
The bundle can then be passed to either `position_from_bundle()`
|
|
or `velocity_from_bundle()` to finish the computation. See the
|
|
package-level documentation for details; most users will simply
|
|
call `position()` or `position_and_velocity()` instead.
|
|
|
|
The barycentric dynamical time `tdb` argument should be a float.
|
|
If there are many dates you want computed, then make `tdb` an
|
|
array, which is more efficient than calling this method multiple
|
|
times; the return values will be arrays providing a value for
|
|
each time in `tdb`.
|
|
|
|
For extra precision, the time can be split into two floats; a
|
|
popular choice is to use `tdb` for the integer or half-integer
|
|
date, and `tdb2` to hold the remaining fraction.
|
|
|
|
Consult the `names` attribute of this ephemeris for the values
|
|
of `name` it supports, such as ``'mars'`` or ``'earthmoon'``.
|
|
|
|
"""
|
|
input_was_scalar = getattr(tdb, 'shape', ()) == ()
|
|
if input_was_scalar:
|
|
tdb = np.array((tdb,))
|
|
# no need to deal with tdb2; numpy broadcast will add fine below.
|
|
|
|
coefficient_sets = self.load(name)
|
|
number_of_sets, axis_count, coefficient_count = coefficient_sets.shape
|
|
|
|
jalpha, jomega = self.jalpha, self.jomega
|
|
days_per_set = (jomega - jalpha) / number_of_sets
|
|
# to keep precision, first subtract, then add
|
|
index, offset = divmod((tdb - jalpha) + tdb2, days_per_set)
|
|
index = index.astype(int)
|
|
|
|
if (index < 0).any() or (number_of_sets < index).any():
|
|
raise DateError('ephemeris %s only covers dates %.1f through %.1f'
|
|
% (self.name, jalpha, jomega))
|
|
|
|
omegas = (index == number_of_sets)
|
|
index[omegas] -= 1
|
|
offset[omegas] += days_per_set
|
|
|
|
coefficients = np.rollaxis(coefficient_sets[index], 1)
|
|
|
|
# Chebyshev recurrence:
|
|
|
|
T = np.empty((coefficient_count, len(index)))
|
|
T[0] = 1.0
|
|
T[1] = t1 = 2.0 * offset / days_per_set - 1.0
|
|
twot1 = t1 + t1
|
|
for i in range(2, coefficient_count):
|
|
T[i] = twot1 * T[i-1] - T[i-2]
|
|
|
|
bundle = coefficients, days_per_set, T, twot1
|
|
return bundle
|
|
|
|
def position_from_bundle(self, bundle):
|
|
"""[DEPRECATED] Return position, given the `coefficient_bundle()` return value."""
|
|
|
|
coefficients, days_per_set, T, twot1 = bundle
|
|
return (T.T * coefficients).sum(axis=2)
|
|
|
|
def velocity_from_bundle(self, bundle):
|
|
"""[DEPRECATED] Return velocity, given the `coefficient_bundle()` return value."""
|
|
|
|
coefficients, days_per_set, T, twot1 = bundle
|
|
coefficient_count = coefficients.shape[2]
|
|
|
|
# Chebyshev derivative:
|
|
|
|
dT = np.empty_like(T)
|
|
dT[0] = 0.0
|
|
dT[1] = 1.0
|
|
dT[2] = twot1 + twot1
|
|
for i in range(3, coefficient_count):
|
|
dT[i] = twot1 * dT[i-1] - dT[i-2] + T[i-1] + T[i-1]
|
|
dT *= 2.0
|
|
dT /= days_per_set
|
|
|
|
return (dT.T * coefficients).sum(axis=2)
|