Specification for Kepler Orbit Propagation
1. Overview
1. Functions
- The
KeplerOrbit
class calculates the satellite position and velocity with the simple two-body problem. We ignored any disturbances and generated forces by the satellite. - This orbit propagation mode provides the simplest and fastest orbit calculation for any orbit(LEO, GEO, Deep Space, and so on.).
2. Files
src/dynamics/orbit/orbit.hpp, cpp
- Definition of
Orbit
base class
- Definition of
src/dynamics/orbit/initialize_orbit.hpp, .cpp
- Make an instance of orbit class.
src/dynamics/orbit/kepler_orbit_propagation.hpp, .cpp
- Libraries
src/library/orbit/kepler_orbit.cpp, hpp
src/library/orbit/orbital_elements.cpp, hpp
3. How to use
- Select
propagate_mode = KEPLER
in the spacecraft's ini file. - Select
initialize_mode
as you want.DEFAULT
: Use default initialize method (RK4
andENCKE
use position and velocity,KEPLER
uses init_mode_kepler)POSITION_VELOCITY_I
: Initialize with position and velocity in the inertial frameORBITAL_ELEMENTS
: Initialize with orbital elements
2. Explanation of Algorithm
1. KeplerOrbit::CalcConstKeplerMotion
function
1. Overview
- This function calculates the following constant values.
- Mean motion
- Frame conversion matrix from in-plane position to the inertial frame position
2. Inputs and outputs
- Input
- Orbital Elements
- $\mu$ : The standard gravitational parameter of the central body
- Output
- $n$ : Mean motion
- $R_{p2i}$ : Frame conversion matrix from in-plane position to the inertial frame position
3. Algorithm
- Mean motion n=√μa3
- Frame conversion matrix from in-plane position to ECI position Rp2i=Rz(−Ω)Rx(−i)Rz(−ω) Rx(θ)=(1000cosθsinθ0−sinθcosθ) Rz(θ)=(cosθsinθ0−sinθcosθ0001)
2. KeplerOrbit::CalcOrbit
function
1. Overview
- This function calculates the position and velocity of the spacecraft in the inertial frame at the designated time.
2. Inputs and outputs
- Input
- $t$ : Time in Julian day
- Orbital Elements
- Constants
- Output
- $\boldsymbol{r}_{i}$ : Position in the inertial frame
- $\boldsymbol{v}_{i}$ : Velocity in the inertial frame
3. Algorithm
- Calculate mean anomaly $l$[rad] l=n∗(t−tepoch)
- Calculate eccentric anomaly $u$[rad] by solving the Kepler Equation
- Details are described in
KeplerOrbit::SolveKeplerFirstOrder
- Details are described in
- Calculate two dimensional position $x^{*}$, $y^{*}$ and velocity $\dot{x}^{*}$, $\dot{y}^{*}$ in the orbital plane x∗=a(cosu−e)y∗=a√1−e2sinu˙x∗=−ansinu1ecosu˙y∗=an√1−e2cosu1ecosu
- Convert to the inertial frame ri=Rp2i[x∗y∗0]vi=Rp2i[˙x∗˙y∗0]
3. KeplerOrbit::SolveKeplerFirstOrder
function
1. Overview
- This function solves the Kepler Equation with the first-order iterative method.
- Note: This method is not suited to the high eccentricity orbit. It is better to use the Newton-Raphson method for such a case.
- From
v6.0.0
we haveSolveKeplerNewtonMethod
and use it. The detail of the method will be write.
- From
2. Inputs and outputs
- Input
- $e$ : eccentricity
- $l$ : mean anomaly [rad]
- $\epsilon$ : threshold for convergence [rad]
- Limit of iteration
- Output
- $u$ : eccentric anomaly [rad]
3. Algorithm
- Set the initial value of eccentric anomaly as follows
- $u_0=l$
- Calculate $u_{n+1}$ with the following equation un+1=l+esinun
- Iterate the calculation until the following conditions are satisfied
- $|u_{n+1} - u_{n}| < \epsilon$
- The iteration number over the limit of iteration
4. OrbitalElements::CalcOeFromPosVel
function
1. Overview
- This function calculates the orbital elements from the initial position and velocity in the inertial frame.
2. Inputs and outputs
- Input
- $\mu$ : The standard gravitational parameter of the central body
- $t$ : Time in Julian day
- $\boldsymbol{r}_{i}$ : Initial position in the inertial frame
- $\boldsymbol{v}_{i}$ : Initial velocity in the inertial frame
- Output
- orbital element
3. Algorithm
- $\boldsymbol{h}_{i}$ : Angular momentum vector of the orbit hi=ri×vi
- $a$ : Semi-major axis a=μ2μr−v2
- $i$ : Inclination i=cos−1hz
- $\Omega$ : Right Ascension of the Ascending Node (RAAN)
- Note: This equation is not support $i = 0$ case. Ω=sin−1(hx√h2x+h2y)
- $x_{p}, y_{p}$ : Position in the orbital plane xp=rixcosΩ+riysinΩ;yp=(−rixsinΩ+riycosΩ)cosi+rizsini
- $\dot{x}_{p}, \dot{y}_{p}$ : Velocity in the orbital plane ˙xp=vixcosΩ+viysinΩ;˙yp=(−vixsinΩ+viycosΩ)cosi+vizsini
- $e$ : Eccentricity c1=hμ˙yp−xprc2=−hμ˙xp−ypre=√c21+c22
- $\omega$ : Argument of Perigee ω=tan−1(c2c1)
- $t_{epoch}$ : Epoch [Julian day] f=tan−1(ypxp)−ωu=tan−1rsinf√1−e2rcosf+ae \[ \begin{align} n &= \sqrt{\frac{\mu}{a^3}}\\ dt &= \frac{u - e\sin{u}}{n}\\ t_{epoch} &= t - \frac{dt}{246060} \end{align} \]
3. Results of verifications
1. Comparison with RK4
1. Overview
- We compared the calculated orbit result between RK4 mode and Kepler mode. In the RK4 mode, all disturbances are disabled since the Kepler mode ignores them.
- In the Kepler mode, we verified the correctness of both initialize modes (
ORBITAL_ELEMENTS
andPOSITION_VELOCITY_I
).
2. Conditions for the verification
- sample_simulation_base.ini
- The following values are modified from the default.
EndTimeSec = 10000 LogOutPutIntervalSec = 5
- The following values are modified from the default.
- SampleDisturbance.ini
- All disturbances are disabled.
- SampleSat.ini
- The following values are modified from the default.
propagate_mode
is changed for each mode.- Orbital elements for Kepler
semi_major_axis_m = 6794500.0 eccentricity = 0.0015 inclination_rad = 0.9012 raan_rad = 0.1411 arg_perigee_rad = 1.7952 epoch_jday = 2.458940966402607e6
- Initial position and velocity (compatible value with the orbital elements)
init_position(0) = 1791860.131 init_position(1) = 4240666.743 init_position(2) = 4985526.129 init_velocity(0) = -7349.913889 init_velocity(1) = 631.6563971 init_velocity(2) = 2095.780148
- The following values are modified from the default.
3. Results
-
The following figure shows the orbit calculation result of Kepler mode with
ORBITAL_ELEMENTS
initialize mode.- The result looks correct.
-
The difference between Kepler mode with
ORBITAL_ELEMENTS
initialize mode and RK4 mode is shown in the following figure.- The error between them is small (less than 10m), and we confirmed that the calculation of Kepler orbit is correct.
-
The following figure shows the difference between Kepler orbit calculation with
ORBITAL_ELEMENTS
andPOSITION_VELOCITY_I
initialize mode.- The error between them is small (less than 10m), and we confirmed that the initializing method is correct.
4. References
- Hiroshi Kinoshita, "Celestial mechanisms and orbital dynamics", ch.2, 1998 (written in Japanese)