diff --git a/Animation_Aerobraking.m b/Animation_Aerobraking.m new file mode 100644 index 0000000000000000000000000000000000000000..7db3c20193a1ce67c39d81ee4e8a9600e81ff2e6 --- /dev/null +++ b/Animation_Aerobraking.m @@ -0,0 +1,167 @@ +%% 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