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

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

4 views (last 30 days)

Show older comments

##### 0 Comments

### Answers (2)

Balavignesh
on 15 Mar 2024

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']);

% 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:

Hope this helps!

Balavignesh

##### 0 Comments

James Tursa
on 19 Mar 2024

Edited: James Tursa
on 19 Mar 2024

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)

##### 0 Comments

### See Also

### Categories

### Products

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!