Animation showing the changes in orbit for a spacecraft performing aerobraking manoeuvres.

Animation_Aerobraking.m

+%% Function Description (Aerobraking Orbit Animation)
+% The script simulates and animates the process of aerobraking for a spacecraft in an elliptical orbit around Mars. 
+% Aerobraking is a technique used to reduce the velocity and altitude of a spacecraft by utilizing atmospheric drag.
+% This script models the gradual change in orbital parameters, specifically the eccentricity, as the spacecraft performs multiple passes through the Martian atmosphere.
+%
+% The script begins by defining the initial orbital parameters such as the periapsis altitude, the initial eccentricity,
+% and the physical properties of the spacecraft and Mars. Using these parameters, it calculates the initial semi-major axis
+% and orbital period of the spacecraft's orbit.
+%
+% The animation is generated by iteratively solving Kepler's equation to determine the spacecraft's position over time in its orbit.
+% As the spacecraft passes through the Martian atmosphere, the drag force reduces its velocity, causing the orbit to become more circular (reducing eccentricity).
+% This is achieved by updating the semi-major axis and eccentricity of the orbit after each pass.
+%
+% Key concepts and calculations:
+% - Semi-major axis (a): Determines the size of the orbit and is recalculated based on the changing eccentricity.
+% - Eccentricity (e): Measures the shape of the orbit. The aerobraking effect reduces the eccentricity over multiple orbits.
+% - Gravitational parameter (mu): Represents the gravitational influence of Mars, affecting the orbital motion.
+% - Drag force (Fd): Calculated based on the atmospheric density, drag coefficient, cross-sectional area, and velocity of the spacecraft.
+%   This force acts opposite to the velocity, reducing the orbital energy and altering the orbit.
+% - Kepler's equation: Used to solve for the eccentric anomaly, which helps in determining the spacecraft's position in its orbit.
+%
+% The script continuously updates the plot to animate the changing orbit, showing the effect of atmospheric drag on the spacecraft's path around Mars.
+% It also calculates the total time spent in aerobraking and the number of orbits required to reach the target eccentricity.
+%
+% Author: Thomas West
+% Date: September 3, 2024
+
+%% Script
+% Clear all variables from the workspace to ensure a fresh start
+clearvars;
+% Close all open figures to clear the plotting space
+close all;
+% Clear the command window for a clean output display
+clc;
+% Set the display format to long with general number formatting for precision
+format longG;
+
+Mars_radius = 3396;  % Radius of Mars in km
+Altitude = 100;
+Apo_Radius = 400;
+ae = Apo_Radius + Mars_radius;
+pe = Altitude + Mars_radius;
+target_eccentricity = (ae - pe) / (ae + pe);  % Target eccentricity to achieve
+
+% Display the result
+fprintf('The final target orbit eccentricity is: %.6f\n', target_eccentricity);
+
+% Constants and Initial Conditions
+Mars_radius = 3396;  % Radius of Mars in km
+rp = Mars_radius + Altitude;  % Periapsis distance from center of Mars, in km
+e = 0.9;  % Initial Eccentricity
+mu = 42828;  % Mars' gravitational parameter in km^3/s^2
+Cd = 1.75;  % Drag coefficient
+A = 15.9;    % Cross-sectional area of the spacecraft in m^2
+m = 1500;  % Mass of the spacecraft in kg
+H = 11.1;  % Scale height of Mars' atmosphere in km
+rho0 = 0.02;  % Reference density at a reference altitude in kg/m^3
+h0 = Altitude;  % Reference altitude above Mars' surface in km
+
+% Calculate initial semi-major axis and orbital period
+a = calculate_semi_major_axis(rp, e);
+P = 2 * pi * sqrt(a^3 / mu);  % Initial orbital period in seconds
+total_time = 0;  % Initialize total time counter
+
+% Setup figure for plotting
+figure;
+hold on;
+[Xe, Ye, Ze] = sphere(50);
+surf(Xe * Mars_radius, Ye * Mars_radius, Ze * Mars_radius, 'EdgeColor', 'none', 'FaceColor', 'red', 'DisplayName', 'Mars');
+xlabel('X (km)');
+ylabel('Y (km)');
+zlabel('Z (km)');
+title('Orbit Adjustment via Aerobraking');
+axis equal;
+grid on;
+view(0, 90);  % Set the view to top-down
+
+% Calculate initial semi-major axis based on rp and eccentricity
+initial_a = rp / (1 - e);
+% Set fixed axes limits based on the initial semi-major axis
+axis_limit = ceil(a + Mars_radius) * 1.1;  % Increase by 10% for margin
+xlim([-20000, 40000]);
+ylim([-5000, 6.5e4]);
+zlim([-axis_limit, axis_limit]);
+
+% Simulation parameters
+orbit_plots = gobjects(0);  % Store plot handles
+orbit_count = 0;
+total_time = 0;
+
+% Orbital parameters
+i = deg2rad(167.90);  % Inclination in radians
+Omega = deg2rad(276.56);  % Right Ascension of Ascending Node (RAAN) in radians
+omega = deg2rad(118.94 + 270);  % Argument of Perigee in radians
+
+% Create a single marker for the spacecraft
+hMarker = plot3(NaN, NaN, NaN, 'ko', 'MarkerFaceColor', 'k');
+
+while e > target_eccentricity
+    orbit_count = orbit_count + 1;
+    
+    % Calculate the orbit for plotting
+    a = rp / (1 - e);  % Update semi-major axis based on current eccentricity
+    P = 2 * pi * sqrt(a^3 / mu);  % Orbital period
+    total_time = total_time + P;  % Add period to total time
+    M_vec = linspace(0, 2 * pi, 1000);  % Mean anomaly vector for full orbit
+    E_vec = M_vec;  % Initialize eccentric anomaly vector
+
+    % Solve Kepler's equation E - e*sin(E) = M for each point
+    for j = 1:length(M_vec)
+        E_vec(j) = kepler(M_vec(j), e);
+    end
+
+    % Convert eccentric anomaly to true anomaly
+    theta_vec = 2 * atan2(sqrt(1 + e) * sin(E_vec / 2), sqrt(1 - e) * cos(E_vec / 2));
+    X_matrix = zeros(3, length(theta_vec));  % Position matrix
+
+    for k = 1:length(theta_vec)
+        % Use provided COE2RV function to calculate position
+        [RIJK, ~] = COE2RV_Elliptic([a, e, i, Omega, omega, theta_vec(k)], mu);
+        X_matrix(:, k) = RIJK;
+    end
+
+    % Marker and trail plot
+    hTrail = plot3(X_matrix(1, 1), X_matrix(2, 1), X_matrix(3, 1), 'm-', 'LineWidth', 1);
+
+    for j = 1:length(theta_vec)
+        set(hTrail, 'XData', X_matrix(1, 1:j), 'YData', X_matrix(2, 1:j), 'ZData', X_matrix(3, 1:j));
+        set(hMarker, 'XData', X_matrix(1, j), 'YData', X_matrix(2, j), 'ZData', X_matrix(3, j));
+        drawnow;
+        pause(0.00001);  % Adjust for slower/faster animation
+    end
+
+    % Keep last 3 orbits
+    orbit_plots = [orbit_plots, hTrail];  % Append new trail
+    if length(orbit_plots) > 20  % Maintain only the last n orbits
+        delete(orbit_plots(1));  % Remove the oldest orbit plot
+        orbit_plots(1) = [];  % Remove handle from array
+    end
+
+    % Aerobraking effects: update eccentricity for the next loop iteration
+    vp = sqrt(mu * (2 / rp - 1 / a));
+    h = rp - Mars_radius;  % Altitude above Mars' surface
+    rho = rho0 * exp(-(h - h0) / H);
+    Fd = -0.5 * Cd * rho * vp^2 * A;
+    delta_v = Fd / m;
+    vp_new = vp + delta_v;
+    a_new = mu / (2 * mu / rp - vp_new^2);
+    e = 1 - rp / a_new;  % Update eccentricity
+end
+
+fprintf('Total time spent in aerobraking: %.2f seconds (%.2f days).\n', total_time, total_time / 86400);
+fprintf('Total number of orbits to achieve target eccentricity: %d\n', orbit_count);
+hold off;
+
+% Kepler's equation solver using Newton's method
+function E = kepler(M, e)
+    E = M;  % Initial guess
+    for i = 1:100  % Maximum 100 iterations
+        f = E - e * sin(E) - M;
+        f_prime = 1 - e * cos(E);
+        E = E - f / f_prime;  % Newton's method update
+        if abs(f) < 1e-10  % Convergence criterion
+            break;
+        end
+    end
+end
\ No newline at end of file