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

8 ビュー (過去 30 日間)
Sathvik 2024 年 2 月 27 日

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.

サインインしてコメントする。

### 回答 (2 件)

Balavignesh 2024 年 3 月 15 日
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:
Hope this helps!
Balavignesh
##### 0 件のコメント-2 件の古いコメントを表示-2 件の古いコメントを非表示

サインインしてコメントする。

James Tursa 2024 年 3 月 19 日

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)

サインインしてコメントする。

### カテゴリ

Help Center および File ExchangeCubeSat and Satellites についてさらに検索

### Community Treasure Hunt

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

Start Hunting!

Translated by