From 6ea360ce0151cf6cbc236a8a8283ce22cdaeebdb Mon Sep 17 00:00:00 2001 From: Joachim Moeyens Date: Wed, 18 Dec 2019 08:02:21 +0000 Subject: [PATCH 01/16] Add formalism clarification comment --- thor/orbits/iod/gauss.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/thor/orbits/iod/gauss.py b/thor/orbits/iod/gauss.py index 8285d0eb..34d50dbe 100644 --- a/thor/orbits/iod/gauss.py +++ b/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 From 83a888ec676982dabd3dacc90dfc44a5b6cac5e4 Mon Sep 17 00:00:00 2001 From: Joachim Moeyens Date: Tue, 14 Jan 2020 11:07:58 -0800 Subject: [PATCH 02/16] Add combinations option to selectObservations --- thor/orbits/iod/iod.py | 57 ++++++++++++++++++++++++++++++++---------- 1 file changed, 44 insertions(+), 13 deletions(-) diff --git a/thor/orbits/iod/iod.py b/thor/orbits/iod/iod.py index 89e1afd9..f88fc7e1 100644 --- a/thor/orbits/iod/iod.py +++ b/thor/orbits/iod/iod.py @@ -1,46 +1,77 @@ import numpy as np +from itertools import combinations from ...config import Config __all__ = ["selectObservations"] -def selectObservations(observations, method="first+middle+last", columnMapping=Config.columnMapping): + +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 + 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] - \ No newline at end of file + return obs_ids[selected_index[keep, :]] + + \ No newline at end of file From fc6dae325a181decaa6d4b93539cb4cea23be2da Mon Sep 17 00:00:00 2001 From: Joachim Moeyens Date: Wed, 15 Jan 2020 12:28:06 -0800 Subject: [PATCH 03/16] Add light time correction to output orbits --- thor/orbits/iod/gauss.py | 27 +++++++++++++++++++++++---- 1 file changed, 23 insertions(+), 4 deletions(-) diff --git a/thor/orbits/iod/gauss.py b/thor/orbits/iod/gauss.py index 34d50dbe..d7eed35b 100644 --- a/thor/orbits/iod/gauss.py +++ b/thor/orbits/iod/gauss.py @@ -127,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 @@ -162,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: @@ -289,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. @@ -305,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. @@ -393,16 +402,26 @@ 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) + orbits.append(np.hstack([epoch, orbit])) return np.array(orbits) \ No newline at end of file From 0afb2e3e3a37381af7c64ee40faaddeea68210d8 Mon Sep 17 00:00:00 2001 From: Joachim Moeyens Date: Wed, 15 Jan 2020 12:28:40 -0800 Subject: [PATCH 04/16] Add propagation dynamical model type to propagateOrbits --- thor/orbits/propagate/pyoorb.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/thor/orbits/propagate/pyoorb.py b/thor/orbits/propagate/pyoorb.py index e3e5905d..acbcfca6 100644 --- a/thor/orbits/propagate/pyoorb.py +++ b/thor/orbits/propagate/pyoorb.py @@ -16,6 +16,7 @@ def propagateOrbits(elements, G=0.15, M1=1, K1=1, + dynamical_model="2", observatoryCode=Config.oorbObservatoryCode): """ Propagate a test particle using its ecliptic coordinates and velocity to @@ -43,7 +44,10 @@ def propagateOrbits(elements, M1 : float or `~np.ndarray` (N), optional K1 : float or `~np.ndarray` (N), optional - + + dynamical_model : str, optional + Propagate orbits using 2-body or n-body integration. + [Default = "2"] observatoryCode : str, optional Observatory from which to measure ephemerides. [Default = `~thor.Config.oorbObservatoryCode`] @@ -113,7 +117,7 @@ def propagateOrbits(elements, ephemeris, err = oo.pyoorb.oorb_ephemeris_full(in_orbits=orbits, in_obscode=observatoryCode, in_date_ephems=epochs, - in_dynmodel='2', + in_dynmodel=dynamical_model, ) columns = ["mjd", "RA_deg", From c7387b5bdccfe3d3d332efd83633722573af1a19 Mon Sep 17 00:00:00 2001 From: Joachim Moeyens Date: Wed, 15 Jan 2020 12:29:18 -0800 Subject: [PATCH 05/16] Add iod function --- thor/orbits/iod/iod.py | 84 +++++++++++++++++++++++++++++++++++++++++- 1 file changed, 83 insertions(+), 1 deletion(-) diff --git a/thor/orbits/iod/iod.py b/thor/orbits/iod/iod.py index f88fc7e1..79aec19d 100644 --- a/thor/orbits/iod/iod.py +++ b/thor/orbits/iod/iod.py @@ -2,8 +2,14 @@ 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"] + +MU = c.G * c.M_SUN def selectObservations(observations, method="combinations", columnMapping=Config.columnMapping): @@ -73,5 +79,81 @@ def selectObservations(observations, method="combinations", columnMapping=Config return np.array([]) 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) + + # 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 \ No newline at end of file From c6e91d92651884231ad0a5ba7b9d8f8b2410d233 Mon Sep 17 00:00:00 2001 From: Joachim Moeyens Date: Wed, 15 Jan 2020 12:56:28 -0800 Subject: [PATCH 06/16] Comment out equality test --- thor/orbits/iod/gauss.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/thor/orbits/iod/gauss.py b/thor/orbits/iod/gauss.py index d7eed35b..ace26194 100644 --- a/thor/orbits/iod/gauss.py +++ b/thor/orbits/iod/gauss.py @@ -388,7 +388,7 @@ def gaussIOD(coords_eq_ang, t, coords_obs, velocity_method="gibbs", light_time=T # 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 From e5fdba0cae732cb3310104f69e67801b3cbf52fc Mon Sep 17 00:00:00 2001 From: Joachim Moeyens Date: Wed, 15 Jan 2020 18:11:44 -0800 Subject: [PATCH 07/16] Fix corner case with not enough observations --- thor/orbits/iod/iod.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/thor/orbits/iod/iod.py b/thor/orbits/iod/iod.py index 79aec19d..398c6e63 100644 --- a/thor/orbits/iod/iod.py +++ b/thor/orbits/iod/iod.py @@ -41,6 +41,9 @@ def selectObservations(observations, method="combinations", columnMapping=Config 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([]) From b26e67499cf56334fabb8239436d91c2703ebb02 Mon Sep 17 00:00:00 2001 From: Joachim Moeyens Date: Thu, 16 Jan 2020 14:12:15 -0800 Subject: [PATCH 08/16] Throwaway orbits more than 300 AU away with velocites faster than 25 AU per day: pyoorb crashes --- thor/orbits/iod/gauss.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/thor/orbits/iod/gauss.py b/thor/orbits/iod/gauss.py index ace26194..2805787e 100644 --- a/thor/orbits/iod/gauss.py +++ b/thor/orbits/iod/gauss.py @@ -358,7 +358,7 @@ def gaussIOD(coords_eq_ang, t, coords_obs, velocity_method="gibbs", light_time=T 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([]) @@ -422,6 +422,10 @@ def gaussIOD(coords_eq_ang, t, coords_obs, velocity_method="gibbs", light_time=T if np.linalg.norm(orbit[3:]) >= C: print("Velocity is greater than speed of light!") + + 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) \ No newline at end of file From 378b19c4a88192ef813776abe921323d632cb88f Mon Sep 17 00:00:00 2001 From: Joachim Moeyens Date: Thu, 16 Jan 2020 14:13:46 -0800 Subject: [PATCH 09/16] Add propagation failure catches to IOD function --- thor/orbits/iod/iod.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/thor/orbits/iod/iod.py b/thor/orbits/iod/iod.py index 398c6e63..2c1fafc9 100644 --- a/thor/orbits/iod/iod.py +++ b/thor/orbits/iod/iod.py @@ -1,3 +1,4 @@ +import sys import numpy as np from itertools import combinations @@ -133,9 +134,13 @@ def iod(observations, # 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) + if np.all(orbits.values) == 0.0: + continue 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 From bf1f13dc189011b58282bf250cc7f4489a14dd4e Mon Sep 17 00:00:00 2001 From: Joachim Moeyens Date: Fri, 17 Jan 2020 09:49:41 -0800 Subject: [PATCH 10/16] Remove 0.0 check after propagate --- thor/orbits/iod/iod.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/thor/orbits/iod/iod.py b/thor/orbits/iod/iod.py index 2c1fafc9..bb841b19 100644 --- a/thor/orbits/iod/iod.py +++ b/thor/orbits/iod/iod.py @@ -139,8 +139,6 @@ def iod(observations, # Propagate initial orbit to all observation times orbits = propagateOrbits(orbits_iod[:, 1:], orbits_iod[:, 0], times_all, **propagatorKwargs) - if np.all(orbits.values) == 0.0: - continue 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 From f017196650ccbabf3f88ab32bc32502a4ff8a3f2 Mon Sep 17 00:00:00 2001 From: Joachim Moeyens Date: Mon, 27 Jan 2020 15:16:41 -0800 Subject: [PATCH 11/16] Add Python 3.6+ badge to README --- README.md | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 3fdbed06..0466babe 100644 --- a/README.md +++ b/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 From 1b3c421f4125673a88f0f7d40bf4284149acede3 Mon Sep 17 00:00:00 2001 From: Joachim Moeyens Date: Wed, 5 Feb 2020 19:26:22 -0800 Subject: [PATCH 12/16] Add setuptools_scm --- MANIFEST.in | 7 +++++++ pyproject.toml | 5 +++++ requirements.txt | 5 ++++- requirements_travis.txt | 5 ++++- setup.py | 10 +++++----- 5 files changed, 25 insertions(+), 7 deletions(-) create mode 100644 MANIFEST.in create mode 100644 pyproject.toml diff --git a/MANIFEST.in b/MANIFEST.in new file mode 100644 index 00000000..5ce142c1 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 00000000..9ccf0a69 --- /dev/null +++ b/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" diff --git a/requirements.txt b/requirements.txt index 4eba8deb..4ff38dcd 100644 --- a/requirements.txt +++ b/requirements.txt @@ -15,4 +15,7 @@ numba spiceypy pytest pytest-cov -coveralls \ No newline at end of file +coveralls +setuptools>=42 +wheel +setuptools_scm>=3.4 diff --git a/requirements_travis.txt b/requirements_travis.txt index d96e60c4..0874bede 100644 --- a/requirements_travis.txt +++ b/requirements_travis.txt @@ -14,4 +14,7 @@ numba spiceypy pytest pytest-cov<2.6.0 -coveralls \ No newline at end of file +coveralls +setuptools>=42 +wheel +setuptools_scm>=3.4 diff --git a/setup.py b/setup.py index 6bd719ad..7b0fbe1a 100644 --- a/setup.py +++ b/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"], ) From d24b598f0a927609c95176b39298a81e861b7845 Mon Sep 17 00:00:00 2001 From: Joachim Moeyens Date: Thu, 13 Feb 2020 08:40:21 -0800 Subject: [PATCH 13/16] Update gauss unit tests --- thor/orbits/iod/tests/test_gauss.py | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/thor/orbits/iod/tests/test_gauss.py b/thor/orbits/iod/tests/test_gauss.py index ea016667..e6b5e5f4 100644 --- a/thor/orbits/iod/tests/test_gauss.py +++ b/thor/orbits/iod/tests/test_gauss.py @@ -87,19 +87,20 @@ def test_gaussIOD(): for i, observation in enumerate(coords_eq_ang): print("\tt [MJD]: {}, RA [Deg]: {}, Dec [Deg]: {}, obs_x [AU]: {}, obs_y [AU]: {}, obs_z [AU]: {}".format(t[i], *observation, *coords_obs[i])) - print("Actual State [AU, AU/d]:\n\t{}".format(states_target[selected_obs[1]])) + state_truth = np.concatenate([np.array([t[1]]), states_target[selected_obs[1]]]) + print("Actual State [MJD, AU, AU/d]:\n\t{}".format(state_truth)) # Using observables, run IOD without iteration - orbits = gaussIOD(coords_eq_ang, t, coords_obs, velocity_method="gibbs", iterate=False, mu=MU, max_iter=100, tol=1e-15) + orbits = gaussIOD(coords_eq_ang, t, coords_obs, velocity_method="gibbs", iterate=False, mu=MU, max_iter=100, tol=1e-15, light_time=False) - print("Predicted States (without iteration) [AU, AU/d]:") + print("Predicted States (without iteration) [MJD, AU, AU/d]:") for orbit in orbits: print("\t{}".format(orbit)) # Using observables, run IOD - orbits = gaussIOD(coords_eq_ang, t, coords_obs, velocity_method="gibbs", iterate=True, mu=MU, max_iter=100, tol=1e-15) + orbits = gaussIOD(coords_eq_ang, t, coords_obs, velocity_method="gibbs", iterate=True, mu=MU, max_iter=100, tol=1e-15, light_time=False) - print("Predicted States (with iteration) [AU, AU/d]:") + print("Predicted States (with iteration) [MJD, AU, AU/d]:") for orbit in orbits: print("\t{}".format(orbit)) @@ -110,8 +111,8 @@ def test_gaussIOD(): closest_v = 1e10 for i, orbit in enumerate(orbits): - r2 = orbit[:3] - v2 = orbit[3:] + r2 = orbit[1:4] + v2 = orbit[4:] r2_mag = np.linalg.norm(r2) v2_mag = np.linalg.norm(v2) @@ -129,7 +130,7 @@ def test_gaussIOD(): print("(Actual - Predicted) / Actual:") for orbit in orbits: - print("\t{}".format((states_target[selected_obs[1]] - orbit) / states_target[selected_obs[1]])) + print("\t{}".format((states_target[selected_obs[1]] - orbit[1:]) / states_target[selected_obs[1]])) print("") # Test position to within 100 meters and velocity to within 10 cm/s From 06a43faecf617598f4dccbf4be2185a105dac386 Mon Sep 17 00:00:00 2001 From: Joachim Moeyens Date: Thu, 13 Feb 2020 09:42:14 -0800 Subject: [PATCH 14/16] Fix GM to be the same as DE430/DE431 --- thor/constants.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/thor/constants.py b/thor/constants.py index d848782e..d7568fef 100644 --- a/thor/constants.py +++ b/thor/constants.py @@ -9,8 +9,8 @@ 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 @@ -18,7 +18,7 @@ class Constants: # 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) From 5d2c4f7bda7de32abe241349b57e692f0824001c Mon Sep 17 00:00:00 2001 From: Joachim Moeyens Date: Thu, 13 Feb 2020 09:45:07 -0800 Subject: [PATCH 15/16] Clean up tests_kepler, make use of to_pandas() --- thor/orbits/tests/test_kepler.py | 195 ++++++++++++++++++------------- 1 file changed, 114 insertions(+), 81 deletions(-) diff --git a/thor/orbits/tests/test_kepler.py b/thor/orbits/tests/test_kepler.py index e98eec95..907cc00a 100644 --- a/thor/orbits/tests/test_kepler.py +++ b/thor/orbits/tests/test_kepler.py @@ -10,114 +10,147 @@ from ..kepler import _convertKeplerianToCartesian MU = c.G * c.M_SUN +TOL = 1e-15 +MAX_ITER = 100 + +EPOCHS = Time( + ["{}-02-02T00:00:00.000".format(i) for i in range(1993, 2050)], + format="isot", + scale="tdb" +) +ISO_EPOCHS = Time( + ["{}-02-02T00:00:00.000".format(i) for i in range(2017, 2022)], + format="isot", + scale="tdb" +) +TARGETS = [ + "Amor", + "Eros", + "Eugenia", + "Ceres", +] +ISO_TARGETS = [ + "1I/2017 U1", # Oumuamua + "C/2019 Q4" # Borisov +] + def test_convertOrbitalElements(): - epochs = Time(["{}-02-02T00:00:00.000".format(i) for i in range(1993, 2050)], format="isot", scale="tdb") - targets = [ - "Amor", - "Eros", - "Eugenia", - "Ceres" - ] - for name in targets: - target = Horizons(id=name, epochs=epochs.mjd) - vectors = target.vectors() - vectors = np.array(vectors["x", "y", "z", "vx", "vy", "vz"]) - vectors = vectors.view("float64").reshape(vectors.shape + (-1,)) + for name in TARGETS: + target = Horizons(id=name, epochs=EPOCHS.mjd) + + vectors = target.vectors().to_pandas() + vectors = vectors[["x", "y", "z", "vx", "vy", "vz"]].values + + elements = target.elements().to_pandas() + elements = elements[["a", "e", "incl", "Omega", "w", "M", "nu"]].values - elements = target.elements() - elements = np.array(elements["a", "e", "incl", "Omega", "w", "M", "nu"]) - elements = elements.view("float64").reshape(elements.shape + (-1,)) + np.testing.assert_allclose( + elements[:, :6], + convertOrbitalElements( + vectors, + "cartesian", + "keplerian", + max_iter=MAX_ITER, + tol=TOL, + ) + ) - np.testing.assert_allclose(elements[:, :6], convertOrbitalElements(vectors, "cartesian", "keplerian")) - np.testing.assert_allclose(vectors, convertOrbitalElements(elements[:, :6], "keplerian", "cartesian")) + np.testing.assert_allclose( + vectors, + convertOrbitalElements( + elements[:, :6], + "keplerian", + "cartesian", + max_iter=MAX_ITER, + tol=TOL, + ) + ) + return def test_convertCartesianToKeplerian_elliptical(): - epochs = Time(["{}-02-02T00:00:00.000".format(i) for i in range(1993, 2050)], format="isot", scale="tdb") - targets = [ - "Amor", - "Eros", - "Eugenia", - "Ceres", - ] - - for name in targets: - target = Horizons(id=name, epochs=epochs.mjd) - vectors = target.vectors() - vectors = np.array(vectors["x", "y", "z", "vx", "vy", "vz"]) - vectors = vectors.view("float64").reshape(vectors.shape + (-1,)) + for name in TARGETS: + target = Horizons(id=name, epochs=EPOCHS.mjd) + vectors = target.vectors().to_pandas() + vectors = vectors[["x", "y", "z", "vx", "vy", "vz"]].values - elements = target.elements() - elements = np.array(elements["a", "q", "e", "incl", "Omega", "w", "M", "nu"]) - elements = elements.view("float64").reshape(elements.shape + (-1,)) + elements = target.elements().to_pandas() + elements = elements[["a", "q", "e", "incl", "Omega", "w", "M", "nu"]].values for v, e in zip(vectors, elements): - np.testing.assert_allclose(_convertCartesianToKeplerian(v.reshape(1, -1), mu=MU), e.reshape(1, -1)) + np.testing.assert_allclose( + _convertCartesianToKeplerian( + v.reshape(1, -1), + mu=MU, + ), + e.reshape(1, -1) + ) + return def test_convertCartesianToKeplerian_parabolic(): warnings.warn("Need to implement and test parabolic conversions!!!") + return def test_convertCartesianToKeplerian_hyperbolic(): - epochs = Time(["{}-02-02T00:00:00.000".format(i) for i in range(2017, 2023)], format="isot", scale="tdb") - iso_targets = [ - "1I/2017 U1", #Oumuamua - "C/2019 Q4" #Borisov - ] - - for name in iso_targets: - target = Horizons(id=name, epochs=epochs.mjd) - vectors = target.vectors() - vectors = np.array(vectors["x", "y", "z", "vx", "vy", "vz"]) - vectors = vectors.view("float64").reshape(vectors.shape + (-1,)) + for name in ISO_TARGETS: + target = Horizons(id=name, epochs=ISO_EPOCHS.mjd) + vectors = target.vectors().to_pandas() + vectors = vectors[["x", "y", "z", "vx", "vy", "vz"]].values - elements = target.elements() - elements = np.array(elements["a", "q", "e", "incl", "Omega", "w", "M", "nu"]) - elements = elements.view("float64").reshape(elements.shape + (-1,)) + elements = target.elements().to_pandas() + elements = elements[["a", "q", "e", "incl", "Omega", "w", "M", "nu"]].values for v, e in zip(vectors, elements): - np.testing.assert_allclose(_convertCartesianToKeplerian(v.reshape(1, -1), mu=MU), e.reshape(1, -1)) + np.testing.assert_allclose( + _convertCartesianToKeplerian( + v.reshape(1, -1), + mu=MU, + ), + e.reshape(1, -1) + ) + return def test_convertKeplerianToCartesian_elliptical(): - epochs = Time(["{}-02-02T00:00:00.000".format(i) for i in range(1993, 2050)], format="isot", scale="tdb") - targets = [ - "Amor", - "Eros", - "Eugenia", - "Ceres" - ] - - for name in targets: - target = Horizons(id=name, epochs=epochs.mjd) - vectors = target.vectors() - vectors = np.array(vectors["x", "y", "z", "vx", "vy", "vz"]) - vectors = vectors.view("float64").reshape(vectors.shape + (-1,)) + for name in TARGETS: + target = Horizons(id=name, epochs=EPOCHS.mjd) + vectors = target.vectors().to_pandas() + vectors = vectors[["x", "y", "z", "vx", "vy", "vz"]].values - elements = target.elements() - elements = np.array(elements["a", "e", "incl", "Omega", "w", "M", "nu"]) - elements = elements.view("float64").reshape(elements.shape + (-1,)) + elements = target.elements().to_pandas() + elements = elements[["a", "e", "incl", "Omega", "w", "M", "nu"]].values for v, e in zip(vectors, elements): - np.testing.assert_allclose(_convertKeplerianToCartesian(e.reshape(1, -1), mu=MU, max_iter=100, tol=1e-15), v.reshape(1, -1)) + np.testing.assert_allclose( + _convertKeplerianToCartesian( + e.reshape(1, -1), + mu=MU, + max_iter=MAX_ITER, + tol=TOL, + ), + v.reshape(1, -1) + ) + return def test_convertKeplerianToCartesian_parabolic(): warnings.warn("Need to implement and test parabolic conversions!!!") def test_convertKeplerianToCartesian_hyperbolic(): - epochs = Time(["{}-02-02T00:00:00.000".format(i) for i in range(2017, 2023)], format="isot", scale="tdb") - iso_targets = [ - "1I/2017 U1", #Oumuamua - "C/2019 Q4" #Borisov - ] - - for name in iso_targets: - target = Horizons(id=name, epochs=epochs.mjd) - vectors = target.vectors() - vectors = np.array(vectors["x", "y", "z", "vx", "vy", "vz"]) - vectors = vectors.view("float64").reshape(vectors.shape + (-1,)) + for name in ISO_TARGETS: + target = Horizons(id=name, epochs=ISO_EPOCHS.mjd) + vectors = target.vectors().to_pandas() + vectors = vectors[["x", "y", "z", "vx", "vy", "vz"]].values - elements = target.elements() - elements = np.array(elements["a", "e", "incl", "Omega", "w", "M"]) - elements = elements.view("float64").reshape(elements.shape + (-1,)) + elements = target.elements().to_pandas() + elements = elements[["a", "e", "incl", "Omega", "w", "M", "nu"]].values for v, e in zip(vectors, elements): - np.testing.assert_allclose(_convertKeplerianToCartesian(e.reshape(1, -1), mu=MU, max_iter=100, tol=1e-15), v.reshape(1, -1)) \ No newline at end of file + np.testing.assert_allclose( + _convertKeplerianToCartesian( + e.reshape(1, -1), + mu=MU, + max_iter=MAX_ITER, + tol=TOL, + ), + v.reshape(1, -1) + ) + return \ No newline at end of file From 749ffedb33c445e644e5002e9bd89f4bce1ae077 Mon Sep 17 00:00:00 2001 From: Joachim Moeyens Date: Thu, 13 Feb 2020 09:56:01 -0800 Subject: [PATCH 16/16] Add additional test epoch in test_propagate.py, make use of to_pandas() --- thor/orbits/propagate/tests/test_propagate.py | 38 +++++++++++-------- 1 file changed, 23 insertions(+), 15 deletions(-) diff --git a/thor/orbits/propagate/tests/test_propagate.py b/thor/orbits/propagate/tests/test_propagate.py index d0f62ddf..b4ba7555 100644 --- a/thor/orbits/propagate/tests/test_propagate.py +++ b/thor/orbits/propagate/tests/test_propagate.py @@ -6,29 +6,30 @@ from ..universal import propagateUniversal MU = c.G * c.M_SUN +MAX_ITER = 100 +TOL = 1e-15 + +TARGETS = [ + "Amor", + "Eros", + "Eugenia", + "C/2019 Q4" # Borisov +] +EPOCHS = [57257.0, 59000.0] def test_propagateUniversal(): """ Using a selection of 4 asteroids, this function queries Horizons for an initial state vector at one epoch, then propagates that state to 1000 different times and compares each propagation to the SPICE 2-body propagator. """ - targets = [ - "Amor", - "Eros", - "Eugenia", - "C/2019 Q4" #Borisov - ] - - epochs = [57257.0] dts = np.linspace(0.01, 500, num=1000) - for name in targets: - for epoch in epochs: + for name in TARGETS: + for epoch in EPOCHS: # Grab vectors from Horizons at epoch target = Horizons(id=name, epochs=epoch, location="@sun") - vectors = target.vectors() - vectors = np.array(vectors["x", "y", "z", "vx", "vy", "vz"]).view("float64") - vectors = vectors.reshape(-1, 6) + vectors = target.vectors().to_pandas() + vectors = vectors[["x", "y", "z", "vx", "vy", "vz"]].values # Propagate vector to each new epoch (epoch + dt) spice_elements = [] @@ -37,7 +38,14 @@ def test_propagateUniversal(): spice_elements = np.array(spice_elements) # Repeat but now using THOR's universal propagator - vectors_new = propagateUniversal(vectors[0:1, :], np.array(epochs), dts + epochs[0], mu=MU, max_iter=1000, tol=1e-15) + vectors_new = propagateUniversal( + vectors[0:1, :], + np.array([epoch]), + dts + epoch, + mu=MU, + max_iter=MAX_ITER, + tol=TOL + ) orbit_id = vectors_new[:, 0] new_epochs = vectors_new[:, 1] @@ -45,7 +53,7 @@ def test_propagateUniversal(): # Make sure the first column is a bunch of 0s since only one orbit was passed np.testing.assert_allclose(orbit_id, np.zeros(len(dts))) # Make sure the second column has all the new epochs - np.testing.assert_allclose(new_epochs, dts + epochs[0]) + np.testing.assert_allclose(new_epochs, dts + epoch) # Extract position and velocity components and compare them r = vectors_new[:, 2:5]