Merge pull request #33 from moeyensj/gauss

Gauss

MANIFEST.in

@@ -0,0 +1,7 @@
+include thor/*.py
+include *.md
+include *.toml
+exclude .travis.yaml
+exclude requirements_travis.txt
+exclude azure-pipelines.yml 
+exclude Dockerfile

README.md

@@ -4,8 +4,9 @@ Tracklet-less Heliocentric Orbit Recovery
 [![Build Status](https://dev.azure.com/moeyensj/thor/_apis/build/status/moeyensj.thor?branchName=master)](https://dev.azure.com/moeyensj/thor/_build/latest?definitionId=2&branchName=master)
 [![Build Status](https://www.travis-ci.com/moeyensj/thor.svg?token=sWjpnqPgpHyuq3j7qPuj&branch=master)](https://www.travis-ci.com/moeyensj/thor)
 [![Coverage Status](https://coveralls.io/repos/github/moeyensj/thor/badge.svg?branch=master&t=pdSkQA)](https://coveralls.io/github/moeyensj/thor?branch=master)
-[![Docker Pulls](https://img.shields.io/docker/pulls/moeyensj/thor)](https://hub.docker.com/r/moeyensj/thor)
-[![License](https://img.shields.io/badge/License-BSD%203--Clause-blue.svg)](https://opensource.org/licenses/BSD-3-Clause)
+[![Docker Pulls](https://img.shields.io/docker/pulls/moeyensj/thor)](https://hub.docker.com/r/moeyensj/thor)  
+[![Python 3.6](https://img.shields.io/badge/Python-3.6%2B-blue)](https://img.shields.io/badge/Python-3.6%2B-blue)
+[![License](https://img.shields.io/badge/License-BSD%203--Clause-blue.svg)](https://opensource.org/licenses/BSD-3-Clause)  
 
 ## Installation
 

pyproject.toml

@@ -0,0 +1,5 @@
+[build-system]
+requires = ["setuptools>=42", "wheel", "setuptools_scm[toml]>=3.4"]
+
+[tool.setuptools_scm]
+write_to = "pkg/version.py"

requirements.txt

@@ -15,4 +15,7 @@ numba
 spiceypy
 pytest
 pytest-cov
-coveralls
+coveralls
+setuptools>=42
+wheel
+setuptools_scm>=3.4

requirements_travis.txt

@@ -14,4 +14,7 @@ numba
 spiceypy
 pytest
 pytest-cov<2.6.0
-coveralls
+coveralls
+setuptools>=42
+wheel
+setuptools_scm>=3.4

setup.py

@@ -2,16 +2,16 @@
 
 setup(
    name="thor",
-   version="0.0.1.dev0",
-   description="Tracklet-less Heliocentric Orbit Recovery",
    license="BSD 3-Clause License",
    author="Joachim Moeyens, Mario Juric",
    author_email="moeyensj@uw.edu",
+   long_description=open("README.md").read(),
+   long_description_content_type="text/markdown",
    url="https://github.com/moeyensj/thor",
    packages=["thor"],
    package_dir={"thor": "thor"},
-   package_data={"thor": ["data/*.orb",
-                             "tests/data/*"]},
-   setup_requires=["pytest-runner"],
+   package_data={"thor": ["tests/*.txt"]},
+   use_scm_version=True,
+   setup_requires=["pytest-runner", "setuptools_scm"],
    tests_require=["pytest"],
 )

thor/constants.py

@@ -9,16 +9,16 @@ class Constants:
     # Speed of light: AU per day (173.14463267424034)
     C = c.c.to(u.au / u.d).value 
 
-    # Gravitational constant: AU**3 / M_sun / d**2 (0.00029591220819207784)
-    G = c.G.to(u.AU**3 / u.M_sun / u.d**2).value 
+    # Gravitational constant:  AU**3 / M_sun / d**2 (0.295912208285591100E-3 -- DE431/DE430)
+    G = 0.295912208285591100E-3
 
     # Solar Mass: M_sun (1.0)
     M_SUN = 1.0
 
     # Earth Mass: M_sun (3.0034893488507934e-06)
     M_EARTH = u.M_earth.to(u.M_sun)
 
-    # Earth Equatorial Radius (6378.1363 km (DE431/DE430))
+    # Earth Equatorial Radius: km (6378.1363 -- DE431/DE430)
     R_EARTH = (6378.1363 * u.km).to(u.AU)
 
     # Mean Obliquity at J2000: radians (0.40909280422232897)

thor/orbits/iod/gauss.py

@@ -71,7 +71,10 @@ def _calcFG(r2_mag, t32, t21, mu=MU):
 
 def _calcM(r0_mag, r_mag, f, g, f_dot, g_dot, c0, c1, c2, c3, c4, c5, alpha, chi, mu=MU):
     # Universal variables will differ between different texts and works in the literature.
-    # c0, c1, c2, c3, c4, c5 are expected to be 
+    # c0, c1, c2, c3, c4, c5 are expected to follow the Battin formalism (adopted by both
+    # Vallado and Curtis in their books). The M matrix is proposed by Shepperd 1985 and follows
+    # the Goodyear formalism. Conversions between the two formalisms can be derived from Table 1 in
+    # Everhart & Pitkin 1982.
     w = chi / np.sqrt(mu)
     alpha_alt = - mu * alpha
     U0 = (1 - alpha_alt * chi**2) * c0
@@ -84,7 +87,7 @@ def _calcM(r0_mag, r_mag, f, g, f_dot, g_dot, c0, c1, c2, c3, c4, c5, alpha, chi
     F = f_dot
     G = g_dot
 
-    # Equations 18 and 19 in Sheppard 1985
+    # Equations 18 and 19 in Shepperd 1985
     U = (U2 * U3 + w * U4 - 3 * U5) / 3
     W = g * U2 + 3 * mu * U
 
@@ -124,7 +127,7 @@ def _calcStateTransitionMatrix(M, r0, v0, f, g, f_dot, g_dot, r, v):
     ])
     return phi
 
-def _iterateGaussIOD(orbit, t21, t32, q1, q2, q3, rho1, rho2, rho3, mu=MU, max_iter=10, tol=1e-15):
+def _iterateGaussIOD(orbit, t21, t32, q1, q2, q3, rho1, rho2, rho3, light_time=True, mu=MU, max_iter=10, tol=1e-15):
     # Iterate over the polynomial solution from Gauss using the universal anomaly 
     # formalism until the solution converges or the maximum number of iterations is reached
 
@@ -159,6 +162,12 @@ def _iterateGaussIOD(orbit, t21, t32, q1, q2, q3, rho1, rho2, rho3, mu=MU, max_i
         # then calculate the Lagrange coefficients and the state for each observation.
         # Use those to calculate the state transition matrix
         for j, dt in enumerate([-t21, t32]):
+            if light_time is True:
+                if j == 1:
+                    dt += (rho2_mag - rho1_mag) / C
+                else:
+                    dt -= (rho3_mag - rho2_mag) / C
+
             # Calculate the universal anomaly 
             # Universal anomaly here is defined in such a way that it satisfies the following
             # differential equation:
@@ -286,7 +295,7 @@ def calcGauss(r1, r2, r3, t1, t2, t3):
     f1, g1, f3, g3 = _calcFG(r2_mag, t32, t21)
     return (1 / (f1 * g3 - f3 * g1)) * (-f3 * r1 + f1 * r3)
 
-def gaussIOD(coords_eq_ang, t, coords_obs, velocity_method="gibbs", iterate=True, mu=MU, max_iter=10, tol=1e-15):
+def gaussIOD(coords_eq_ang, t, coords_obs, velocity_method="gibbs", light_time=True, iterate=True, mu=MU, max_iter=10, tol=1e-15):
     """
     Compute up to three intial orbits using three observations in angular equatorial
     coordinates. 
@@ -302,6 +311,9 @@ def gaussIOD(coords_eq_ang, t, coords_obs, velocity_method="gibbs", iterate=True
     velocity_method : {'gauss', gibbs', 'herrick+gibbs'}, optional
         Which method to use for calculating the velocity at the second observation.
         [Default = 'gibbs']
+    light_time : bool, optional
+        Correct for light travel time. 
+        [Default = True]
     iterate : bool, optional
         Iterate initial orbit using universal anomaly to better approximate the 
         Lagrange coefficients. 
@@ -346,7 +358,7 @@ def gaussIOD(coords_eq_ang, t, coords_obs, velocity_method="gibbs", iterate=True
     coseps2 = np.dot(q2, rho2_hat) / q2_mag
     C0 = V * t31 * q2_mag**4 / B
     h0 = - A / B
-    
+
     if np.isnan(C0) or np.isnan(h0):
         return np.array([])
 
@@ -376,7 +388,7 @@ def gaussIOD(coords_eq_ang, t, coords_obs, velocity_method="gibbs", iterate=True
 
         # Test if we get the same rho2 as using equation 22 in Milani et al. 2008
         rho2_mag = (h0 - q2_mag**3 / r2_mag**3) * q2_mag / C0
-        np.testing.assert_almost_equal(np.dot(rho2_mag, rho2_hat), rho2)
+        #np.testing.assert_almost_equal(np.dot(rho2_mag, rho2_hat), rho2)
 
         r1 = q1 + rho1
         r2 = q2 + rho2
@@ -390,16 +402,30 @@ def gaussIOD(coords_eq_ang, t, coords_obs, velocity_method="gibbs", iterate=True
             v2 = calcHerrickGibbs(r1, r2, r3, t1, t2, t3)
         else:
             raise ValueError("velocity_method should be one of {'gauss', 'gibbs', 'herrick+gibbs'}")
+
+        epoch = t2
         orbit = np.concatenate([r2, v2])
 
         if iterate == True:
             orbit = _iterateGaussIOD(orbit, t21, t32, 
                                      q1, q2, q3, 
                                      rho1, rho2, rho3,
-                                     mu=mu, max_iter=max_iter, tol=tol)
+                                     light_time=light_time,
+                                     mu=mu, 
+                                     max_iter=max_iter,
+                                     tol=tol)
+
+        if light_time == True:
+            rho2_mag = np.linalg.norm(orbit[:3] - q2)
+            lt = rho2_mag / C
+            epoch -= lt
 
         if np.linalg.norm(orbit[3:]) >= C:
             print("Velocity is greater than speed of light!")
-        orbits.append(orbit)
+
+        if (np.linalg.norm(orbit[:3]) > 300.) and (np.linalg.norm(orbit[3:]) > 25.):
+            continue
+
+        orbits.append(np.hstack([epoch, orbit]))
 
     return np.array(orbits)

thor/orbits/iod/iod.py

@@ -1,46 +1,165 @@
+import sys
 import numpy as np
+from itertools import combinations
 
 from ...config import Config
+from ...constants import Constants as c
+from .gauss import gaussIOD
+from ..propagate import propagateOrbits
 
-__all__ = ["selectObservations"]
+__all__ = ["selectObservations",
+           "iod"]
 
-def selectObservations(observations, method="first+middle+last", columnMapping=Config.columnMapping):
+MU = c.G * c.M_SUN
+
+
+def selectObservations(observations, method="combinations", columnMapping=Config.columnMapping):
     """
     Selects which three observations to use for IOD depending on the method. 
     
     Methods:
         'first+middle+last' : Grab the first, middle and last observations in time. 
         'thirds' : Grab the middle observation in the first third, second third, and final third. 
+        'combinations' : Return the observation IDs corresponding to every possible combination of three observations with
+            non-coinciding observation times.
     
     Parameters
     ----------
     observations : `~pandas.DataFrame`
         Pandas DataFrame containing observations with at least a column of observation IDs and a column
         of exposure times. 
-    method : {'first+middle+last', 'thirds}, optional
+    method : {'first+middle+last', 'thirds', 'combinations'}, optional
         Which method to use to select observations. 
-        [Default = `first+middle+last`]
+        [Default = 'combinations']
     columnMapping : dict, optional
         Column name mapping of observations to internally used column names. 
         [Default = `~thor.Config.columnMapping`]
         
     Returns
     -------
-    obs_id : `~numpy.ndarray' (3 or 0)
+    obs_id : `~numpy.ndarray' (N, 3 or 0)
         An array of selected observation IDs. If three unique observations could 
-        not be selected than returns an empty array. 
+        not be selected then returns an empty array. 
     """
+    obs_ids = observations[columnMapping["obs_id"]].values
+    if len(obs_ids) < 3:
+        return np.array([])
+
+    indexes = np.arange(0, len(obs_ids))
+    times = observations[columnMapping["exp_mjd"]].values
+    selected = np.array([])
+
     if method == "first+middle+last":
-        selected = np.percentile(observations[columnMapping["exp_mjd"]].values, [0, 50, 100], interpolation="nearest")
+        selected_times = np.percentile(times, 
+                                 [0, 50, 100], 
+                                 interpolation="nearest")
+        selected_index = np.intersect1d(times, selected_times, return_indices=True)[1]
+        selected_index = np.array([selected_index])
+
     elif method == "thirds":
-        selected = np.percentile(observations[columnMapping["exp_mjd"]].values, [1/6 * 100, 50, 5/6*100], interpolation="nearest")
+        selected_times = np.percentile(times, 
+                                 [1/6*100, 50, 5/6*100], 
+                                 interpolation="nearest")
+        selected_index = np.intersect1d(times, selected_times, return_indices=True)[1]
+        selected_index = np.array([selected_index])
+
+    elif method == "combinations":
+        # Make all possible combinations of 3 observations
+        selected_index = np.array([np.array(index) for index in combinations(indexes, 3)])
+
     else:
         raise ValueError("method should be one of {'first+middle+last', 'thirds'}")
-
-    if len(np.unique(selected)) != 3:
-        print("Could not find three observations that satisfy the criteria.")
+
+    # Make sure each returned combination of observation ids have at least 3 unique
+    # times 
+    keep = []
+    for i, comb in enumerate(times[selected_index]):
+        if len(np.unique(comb)) == 3:
+            keep.append(i)
+    keep = np.array(keep)
+
+    # Return an empty array if no observations satisfy the criteria
+    if len(keep) == 0:
         return np.array([])
 
-    index = np.intersect1d(observations[columnMapping["exp_mjd"]].values, selected, return_indices=True)[1]
-    return observations[columnMapping["obs_id"]].values[index]
-
+    return obs_ids[selected_index[keep, :]]
+
+
+def iod(observations,
+        observation_selection_method="combinations",
+        iterate=True, 
+        light_time=True,
+        max_iter=50, 
+        tol=1e-15, 
+        propagatorKwargs={
+            "observatoryCode" : "I11",
+            "mjdScale" : "UTC",
+            "dynamical_model" : "2",
+        },
+        mu=MU, 
+        columnMapping=Config.columnMapping):
+
+    # Extract column names
+    obs_id_col = columnMapping["obs_id"]
+    ra_col = columnMapping["RA_deg"]
+    dec_col = columnMapping["Dec_deg"]
+    time_col = columnMapping["exp_mjd"]
+    ra_err_col = columnMapping["RA_sigma_deg"]
+    dec_err_col = columnMapping["Dec_sigma_deg"]
+    obs_x_col = columnMapping["obs_x_au"]
+    obs_y_col = columnMapping["obs_y_au"]
+    obs_z_col = columnMapping["obs_z_au"]
+
+    # Extract observation IDs, sky-plane positions, sky-plane position uncertainties, times of observation,
+    # and the location of the observer at each time
+    obs_ids_all = observations[obs_id_col].values
+    coords_eq_ang_all = observations[observations[obs_id_col].isin(obs_ids_all)][[ra_col, dec_col]].values
+    coords_eq_ang_err_all = observations[observations[obs_id_col].isin(obs_ids_all)][[ra_err_col, dec_err_col]].values
+    coords_obs_all = observations[observations[obs_id_col].isin(obs_ids_all)][[obs_x_col, obs_y_col, obs_z_col]].values
+    times_all = observations[observations[obs_id_col].isin(obs_ids_all)][time_col].values
+
+    # Select observation IDs to use for IOD
+    obs_ids = selectObservations(observations, method=observation_selection_method, columnMapping=columnMapping)
+
+    min_chi2 = 1e10
+    best_orbit = None
+    best_obs_ids = None
+
+    for ids in obs_ids:
+        # Grab sky-plane positions of the selected observations, the heliocentric ecliptic position of the observer,
+        # and the times at which the observations occur
+        coords_eq_ang = observations[observations[obs_id_col].isin(ids)][[ra_col, dec_col]].values
+        coords_obs = observations[observations[obs_id_col].isin(ids)][[obs_x_col, obs_y_col, obs_z_col]].values
+        times = observations[observations[obs_id_col].isin(ids)][time_col].values
+
+        # Run IOD 
+        orbits_iod = gaussIOD(coords_eq_ang, times, coords_obs, light_time=light_time, iterate=iterate, max_iter=max_iter, tol=tol)
+        if np.all(np.isnan(orbits_iod)) == True:
+            continue
+
+        # Propagate initial orbit to all observation times
+        orbits = propagateOrbits(orbits_iod[:, 1:], orbits_iod[:, 0], times_all, **propagatorKwargs)
+        orbits = orbits[['orbit_id', 'mjd', 'RA_deg', 'Dec_deg', 
+                         'HEclObj_X_au', 'HEclObj_Y_au', 'HEclObj_Z_au',
+                         'HEclObj_dX/dt_au_p_day', 'HEclObj_dY/dt_au_p_day', 'HEclObj_dZ/dt_au_p_day']].values
+
+        # For each unique initial orbit calculate residuals and chi-squared
+        # Find the orbit which yields the lowest chi-squared
+        orbit_ids = np.unique(orbits[:, 0])
+        for i, orbit_id in enumerate(orbit_ids):
+            orbit = orbits[np.where(orbits[:, 0] == orbit_id)]
+
+            pred_dec = np.radians(orbit[:, 3])
+            residual_ra = (coords_eq_ang_all[:, 0] - orbit[:, 2]) * np.cos(pred_dec)
+            residual_dec =  (coords_eq_ang_all[:, 1] - orbit[:, 3])
+
+            chi2 = np.sum(residual_ra**2 / coords_eq_ang_err_all[:, 0]**2 + residual_dec**2 / coords_eq_ang_err_all[:, 1]**2) / (2 * len(residual_ra) - 6)
+
+            if chi2 < min_chi2:
+                best_orbit = orbits_iod[i, :]
+                best_obs_ids = ids
+                min_chi2 = chi2
+
+    return best_orbit, best_obs_ids, min_chi2
+
+