Skip to content

Commit 56fb657

Browse files
authored
Merge pull request #33 from moeyensj/gauss
Gauss
2 parents edaa3d1 + 749ffed commit 56fb657

13 files changed

+350
-140
lines changed

MANIFEST.in

+7
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,7 @@
1+
include thor/*.py
2+
include *.md
3+
include *.toml
4+
exclude .travis.yaml
5+
exclude requirements_travis.txt
6+
exclude azure-pipelines.yml
7+
exclude Dockerfile

README.md

+3-2
Original file line numberDiff line numberDiff line change
@@ -4,8 +4,9 @@ Tracklet-less Heliocentric Orbit Recovery
44
[![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)
55
[![Build Status](https://www.travis-ci.com/moeyensj/thor.svg?token=sWjpnqPgpHyuq3j7qPuj&branch=master)](https://www.travis-ci.com/moeyensj/thor)
66
[![Coverage Status](https://coveralls.io/repos/github/moeyensj/thor/badge.svg?branch=master&t=pdSkQA)](https://coveralls.io/github/moeyensj/thor?branch=master)
7-
[![Docker Pulls](https://img.shields.io/docker/pulls/moeyensj/thor)](https://hub.docker.com/r/moeyensj/thor)
8-
[![License](https://img.shields.io/badge/License-BSD%203--Clause-blue.svg)](https://opensource.org/licenses/BSD-3-Clause)
7+
[![Docker Pulls](https://img.shields.io/docker/pulls/moeyensj/thor)](https://hub.docker.com/r/moeyensj/thor)
8+
[![Python 3.6](https://img.shields.io/badge/Python-3.6%2B-blue)](https://img.shields.io/badge/Python-3.6%2B-blue)
9+
[![License](https://img.shields.io/badge/License-BSD%203--Clause-blue.svg)](https://opensource.org/licenses/BSD-3-Clause)
910

1011
## Installation
1112

pyproject.toml

+5
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,5 @@
1+
[build-system]
2+
requires = ["setuptools>=42", "wheel", "setuptools_scm[toml]>=3.4"]
3+
4+
[tool.setuptools_scm]
5+
write_to = "pkg/version.py"

requirements.txt

+4-1
Original file line numberDiff line numberDiff line change
@@ -15,4 +15,7 @@ numba
1515
spiceypy
1616
pytest
1717
pytest-cov
18-
coveralls
18+
coveralls
19+
setuptools>=42
20+
wheel
21+
setuptools_scm>=3.4

requirements_travis.txt

+4-1
Original file line numberDiff line numberDiff line change
@@ -14,4 +14,7 @@ numba
1414
spiceypy
1515
pytest
1616
pytest-cov<2.6.0
17-
coveralls
17+
coveralls
18+
setuptools>=42
19+
wheel
20+
setuptools_scm>=3.4

setup.py

+5-5
Original file line numberDiff line numberDiff line change
@@ -2,16 +2,16 @@
22

33
setup(
44
name="thor",
5-
version="0.0.1.dev0",
6-
description="Tracklet-less Heliocentric Orbit Recovery",
75
license="BSD 3-Clause License",
86
author="Joachim Moeyens, Mario Juric",
97
author_email="[email protected]",
8+
long_description=open("README.md").read(),
9+
long_description_content_type="text/markdown",
1010
url="https://github.com/moeyensj/thor",
1111
packages=["thor"],
1212
package_dir={"thor": "thor"},
13-
package_data={"thor": ["data/*.orb",
14-
"tests/data/*"]},
15-
setup_requires=["pytest-runner"],
13+
package_data={"thor": ["tests/*.txt"]},
14+
use_scm_version=True,
15+
setup_requires=["pytest-runner", "setuptools_scm"],
1616
tests_require=["pytest"],
1717
)

thor/constants.py

+3-3
Original file line numberDiff line numberDiff line change
@@ -9,16 +9,16 @@ class Constants:
99
# Speed of light: AU per day (173.14463267424034)
1010
C = c.c.to(u.au / u.d).value
1111

12-
# Gravitational constant: AU**3 / M_sun / d**2 (0.00029591220819207784)
13-
G = c.G.to(u.AU**3 / u.M_sun / u.d**2).value
12+
# Gravitational constant: AU**3 / M_sun / d**2 (0.295912208285591100E-3 -- DE431/DE430)
13+
G = 0.295912208285591100E-3
1414

1515
# Solar Mass: M_sun (1.0)
1616
M_SUN = 1.0
1717

1818
# Earth Mass: M_sun (3.0034893488507934e-06)
1919
M_EARTH = u.M_earth.to(u.M_sun)
2020

21-
# Earth Equatorial Radius (6378.1363 km (DE431/DE430))
21+
# Earth Equatorial Radius: km (6378.1363 -- DE431/DE430)
2222
R_EARTH = (6378.1363 * u.km).to(u.AU)
2323

2424
# Mean Obliquity at J2000: radians (0.40909280422232897)

thor/orbits/iod/gauss.py

+34-8
Original file line numberDiff line numberDiff line change
@@ -71,7 +71,10 @@ def _calcFG(r2_mag, t32, t21, mu=MU):
7171

7272
def _calcM(r0_mag, r_mag, f, g, f_dot, g_dot, c0, c1, c2, c3, c4, c5, alpha, chi, mu=MU):
7373
# Universal variables will differ between different texts and works in the literature.
74-
# c0, c1, c2, c3, c4, c5 are expected to be
74+
# c0, c1, c2, c3, c4, c5 are expected to follow the Battin formalism (adopted by both
75+
# Vallado and Curtis in their books). The M matrix is proposed by Shepperd 1985 and follows
76+
# the Goodyear formalism. Conversions between the two formalisms can be derived from Table 1 in
77+
# Everhart & Pitkin 1982.
7578
w = chi / np.sqrt(mu)
7679
alpha_alt = - mu * alpha
7780
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
8487
F = f_dot
8588
G = g_dot
8689

87-
# Equations 18 and 19 in Sheppard 1985
90+
# Equations 18 and 19 in Shepperd 1985
8891
U = (U2 * U3 + w * U4 - 3 * U5) / 3
8992
W = g * U2 + 3 * mu * U
9093

@@ -124,7 +127,7 @@ def _calcStateTransitionMatrix(M, r0, v0, f, g, f_dot, g_dot, r, v):
124127
])
125128
return phi
126129

127-
def _iterateGaussIOD(orbit, t21, t32, q1, q2, q3, rho1, rho2, rho3, mu=MU, max_iter=10, tol=1e-15):
130+
def _iterateGaussIOD(orbit, t21, t32, q1, q2, q3, rho1, rho2, rho3, light_time=True, mu=MU, max_iter=10, tol=1e-15):
128131
# Iterate over the polynomial solution from Gauss using the universal anomaly
129132
# formalism until the solution converges or the maximum number of iterations is reached
130133

@@ -159,6 +162,12 @@ def _iterateGaussIOD(orbit, t21, t32, q1, q2, q3, rho1, rho2, rho3, mu=MU, max_i
159162
# then calculate the Lagrange coefficients and the state for each observation.
160163
# Use those to calculate the state transition matrix
161164
for j, dt in enumerate([-t21, t32]):
165+
if light_time is True:
166+
if j == 1:
167+
dt += (rho2_mag - rho1_mag) / C
168+
else:
169+
dt -= (rho3_mag - rho2_mag) / C
170+
162171
# Calculate the universal anomaly
163172
# Universal anomaly here is defined in such a way that it satisfies the following
164173
# differential equation:
@@ -286,7 +295,7 @@ def calcGauss(r1, r2, r3, t1, t2, t3):
286295
f1, g1, f3, g3 = _calcFG(r2_mag, t32, t21)
287296
return (1 / (f1 * g3 - f3 * g1)) * (-f3 * r1 + f1 * r3)
288297

289-
def gaussIOD(coords_eq_ang, t, coords_obs, velocity_method="gibbs", iterate=True, mu=MU, max_iter=10, tol=1e-15):
298+
def gaussIOD(coords_eq_ang, t, coords_obs, velocity_method="gibbs", light_time=True, iterate=True, mu=MU, max_iter=10, tol=1e-15):
290299
"""
291300
Compute up to three intial orbits using three observations in angular equatorial
292301
coordinates.
@@ -302,6 +311,9 @@ def gaussIOD(coords_eq_ang, t, coords_obs, velocity_method="gibbs", iterate=True
302311
velocity_method : {'gauss', gibbs', 'herrick+gibbs'}, optional
303312
Which method to use for calculating the velocity at the second observation.
304313
[Default = 'gibbs']
314+
light_time : bool, optional
315+
Correct for light travel time.
316+
[Default = True]
305317
iterate : bool, optional
306318
Iterate initial orbit using universal anomaly to better approximate the
307319
Lagrange coefficients.
@@ -346,7 +358,7 @@ def gaussIOD(coords_eq_ang, t, coords_obs, velocity_method="gibbs", iterate=True
346358
coseps2 = np.dot(q2, rho2_hat) / q2_mag
347359
C0 = V * t31 * q2_mag**4 / B
348360
h0 = - A / B
349-
361+
350362
if np.isnan(C0) or np.isnan(h0):
351363
return np.array([])
352364

@@ -376,7 +388,7 @@ def gaussIOD(coords_eq_ang, t, coords_obs, velocity_method="gibbs", iterate=True
376388

377389
# Test if we get the same rho2 as using equation 22 in Milani et al. 2008
378390
rho2_mag = (h0 - q2_mag**3 / r2_mag**3) * q2_mag / C0
379-
np.testing.assert_almost_equal(np.dot(rho2_mag, rho2_hat), rho2)
391+
#np.testing.assert_almost_equal(np.dot(rho2_mag, rho2_hat), rho2)
380392

381393
r1 = q1 + rho1
382394
r2 = q2 + rho2
@@ -390,16 +402,30 @@ def gaussIOD(coords_eq_ang, t, coords_obs, velocity_method="gibbs", iterate=True
390402
v2 = calcHerrickGibbs(r1, r2, r3, t1, t2, t3)
391403
else:
392404
raise ValueError("velocity_method should be one of {'gauss', 'gibbs', 'herrick+gibbs'}")
405+
406+
epoch = t2
393407
orbit = np.concatenate([r2, v2])
394408

395409
if iterate == True:
396410
orbit = _iterateGaussIOD(orbit, t21, t32,
397411
q1, q2, q3,
398412
rho1, rho2, rho3,
399-
mu=mu, max_iter=max_iter, tol=tol)
413+
light_time=light_time,
414+
mu=mu,
415+
max_iter=max_iter,
416+
tol=tol)
417+
418+
if light_time == True:
419+
rho2_mag = np.linalg.norm(orbit[:3] - q2)
420+
lt = rho2_mag / C
421+
epoch -= lt
400422

401423
if np.linalg.norm(orbit[3:]) >= C:
402424
print("Velocity is greater than speed of light!")
403-
orbits.append(orbit)
425+
426+
if (np.linalg.norm(orbit[:3]) > 300.) and (np.linalg.norm(orbit[3:]) > 25.):
427+
continue
428+
429+
orbits.append(np.hstack([epoch, orbit]))
404430

405431
return np.array(orbits)

thor/orbits/iod/iod.py

+133-14
Original file line numberDiff line numberDiff line change
@@ -1,46 +1,165 @@
1+
import sys
12
import numpy as np
3+
from itertools import combinations
24

35
from ...config import Config
6+
from ...constants import Constants as c
7+
from .gauss import gaussIOD
8+
from ..propagate import propagateOrbits
49

5-
__all__ = ["selectObservations"]
10+
__all__ = ["selectObservations",
11+
"iod"]
612

7-
def selectObservations(observations, method="first+middle+last", columnMapping=Config.columnMapping):
13+
MU = c.G * c.M_SUN
14+
15+
16+
def selectObservations(observations, method="combinations", columnMapping=Config.columnMapping):
817
"""
918
Selects which three observations to use for IOD depending on the method.
1019
1120
Methods:
1221
'first+middle+last' : Grab the first, middle and last observations in time.
1322
'thirds' : Grab the middle observation in the first third, second third, and final third.
23+
'combinations' : Return the observation IDs corresponding to every possible combination of three observations with
24+
non-coinciding observation times.
1425
1526
Parameters
1627
----------
1728
observations : `~pandas.DataFrame`
1829
Pandas DataFrame containing observations with at least a column of observation IDs and a column
1930
of exposure times.
20-
method : {'first+middle+last', 'thirds}, optional
31+
method : {'first+middle+last', 'thirds', 'combinations'}, optional
2132
Which method to use to select observations.
22-
[Default = `first+middle+last`]
33+
[Default = 'combinations']
2334
columnMapping : dict, optional
2435
Column name mapping of observations to internally used column names.
2536
[Default = `~thor.Config.columnMapping`]
2637
2738
Returns
2839
-------
29-
obs_id : `~numpy.ndarray' (3 or 0)
40+
obs_id : `~numpy.ndarray' (N, 3 or 0)
3041
An array of selected observation IDs. If three unique observations could
31-
not be selected than returns an empty array.
42+
not be selected then returns an empty array.
3243
"""
44+
obs_ids = observations[columnMapping["obs_id"]].values
45+
if len(obs_ids) < 3:
46+
return np.array([])
47+
48+
indexes = np.arange(0, len(obs_ids))
49+
times = observations[columnMapping["exp_mjd"]].values
50+
selected = np.array([])
51+
3352
if method == "first+middle+last":
34-
selected = np.percentile(observations[columnMapping["exp_mjd"]].values, [0, 50, 100], interpolation="nearest")
53+
selected_times = np.percentile(times,
54+
[0, 50, 100],
55+
interpolation="nearest")
56+
selected_index = np.intersect1d(times, selected_times, return_indices=True)[1]
57+
selected_index = np.array([selected_index])
58+
3559
elif method == "thirds":
36-
selected = np.percentile(observations[columnMapping["exp_mjd"]].values, [1/6 * 100, 50, 5/6*100], interpolation="nearest")
60+
selected_times = np.percentile(times,
61+
[1/6*100, 50, 5/6*100],
62+
interpolation="nearest")
63+
selected_index = np.intersect1d(times, selected_times, return_indices=True)[1]
64+
selected_index = np.array([selected_index])
65+
66+
elif method == "combinations":
67+
# Make all possible combinations of 3 observations
68+
selected_index = np.array([np.array(index) for index in combinations(indexes, 3)])
69+
3770
else:
3871
raise ValueError("method should be one of {'first+middle+last', 'thirds'}")
39-
40-
if len(np.unique(selected)) != 3:
41-
print("Could not find three observations that satisfy the criteria.")
72+
73+
# Make sure each returned combination of observation ids have at least 3 unique
74+
# times
75+
keep = []
76+
for i, comb in enumerate(times[selected_index]):
77+
if len(np.unique(comb)) == 3:
78+
keep.append(i)
79+
keep = np.array(keep)
80+
81+
# Return an empty array if no observations satisfy the criteria
82+
if len(keep) == 0:
4283
return np.array([])
4384

44-
index = np.intersect1d(observations[columnMapping["exp_mjd"]].values, selected, return_indices=True)[1]
45-
return observations[columnMapping["obs_id"]].values[index]
46-
85+
return obs_ids[selected_index[keep, :]]
86+
87+
88+
def iod(observations,
89+
observation_selection_method="combinations",
90+
iterate=True,
91+
light_time=True,
92+
max_iter=50,
93+
tol=1e-15,
94+
propagatorKwargs={
95+
"observatoryCode" : "I11",
96+
"mjdScale" : "UTC",
97+
"dynamical_model" : "2",
98+
},
99+
mu=MU,
100+
columnMapping=Config.columnMapping):
101+
102+
# Extract column names
103+
obs_id_col = columnMapping["obs_id"]
104+
ra_col = columnMapping["RA_deg"]
105+
dec_col = columnMapping["Dec_deg"]
106+
time_col = columnMapping["exp_mjd"]
107+
ra_err_col = columnMapping["RA_sigma_deg"]
108+
dec_err_col = columnMapping["Dec_sigma_deg"]
109+
obs_x_col = columnMapping["obs_x_au"]
110+
obs_y_col = columnMapping["obs_y_au"]
111+
obs_z_col = columnMapping["obs_z_au"]
112+
113+
# Extract observation IDs, sky-plane positions, sky-plane position uncertainties, times of observation,
114+
# and the location of the observer at each time
115+
obs_ids_all = observations[obs_id_col].values
116+
coords_eq_ang_all = observations[observations[obs_id_col].isin(obs_ids_all)][[ra_col, dec_col]].values
117+
coords_eq_ang_err_all = observations[observations[obs_id_col].isin(obs_ids_all)][[ra_err_col, dec_err_col]].values
118+
coords_obs_all = observations[observations[obs_id_col].isin(obs_ids_all)][[obs_x_col, obs_y_col, obs_z_col]].values
119+
times_all = observations[observations[obs_id_col].isin(obs_ids_all)][time_col].values
120+
121+
# Select observation IDs to use for IOD
122+
obs_ids = selectObservations(observations, method=observation_selection_method, columnMapping=columnMapping)
123+
124+
min_chi2 = 1e10
125+
best_orbit = None
126+
best_obs_ids = None
127+
128+
for ids in obs_ids:
129+
# Grab sky-plane positions of the selected observations, the heliocentric ecliptic position of the observer,
130+
# and the times at which the observations occur
131+
coords_eq_ang = observations[observations[obs_id_col].isin(ids)][[ra_col, dec_col]].values
132+
coords_obs = observations[observations[obs_id_col].isin(ids)][[obs_x_col, obs_y_col, obs_z_col]].values
133+
times = observations[observations[obs_id_col].isin(ids)][time_col].values
134+
135+
# Run IOD
136+
orbits_iod = gaussIOD(coords_eq_ang, times, coords_obs, light_time=light_time, iterate=iterate, max_iter=max_iter, tol=tol)
137+
if np.all(np.isnan(orbits_iod)) == True:
138+
continue
139+
140+
# Propagate initial orbit to all observation times
141+
orbits = propagateOrbits(orbits_iod[:, 1:], orbits_iod[:, 0], times_all, **propagatorKwargs)
142+
orbits = orbits[['orbit_id', 'mjd', 'RA_deg', 'Dec_deg',
143+
'HEclObj_X_au', 'HEclObj_Y_au', 'HEclObj_Z_au',
144+
'HEclObj_dX/dt_au_p_day', 'HEclObj_dY/dt_au_p_day', 'HEclObj_dZ/dt_au_p_day']].values
145+
146+
# For each unique initial orbit calculate residuals and chi-squared
147+
# Find the orbit which yields the lowest chi-squared
148+
orbit_ids = np.unique(orbits[:, 0])
149+
for i, orbit_id in enumerate(orbit_ids):
150+
orbit = orbits[np.where(orbits[:, 0] == orbit_id)]
151+
152+
pred_dec = np.radians(orbit[:, 3])
153+
residual_ra = (coords_eq_ang_all[:, 0] - orbit[:, 2]) * np.cos(pred_dec)
154+
residual_dec = (coords_eq_ang_all[:, 1] - orbit[:, 3])
155+
156+
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)
157+
158+
if chi2 < min_chi2:
159+
best_orbit = orbits_iod[i, :]
160+
best_obs_ids = ids
161+
min_chi2 = chi2
162+
163+
return best_orbit, best_obs_ids, min_chi2
164+
165+

0 commit comments

Comments
 (0)