Add origin to getPerturberState

moeyensj · Jun 12, 2020 · 2813794 · 2813794
1 parent b1561da
commit 2813794
Show file tree

Hide file tree

Showing 3 changed files with 24 additions and 12 deletions.
diff --git a/thor/observatories/state.py b/thor/observatories/state.py
@@ -5,7 +5,7 @@
 from ..constants import Constants as c
 from ..utils import _checkTime
 from ..utils import setupSPICE
-from ..orbits.ephemeris import getMajorBodyState
+from ..orbits.ephemeris import getPerturberState
 from .codes import readMPCObsCodeFile
 
 __all__ = ["getObserverState"]
@@ -81,15 +81,15 @@ def getObserverState(observatory_codes, observation_times):
         o_vec_ITRF93 = np.dot(R_EARTH, o_hat_ITRF93)
 
         # Grab earth state vector
-        state = getMajorBodyState("earth", observation_times)
+        state = getPerturberState("earth", observation_times, origin="heliocenter")
 
         # Convert MJD epochs in UTC to ET in TDB
         epochs_utc = observation_times.utc
         epochs_et = np.array([sp.str2et("JD {:.16f} UTC".format(i)) for i in epochs_utc.jd])
 
         # Grab rotaton matrices from ITRF93 to ecliptic J2000
         # The ITRF93 high accuracy Earth rotation model takes into account:
-        # Precession:  1976 IAU model due to Lieske.
+        # Precession:  1976 IAU model from Lieske.
         # Nutation:  1980 IAU model, with IERS corrections due to Herring et al.
         # True sidereal time using accurate values of TAI-UT1
         # Polar motion

diff --git a/thor/orbits/ephemeris/spice.py b/thor/orbits/ephemeris/spice.py
@@ -9,6 +9,8 @@
 S_TO_DAY = c.S_TO_DAY
 
 NAIF_MAPPING = {
+    "solar system barycenter" : 0,
+    "sun" : 10,
     "mercury" : 199,
     "venus" : 299,
     "earth" : 399,
@@ -19,12 +21,14 @@
     "neptune" : 899,
 }
 
-def getMajorBodyState(body_name, times):
+__all__ = ["getPerturberState"]
+
+def getPerturberState(body_name, times, origin="heliocenter"):
     """
-    Query the JPL ephemeris files loaded in SPICE for the state vectors of the desired
-    major planet. 
+    Query the JPL ephemeris files loaded in SPICE for the state vectors of desired perturbers.
     
-    Planets available: 
+    Major bodies and dynamical centers available: 
+        'solar system barycenter', 'sun',
         'mercury', 'venus', 'earth', 
         'mars', 'jupiter', 'saturn', 
         'uranus', 'neptune'
@@ -35,13 +39,23 @@ def getMajorBodyState(body_name, times):
         Name of major body.
     times : `~astropy.time.core.Time` (N)
         Times at which to get state vectors.
+    origin : {'barycenter', 'heliocenter'}
+        Coordinate system origin location.
     
     Returns
     -------
     states : `~numpy.ndarray` (N, 6)
         Heliocentric ecliptic J2000 state vector with postion in AU 
         and velocity in AU per day.
     """
+    if origin == "barycenter":
+        center = 0 # Solar System Barycenter
+    elif origin == "heliocenter":
+        center = 10 # Heliocenter
+    else: 
+        err = ("origin should be one of 'heliocenter' or 'barycenter'")
+        raise ValueError(err)
+
     # Make sure SPICE is ready to roll
     setupSPICE(verbose=False)
 
@@ -60,7 +74,7 @@ def getMajorBodyState(body_name, times):
             epoch, 
             'ECLIPJ2000',
             'NONE',
-            10 # Sun
+            center
         )
         states.append(state)
     states = np.vstack(states)

diff --git a/thor/utils/spice.py b/thor/utils/spice.py
@@ -49,7 +49,6 @@
 
 def getSPICEKernels(kernels=["LSK - Latest",  
                              "Planetary Constants", 
-                             "Earth PCK - Long Term Predict Low Accuracy",
                              "Earth PCK - Historical High Accuracy", 
                              "Earth PCK - Latest High Accuracy", 
                              "Planetary SPK"]):
@@ -83,17 +82,16 @@ def getSPICEKernels(kernels=["LSK - Latest",
     -------
     None
     """
-    for kernel, info in KERNELS.items():
+    for kernel in kernels:
         print("Checking for {} kernel...".format(kernel))
-        _downloadFile(os.path.join(os.path.dirname(__file__), "data", info[0]), info[1], update=info[2])
+        _downloadFile(os.path.join(os.path.dirname(__file__), "data", KERNELS[kernel][0]), KERNELS[kernel][1], update=KERNELS[kernel][2])
         print("")
     return
 
 def setupSPICE(
         kernels=[
             "LSK - Latest",  
             "Planetary Constants", 
-            "Earth PCK - Long Term Predict Low Accuracy",
             "Earth PCK - Historical High Accuracy", 
             "Earth PCK - Latest High Accuracy", 
             "Planetary SPK"