Ideal Two Body Problem

The study of the ideal two-body problem forms the core of our codebase, serving as the foundation for more complex algorithms and calculations in orbital dynamics. This module introduces two key types: the state vector and the coe (classical orbital elements); and provides a comprehensive set of functions that enable precise analysis and prediction of celestial motion.

The ideal two-body problem focuses on the motion of two point masses, assuming no external influences from other celestial bodies. Although this simplification may appear trivial, it serves as a fundamental building block for understanding and modeling more intricate orbital dynamics scenarios.

Within the idealTwoBodyProblem.jl file, a range of functions has been developed to address various aspects of the ideal two-body problem. These functions form the backbone of orbital computations, providing essential tools for the calculation of critical parameters and transformations.

The figure below shows the two new data types defined in the ideal two-body problem file: MyStateVector and MyCOE. Below each appears their store objects and their type (indicated as fieldName::Type). These two data types are crucial for astrodynamics, as they allow to completely describe the motion of a body in space.

I2BPTypes

It is important to understand what each of these types entails and how to construct them:

  1. MyStateVector: stores the position vector (r) and the velocity vector (v) of a body in space. Although there is no general imposition on the units, because of the nature of distances in space, the position of the satellite/spacecraft is usually given in km, while the velocity is given in km/s. Providing the coordinates in these units ensures the correct functioning of all the functions in the code. There is only one possible way to construct a MyStateVector object, and it is as follows:
    • MyStateVector(r, v), where both r and v must be MyVector objects.
  1. MyCOE: stores all six classical orbital elements (Ω, i, ω, a, e, θ), and p, so that the orbit can always be well defined. This is because of the peculiarity of the parabolic orbit, where a = ∞. However, for the construction, only six of the elements are given, where the user can choose to either provide the semi-major axis a or the semi-latus rectum p, but never both. However, this seventh element can always be accessed as it is computed internally. It is very important to note that all the angles must be given in radians, which is not optional. Additionally, all the parameters must be of type Float64. Examples for each constructor are provided below:
    • MyCOE(RAAN = 0.18, inc = 0.26, omega = 0.30, a = 12000.0, e = 0.54, theta = π), where a is provided and p is computed internally.
    • MyCOE(RAAN = 0.52, inc = 0.0, omega = 0.92, p = 8000.0, e = 1.0, theta = 2.09), where p is provided and a is computed internally.

This first table presents an overview of the functions available within the ideal two-body problem file. These functions encompass a wide range of functionalities, including the computation of the angular momentum vector, the conversion between state vectors and classical orbital elements, and the propagation of state vectors over time. With a total of 26 functions, this comprehensive set ensures that users have the necessary tools to accurately analyze and simulate the motion of celestial objects.

Function NameInputsOutputsFunctionality
angularMomentumVectorstateVector::MyStateVectorh::MyVectorObtain h vector
mechanicalEnergystateVector::MyStateVector, mu::Float64xi::Float64Obtain ξ
eccentricityVectorstateVector::MyStateVector, mu::Float64e::MyVectorObtain e vector
semilatusRectumstateVector::MyStateVector, mu::Float64p::Float64Obtain p
periodOrbitstateVector::MyStateVector, mu::Float64, OrbitType::Stringperiod::Float64Obtain period τ [s]
perifocalBasisstateVector::MyStateVector, mu::Float64perifocalBasis::MyBasisPerifocal Basis
periapsisstateVector::MyStateVector, mu::Float64r_p::Float64Obtain r_p
apoapsisstateVector::MyStateVector, mu::Float64r_a::Float64Obtain r_a
trajectoryEquationstateVector::MyStateVector, mu::Float64, theta::Float64trajectoryEquation::Float64Obtain r
inclinationstateVector::MyStateVectorinc::Float64Obtain i
semiMajorAxisstateVector::MyStateVector, mu::Float64a::Float64Obtain a
orbitTypestateVector::MyStateVector, mu::Float64orbitType::StringObtain orbit type: Elliptic, Hyperbolic, Parabolic, Circular
trueAnomalystateVector::MyStateVector, mu::Float64theta::Float64Obtain θ [rad]
stateVector_to_COEstateVector::MyStateVector, mu::Float64coe::MyCOE(r, v) → COE
COE_to_stateVectorcoe::MyCOE, mu::Float64stateVector::MyStateVectorCOE → (r, v)
escapeVelocityr::MyVector, mu::Float64v_esc::Float64Obtain v_e
circularVelocity2 OPTIONS: 1.r::MyVector, mu::Float64 2.r_norm::Float64, mu::Float64v_circ::Float64Obtain v_c
velocityPeriapsisstateVector::MyStateVector, mu::Float64v_p::Float64Obtain v_p
velocityApoapsisstateVector::MyStateVector, mu::Float64v_a::Float64Obtain v_a
velocityHyperbolicAsymptotestateVector::MyStateVector, mu::Float64v_h::Float64Obtain v_h
speed_visVivar_norm::Float64, a::Float64, mu::Float64v::Float64Obtain v
rightAscensionOfTheAscendingNodestateVector::MyStateVectorRAAN::Float64Obtain Ω
argumentOfPeriapsisstateVector::MyStateVector, mu::Float64omega::Float64Obtain ω
vector_in_ECI_to_perifocalECI_vector::Vector, coe::MyCOEperifocalVector::MyVectorr_ECI → r_PQW
vector_in_perifocal_to_ECIperifocalVector::Vector, coe::MyCOEECI_vector::MyVectorr_PQW → r_ECI
propagate_StateVectorstateVector::MyStateVector, mu::Float64, delta_t::Float64stateVector2::MyStateVectorPropagate (r, v) by Δt

Additionally, as displayed on the second table, the file incorporates functions specific to Kepler's problem, which focuses on the motion of bodies under the influence of a central gravitational force. The aim of this part is to be able to relate time and position. In order to achieve that, it is necessary to be able to convert between the different anomalies at ease. These conversions differ depending on the orbit type (elliptic, parabolic or hyperbolic), and obtaining them can become quite tricky. Because for the path from the time to the auxiliary anomaly, an iterative numerical method is needed. In our code, the Newton-Raphson method was used, provided by the Roots package.

However, in the code not only the specific cases are developed. We also incorporated the universal formulation of Kepler's equation to solve for the universal anomaly $\chi$ which applied an Stumpff function and the numerical solver.

Function NameInputsOutputsFunctionality
timeSincePeriapsisstateVector::MyStateVector, mu::Float64, theta::Float64time::Float64$\theta$ [rad] $\rightarrow \Delta t$ [s]
timeSinceTrueAnomalystateVector::MyStateVector, mu::Float64, theta1::Float64, theta2::Float64time::Float64$\Delta t$ between $\theta_1$ & $\theta_2$
stumpffFunctionz::Float64S::Float64Stumpff function
universalAnomalystateVector::MyStateVector, mu::Float64, time::Float64X::Float64$\Delta t \rightarrow \chi$
universalAnomaly_to_trueAnomalystateVector::MyStateVector, mu::Float64, X::Float64theta::Float64$\chi \rightarrow \theta$
trueAnomaly_to_universalAnomalystateVector::MyStateVector, mu::Float64, theta::Float64X::Float64$\theta \rightarrow \chi$
timeSincePeriapsis_to_trueAnomalystateVector::MyStateVector, mu::Float64, time::Float64theta::Float64$\Delta t \rightarrow \theta$
trueAnomaly_to_eccentricAnomalystateVector::MyStateVector, mu::Float64, theta::Float64E::Float64$\theta \rightarrow E$
eccentricAnomaly_to_trueAnomalystateVector::MyStateVector, mu::Float64, E::Float64theta::Float64$E \rightarrow \theta$
eccentricAnomaly_to_meanAnomalystateVector::MyStateVector, mu::Float64, E::Float64Me::Float64$E \rightarrow M_e$
meanAnomaly_to_eccentricAnomalystateVector::MyStateVector, mu::Float64, Me::Float64E::Float64$M_e \rightarrow E$
trueAnomaly_to_meanAnomaly_EstateVector::MyStateVector, mu::Float64, theta::Float64Me::Float64$\theta \rightarrow M_e$
meanAnomaly_to_trueAnomaly_EstateVector::MyStateVector, mu::Float64, Me::Float64theta::Float64$M_e \rightarrow \theta$
universalAnomaly_to_eccentricAnomalystateVector::MyStateVector, mu::Float64, X::Float64E::Float64$\chi \rightarrow E$
eccentricAnomaly_to_universalAnomalystateVector::MyStateVector, mu::Float64, E::Float64X::Float64$E \rightarrow \chi$
timeSincePeriapsis_to_MestateVector::MyStateVector, mu::Float64, time::Float64Me::Float64$\Delta t \rightarrow M_e$
Me_to_timeSincePeriapsisstateVector::MyStateVector, mu::Float64, Me::Float64time::Float64$M_e \rightarrow \Delta t$
trueAnomaly_to_parabolicAnomalystateVector::MyStateVector, mu::Float64, theta::Float64P::Float64$\theta \rightarrow P$
parabolicAnomaly_to_trueAnomalystateVector::MyStateVector, mu::Float64, P::Float64theta::Float64$P \rightarrow \theta$
parabolicAnomaly_to_meanAnomalystateVector::MyStateVector, mu::Float64, P::Float64Mp::Float64$P \rightarrow M_p$
meanAnomaly_to_parabolicAnomalystateVector::MyStateVector, mu::Float64, Mp::Float64P::Float64$M_p \rightarrow P$
trueAnomaly_to_meanAnomaly_PstateVector::MyStateVector, mu::Float64, theta::Float64Mp::Float64$\theta \rightarrow M_p$
meanAnomaly_to_trueAnomaly_PstateVector::MyStateVector, mu::Float64, Mp::Float64theta::Float64$M_p \rightarrow \theta$
universalAnomaly_to_hyperbolicAnomalystateVector::MyStateVector, mu::Float64, X::Float64H::Float64$\chi \rightarrow H$
hyperbolicAnomaly_to_universalAnomalystateVector::MyStateVector, mu::Float64, H::Float64X::Float64$H \rightarrow \chi$
timeSincePeriapsis_to_MhstateVector::MyStateVector, mu::Float64, time::Float64Mh::Float64$\Delta t \rightarrow M_h$
Mh_to_timeSincePeriapsisstateVector::MyStateVector, mu::Float64, Mh::Float64time::Float64$M_h \rightarrow \Delta t$
trueAnomaly_to_meanAnomaly_HstateVector::MyStateVector, mu::Float64, theta::Float64Mh::Float64$\theta \rightarrow M_h$
meanAnomaly_to_trueAnomaly_HstateVector::MyStateVector, mu::Float64, Mh::Float64theta::Float64$M_h \rightarrow \theta$
universalAnomaly_to_hyperbolicAnomalystateVector::MyStateVector, mu::Float64, X::Float64H::Float64$\chi \rightarrow H$
hyperbolicAnomaly_to_universalAnomalystateVector::MyStateVector, mu::Float64, H::Float64X::Float64$H \rightarrow \chi$
timeSincePeriapsis_to_MhstateVector::MyStateVector, mu::Float64, time::Float64Mh::Float64$\Delta t \rightarrow M_h$
Mh_to_timeSincePeriapsisstateVector::MyStateVector, mu::Float64, Mh::Float64time::Float64$M_h \rightarrow \Delta t$

To better understand all the possible directions of conversion between the different anomalies and time, the diagram below can be consulted.

AnomaliesDiagram

Inside the file, the initial orbit determination is also tackled through the Lambert's problem. The functions written in this third table facilitate the determination of the time of flight, conic section parameters, and the solution of Lambert's problem itself.

Function NameInputsOutputsFunctionality
Lambert_solver1::MyVector, r2::MyVector, tF::Float64, mu::Float64, k::MyVectorcoe1::MyCOE, coe2::MyCOEGiven two positions $\mathbf{r_1}$ and $\mathbf{r_2}$ and the time of flight between them, find the orbit that connects them
Lambert_conicr1::MyVector, r2::MyVector, eT::Float64, k::MyVectorcoe1::MyCOE, coe2::MyCOEObtain the conic section for a given transverse eccentricity $e_T$ using the Avanzini algorithm
timeOfFlightcoe1::MyCOE, coe2::MyCOE, mu::Float64, theta2::Float64, mu::Float64delta_t::Float64Time of flight between two true anomalies ($\theta_1$ and $\theta_2$)