Project Didymos
Earth to NEO 65803 Preliminary Analysis
Porkchop Plot: 1000-Day Launch Window Iteration
Objective
To find the optimal transfer time to reach 65803 Didymos, a Near-Earth Object (NEO), in a 1000 day interval.
The provided orbital elements of Earth and Didymos at some reference time (tref) are:
Orbital Elements {a, e, i, Ω, ω, θ(tref)}
- 1 AU = 149.6 * 106 km
- 1 day = 86400s
- μ = 132712 * 106 km3/s2
65803 Didymos (Credit: Wikipedia)
Original Positions & OrbitAssumptions
- Zero-sphere of influence at both Earth and Didymos
- Restricted two body equations are valid (Justified by the fact that MSun >> MDidymos and MSun >> MEarth)
- No perturbations (Negligible drag, J2 perturbations, 3rd bodies, etc)
- No orbital decay (The only Keplerian Element that changes over time is true anomaly)
Why MATLAB?
MATLAB provides the ode78 suite, which is essential for high-order propagation required in orbital mechanics. Implementing a comparable variable time-step ODE solver in a lower-level language is a significant undertaking on its own.
Furthermore, the native graphing and data visualization tools in MATLAB are significantly more streamlined than system-level alternatives like C++. This allows for rapid iteration of complex Porkchop plots and 3D trajectory renders, which were critical for the preliminary analysis phase of this project.
Restricted 2-Body Problem (R2BP)
Every body is pulled by every other body, but the N-body problem has no analytical solution and is expensive to integrate. I used the restricted two-body problem instead, where the Sun is the only attracting mass and the orbiting body (Earth, Didymos or the spacecraft) is too light to move it; that holds here because MSun is far larger than MEarth or MDidymos, and it is what gives Kepler’s problem below a semi-analytical solution.
The situation I applied the R2BP to(Credit: Wiki)
Keplerian Element Pic(Credit: Sun, He & Bouman, Katherine & Tiede, Paul & Wang, et al (2022))
Orbital Elements
A state in the R2BP has six degrees of freedom, written either as a Cartesian position and velocity vector or as the Keplerian elements {a, e, i, Ω, ω, ν}in the perifocal frame. With no perturbations, true anomaly (ν) is the only Keplerian element that changes with time, which is why I propagate in elements, and converting between the two forms is a closed-form rotation (“kep2rv.m” and “rv2kep.m”).
Kepler’s Problem
Kepler’s problem finds the state after a known time of flight from a known starting state. The only numerical step is Kepler’s equation, M = E - e sin(E), which I solve for the eccentric anomaly E with Newton-Raphson to a tolerance of 10-12rad; the rest is algebra, so it costs less than stepping the orbit forward in time and does not build up the error a time-step integrator gathers over a long time of flight (“propagate_elliptical_ta.m” and “propagate_elliptical.m”).
Kep Problem Reference Circle
Position Boundary Value ProblemLambert Solver / Gauss Problem
The cost of a transfer is its ΔV, the velocity change that takes the spacecraft from its current orbit onto the transfer orbit plus the change that takes it from the transfer orbit onto the final orbit, and the rocket equation turns that into propellant mass for any given spacecraft. Getting both velocities means finding the orbit that joins two known positions in a set time of flight, which is Lambert’s (Gauss’s) boundary value problem; I solve it with a universal-variable Lambert solver (“lambert_universal.m”), which is semi-analytical like Kepler’s problem and can fail when the two position vectors are linearly dependent.
How the Analysis Fits Together
For each departure and arrival day, I propagate Earth to the departure day and Didymos to the arrival day with Kepler’s problem, solve Lambert’s problem between the two positions for that time of flight, and add the departure and arrival ΔV. Repeating that over the whole window gives the porkchop plot, and I check a chosen transfer by integrating Earth, Didymos and the transfer orbit about the Sun with ode78 and a two-body right-hand side (“rhs_2body.m”).
Porkchop Plot (Prograde) - [tref, tref + 1000 days]
Specific Time Interval
Analysis Window:
- t₀ ≥ t_ref
- t_f ≤ t_ref + 1000 days
- [t_ref, t_ref + 1000 days]
Mission Results
- Departure (t₀): t_ref + 600 days
- Arrival (t_f): t_ref + 1000 days
- TOF: 400 days
- ΔV at departure: 3.47 km/s
- ΔV at arrival: 4.39 km/s
- Local Optimal ΔV 7.86 km/s
Constraints & Parameters:
- I skip any TOF under 100 days to cut computing time, since transfers that short need too much ΔV anyway.
- Data calculated at 5-day intervals for optimization.
Optimal Earth to Didymos Transfer TrajectoryWhat the Answer Means
The 7.86 km/s adds two velocity changes, 3.47 km/s to leave Earth’s orbital velocity for the transfer orbit and 4.39 km/s to match Didymos’s velocity at arrival, so it is the cost of a prograde rendezvous, not a flyby. With a zero sphere of influence, both changes are measured against each body’s velocity around the Sun, so the number leaves out escaping Earth’s gravity and is not the ΔV needed from a parking orbit. The best arrival, tref + 1000 days, is the last day I searched, so 7.86 km/s is the lowest value inside the window, not a proven minimum, and on the 5-day grid the true minimum can sit up to 5 days from the departure and arrival days I report.