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