How to implement the Central Difference method integrator for a two body problem?

Considering I have the perigee altitude (measured from the center of the earth), entry tangential velocity at the periapsis, and the orbital inclination is 0, how do i implement this into MatLab to calculate the numerical orbit time period. Define a stopping criterion for the time integration loop in MatLab and save the time variable upon meeting the stopping criterion.

Answer 1

[Balavignesh]

Hi Sathvik,

It is my understanding that you would like to calculate the numerical orbit time period for a given set of input values. One effective approach is to integrate the orbit's equations of motion using a numerical method, such as the Runge-Kutta method.

In MATLAB, the ode45 function is a versatile solver for ordinary differential equations. For this example, we'll employ the equations of the two-body problem in polar coordinates, assuming the orbit is either circular or elliptical (without perturbations) and that Earth is the central body.

The initial step involves calculating the semi-major axis and the initial angular momentum, followed by setting up the differential equations for solving. You can utilize ode45 to perform the integration over time and define a stopping criterion for the integration process.

Below is an example MATLAB code that illustrates this concept:

% Constants
mu = 3.986004418e5; % Earth's gravitational parameter in km^3/s^2

% Inputs
r_perigee = 6571; % Example: perigee altitude from the center of the Earth in km (Earth's radius + 200 km)
v_perigee = 7.8; % Example: tangential velocity at perigee in km/s

% Initial conditions
r0 = r_perigee; % Initial radius
v0 = v_perigee; % Initial tangential velocity

% Equations of motion for the two-body problem in polar coordinates
% dr/dt = vr
% dvr/dt = vtheta^2 / r - mu / r^2
% dtheta/dt = vtheta / r
% dvtheta/dt = - vr * vtheta / r
% Where vr is radial velocity (initially 0 for perigee) and vtheta is tangential velocity

% Define the ODE system
odefun = @(t, y) [y(2); y(4)^2 / y(1) - mu / y(1)^2; y(4) / y(1); -y(2) * y(4) / y(1)];

% Initial state vector: [r; vr; theta; vtheta]
y0 = [r0; 0; 0; v0];

% Time span: start at 0, we don't define an end time, letting the event function handle stopping
tspan = [0 Inf];

% Define an event function to stop the integration at one orbit (theta reaches 2*pi)
options = odeset('Events', @(t, y) event_OneOrbit(t, y));

% Solve the ODE
[t, y, te, ye, ie] = ode45(odefun, tspan, y0, options);

% The orbit time period is the time at the event
orbitTimePeriod = te;

% Display the result
disp(['Orbit Time Period: ', num2str(orbitTimePeriod), ' seconds']);

Orbit Time Period: 5324.6091 seconds

% Event function to stop integration when theta reaches 2*pi
function [value,isterminal,direction] = event_OneOrbit(t, y)
    value = y(3) - 2*pi; % Detect theta - 2*pi
    isterminal = 1; % Stop the integration
    direction = 1; % Positive direction only
end

Kindly have a look at the following documentation links to have more information on:

• 'odeset' function: https://www.mathworks.com/help/matlab/ref/odeset.html
• 'ode45' function: https://www.mathworks.com/help/matlab/ref/ode45.html

Hope this helps!

Balavignesh

Answer 2

[James Tursa]

You have:

altitude

tangential velocity

If you just want the orbital period, you can calculate as follows using the Vis-Viva equation and Kepler's 3rd Law:

Re = 6378137; % Equatorial radius of Earth (m)
mu = 3.986e14; % Gravitational parameter of Earth (m^3/s^2)
alt = _____ % Altitude of satellite (m)  <-- You fill this in
vt = _____ % Tangential velocity of satellite (m/s)  <--- You fill this in
r = Re + alt; % Radius of satellite (m)
a = 1 / (2/r - vt^2/mu); % Vis-Viva solved for semi-major axis (m)
n = sqrt(mu/a^3); % Kepler's 3rd Law solved for mean motion (rad/s)
P = 2*pi/n; % Period (s)