I understand that your code is getting stuck while solving the function using the ode45 solver.
To address the issue, please follow these suggestions:
- Using the gravitational constant (G) and the Earth's mass (M) separately can lead to computational inefficiencies and precision errors. Instead, replace
with the gravitational parameter
for a more accurate and simplified calculation.
% Gravitational parameter and Earth's radius
mu = 398600; % Gravitational parameter in km^3/s^2
R = 6378; % Earth's radius in km
% Initial conditions
r0 = [3207, 5459, 2714]; % Initial position vector in km
v0 = [-6.532, 0.7835, 6.142]; % Initial velocity vector in km/s
t0 = 0; % Initial time in seconds
- Modify the odefunc to use μ directly.
% Define the ODE function for orbital motion
odeFunc = @(t, y) [
y(4);
y(5);
y(6);
-mu * y(1) / (norm(y(1:3))^3);
-mu * y(2) / (norm(y(1:3))^3);
-mu * y(3) / (norm(y(1:3))^3)
];
% Set options for the ode45 solver
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-6);
- The initial conditions vectors were concatenated using a semicolon [r0; v0], causing the solver to misinterpret the input vector. Change [r0; v0] to [r0, v0] to ensure the initial state vector is correctly formatted as a single row, compatible with ode45.
% Solve the ODEs numerically over the specified time span
[t, y] = ode45(odeFunc, [t0, 1.66 * 3600], [r0, v0], options);
% Calculate altitude above Earth's surface
altitudes = sqrt(sum(y(:, 1:3).^2, 2)) - R;
% Find the maximum altitude and the time at which it occurs
[maxAltitude, maxAltitudeIndex] = max(altitudes);
timeOfMaxAltitude = t(maxAltitudeIndex);
Applying these changes will ensure that the implementation is accurate as demonstrated here:
% Display the results
fprintf('Maximum Altitude: %.2f km\n', maxAltitude);
fprintf('Time of Maximum Altitude: %.2f hours\n', timeOfMaxAltitude / 3600);
Refer to the following documentation for more details about the function:
Hope this helps.
