From ba639cce96d463c6b9b4d62b37e6601f39f7d33c Mon Sep 17 00:00:00 2001 From: Dar Dahlen Date: Fri, 27 Sep 2024 15:19:36 -0700 Subject: [PATCH] Add 2010 TK7 trojan analysis example (#126) * Add 2010 TK7 example * fix citations --- CHANGELOG.md | 4 ++ src/examples/earth_trojan.py | 96 ++++++++++++++++++++++++++++++++++++ src/kete/horizons.py | 12 ++--- 3 files changed, 103 insertions(+), 9 deletions(-) create mode 100644 src/examples/earth_trojan.py diff --git a/CHANGELOG.md b/CHANGELOG.md index d0ebeb1..b39eab1 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +### Added + +- Added worked example of calculating if an object is an Earth Trojan. + ### Fixed - Epoch and Perihelion time conversion when loading JPL Horizons covariance matrices was diff --git a/src/examples/earth_trojan.py b/src/examples/earth_trojan.py new file mode 100644 index 0000000..c258a9d --- /dev/null +++ b/src/examples/earth_trojan.py @@ -0,0 +1,96 @@ +""" +Earth Trojan 2010 TK7 +===================== + +Following some of the analysis from the following papers: + +- `"2002 AA29: Earth's recurrent quasi-satellite?" - Paweł Wajer - Icarus - 2009 `_ +- `"Dynamical evolution of Earth’s quasi-satellites: 2004 GU9 and 2006 FV35" - Paweł Wajer - Icarus - 2010 `_ + +We can perform some dynamical analysis to see if the asteroid 2010 TK7 appears +to be an Earth trojan. + +The end result of this is a plot which shows the current energy level of the +orbit, if the current position of the object is in a gravitational "low", and +the energy line (in red) is below the height of the walls of that low point, then +the object is at least temporarily captured. + +""" + +import numpy as np +import matplotlib.pyplot as plt +import kete + + +# Grab 2010 TK7 and propagate it to 2010 +state_a = kete.HorizonsProperties.fetch("2010 TK7", update_name=False).state +state_a = kete.propagate_n_body([state_a], kete.Time.from_ymd(2010, 1, 1))[0] + +state_b = kete.spice.get_state("Earth", state_a.jd) +# Earth mass in Solar masses +mass_planet = 3.0435727639549377e-06 + +# n_sigma is the number of step on the x axis in the final plot +n_sigma = 500 + +# m_steps is the number of steps in the numerical integration +m_steps = 100 + + +jd = state_a.jd + +elem_a = state_a.elements +elem_b = state_b.elements + +period_a = elem_a.orbital_period +period_b = elem_b.orbital_period + +sigma_steps = np.linspace(0, 1, n_sigma + 1)[:-1] + +energy_vals = [] +for shift in sigma_steps: + + total = 0 + for frac in np.linspace(0, 1, m_steps + 1)[:-1]: + state_a = kete.propagate_two_body([state_a], jd + (frac + shift) * period_a)[0] + state_b = kete.propagate_two_body([state_b], jd + frac * period_b)[0] + pos_a = state_a.pos + pos_b = state_b.pos + total += 1 / (pos_b - pos_a).r - np.dot(pos_b, pos_a) / pos_b.r**3 + + energy_vals.append(total / m_steps) + +cur_energy = 3 * (elem_a.semi_major - 1) ** 2 / (8.0 * mass_planet) + energy_vals[0] + +# energy vals start from the current positions and go a full orbit. +# These may be unrolled so that the Earth is in the middle, the following +# is a series of index manipulations to find the correct position for that. + +longitude_b = elem_b.mean_anomaly + elem_b.peri_arg + elem_b.lon_of_ascending +longitude_a = elem_a.mean_anomaly + elem_a.peri_arg + elem_a.lon_of_ascending +i = np.digitize(((longitude_a - longitude_b) % 360) / 360, sigma_steps) + +# here are the unrolled values, which may be plotted +vals_unrolled = np.roll(energy_vals, i + n_sigma // 2) +x_unrolled = np.linspace(-180, 180, n_sigma + 1)[:-1] + +# plot the energy curve +plt.plot(x_unrolled, vals_unrolled) + +# plot the current position of the object +plt.scatter(x_unrolled[(i + n_sigma // 2) % (n_sigma)], energy_vals[0]) + +# plot the current energy of the object +plt.axhline(cur_energy, c="r") +plt.axvline(0, c="k") + +# We can then find and plot the inflection points +local_max = np.argwhere(np.diff(np.sign(np.diff(vals_unrolled))) < 0).ravel() + 1 +for x in vals_unrolled[local_max]: + plt.axhline(x, c="k", ls="--") + + +plt.xlabel(r"$\sigma$ (deg)") +plt.ylabel(r"R($\sigma)$") +plt.title(f"{elem_a.desig}") +plt.show() diff --git a/src/kete/horizons.py b/src/kete/horizons.py index 545cbb0..f49abf2 100644 --- a/src/kete/horizons.py +++ b/src/kete/horizons.py @@ -91,20 +91,15 @@ def fetch(name, update_name=True, cache=True, update_cache=False, exact_name=Fal ] if len(val) == 0: continue - if kete_v == "peri_time": - phys[kete_v] = Time(val[0], scaling="utc").jd - else: - phys[kete_v] = val[0] + phys[kete_v] = val[0] - phys["epoch"] = Time(float(props["orbit"]["epoch"]), scaling="utc").jd + phys["epoch"] = float(props["orbit"]["epoch"]) if "moid" in props["orbit"] and props["orbit"]["moid"] is not None: phys["moid"] = float(props["orbit"]["moid"]) if "data_arc" in props["orbit"] and props["orbit"]["data_arc"] is not None: phys["arc_len"] = float(props["orbit"]["data_arc"]) if props["orbit"]["covariance"] is not None: - cov_epoch = Time( - float(props["orbit"]["covariance"]["epoch"]), scaling="utc" - ).jd + cov_epoch = float(props["orbit"]["covariance"]["epoch"]) mat = np.nan_to_num( np.array(props["orbit"]["covariance"]["data"], dtype=float) ) @@ -118,7 +113,6 @@ def fetch(name, update_name=True, cache=True, update_cache=False, exact_name=Fal elements = { lab: phys[lookup_rev[lab]] for lab in labels if lab in lookup_rev } - elements["tp"] = Time(elements.get("tp"), scaling="utc").jd params = [(lookup_rev.get(x, x), elements.get(x, np.nan)) for x in labels] phys["covariance"] = Covariance(name, cov_epoch, params, mat) else: