From 88f3a68c73d1987f76f6f5cbad2e4da248438ea5 Mon Sep 17 00:00:00 2001
From: Thomas West <tw00956@surrey.ac.uk>
Date: Tue, 3 Sep 2024 16:24:06 +0000
Subject: [PATCH] MATLAB code for aerobraking calculations.
---
Aerobraking.m | 271 ++++++++++++++++++++++++++++++++++++++++++++++++++
1 file changed, 271 insertions(+)
create mode 100644 Aerobraking.m
diff --git a/Aerobraking.m b/Aerobraking.m
new file mode 100644
index 0000000..0ce7075
--- /dev/null
+++ b/Aerobraking.m
@@ -0,0 +1,271 @@
+%% Function Description (Aerobraking Calculation for Mars Orbit)
+% This script calculates the aerobraking maneuver required to adjust a spacecraft's orbit around Mars.
+% Aerobraking is a technique used to reduce the apoapsis of an orbit by passing through the atmosphere of
+% a planet, where drag forces slow down the spacecraft, causing it to lose energy and gradually lower its orbit.
+
+% The script begins by clearing the workspace, closing figures, and setting the output format to ensure
+% a fresh environment. It defines constants and initial conditions for the spacecraft's orbit around Mars,
+% such as the semi-major axis, eccentricity, gravitational parameter of Mars, and atmospheric properties
+% needed for drag calculations.
+
+% Key steps and calculations in this script are:
+% - Defining the initial orbital parameters of the spacecraft, including periapsis distance,
+% initial eccentricity, drag coefficient, cross-sectional area, and mass.
+% - Calculating the target eccentricity and semi-major axis for the orbit after aerobraking using the
+% desired apoapsis and periapsis distances.
+% - Setting up a simulation loop to iteratively compute the effects of atmospheric drag on the spacecraft
+% during each pass through Mars' atmosphere. The drag force is calculated using the atmospheric density
+% model (exponential atmosphere) and the drag equation.
+% - Updating the spacecraft's velocity, semi-major axis, and eccentricity after each orbit to account for
+% the effects of drag. The loop continues until the target eccentricity is achieved or the apoapsis
+% radius is reduced below a specified limit.
+% - Propagating the orbit using numerical integration (ode45) to simulate the spacecraft's path around Mars
+% and calculate the position and velocity vectors over time. The output state vectors are used to plot
+% the spacecraft's trajectory and visualize the orbit adjustment process.
+% - Plotting the results, including Mars' position, the spacecraft's trajectory, and other relevant elements,
+% to provide a visual representation of the aerobraking maneuver and the final orbit around Mars.
+
+% The script is designed to be flexible, allowing for adjustments to various parameters, such as initial
+% orbital elements, atmospheric properties, and drag coefficients. It provides a comprehensive approach
+% to computing aerobraking maneuvers for spacecraft around Mars, taking into account the atmospheric
+% drag effects and orbital dynamics.
+
+% 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;
+
+% Inclination, RAAN, and Argument of Periapsis (converted to radians)
+i = deg2rad(167.90); % Inclination in radians
+Omega = deg2rad(276.56); % Right Ascension of the Ascending Node (RAAN) in radians
+omega = deg2rad(118.94 + 270); % Argument of periapsis in radians
+
+Mars_radius = 3396; % Radius of Mars in km
+Altitude = 150;
+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
+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 chrome-extension://efaidnbmnnnibpcajpcglclefindmkaj/https://ntrs.nasa.gov/api/citations/20210023957/downloads/Mars-GRAM%20User%20Guide.pdf
+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;
+xlim([-2e4, 5.5e4]); % Set limits for the x-axis
+ylim([-1e4, 7e4]); % Set limits for the y-axis
+
+% Initialize variables
+orbit_count = 0;
+total_time = 0; % Total time in seconds
+
+% Simulation loop
+while e > target_eccentricity
+ orbit_count = orbit_count + 1;
+
+ % Update velocity and orbital elements
+ vp = sqrt(mu * (2 / rp - 1 / (rp / (1 - e))));
+ 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 = mu / (2 * mu / rp - vp_new^2);
+ e = 1 - rp / a; % Update eccentricity
+
+ % Recalculate orbital period and update total time
+ P = 2 * pi * sqrt(a^3 / mu);
+ total_time = total_time + P; % Update total time
+
+ % Plot initial orbit and every 20 orbits or the final orbit
+ if orbit_count == 1 || mod(orbit_count, 20) == 0 || e <= target_eccentricity
+ t_vec = linspace(0, P, 1000);
+ theta_vec = linspace(0, 2 * pi, length(t_vec));
+ X_matrix = zeros(6, length(t_vec));
+ for k = 1:length(t_vec)
+ [RIJK, VIJK] = COE2RV_Elliptic([a, e, i, Omega, omega, theta_vec(k)], mu);
+ X_matrix(:, k) = [RIJK; VIJK]; % Storing position and velocity vectors
+ end
+ days_elapsed = total_time / 86400; % Convert seconds to days
+ plot3(X_matrix(1, :), X_matrix(2, :), X_matrix(3, :), 'LineWidth', 1.5, ...
+ 'DisplayName', sprintf('Orbit %d - %.2f days', orbit_count, days_elapsed));
+ end
+end
+
+% Print the total time spent in orbit adjustment
+fprintf('Total time spent in aerobraking: %.2f seconds, which is approximately %.2f days.\n', total_time, total_time / 86400);
+fprintf('Total number of orbits to achieve target eccentricity: %d\n', orbit_count);
+
+% Add legend and finalize plot
+legend show;
+hold off;
+
+% 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
+
+% Calculate initial semi-major axis
+a = calculate_semi_major_axis(rp, e);
+
+%% Specific Apoapsis Radius
+
+% Initial calculations for apoapsis radius and check condition
+ra_initial = a * (1 + e) - Mars_radius;
+if ra_initial < Apo_Radius
+ fprintf('Initial orbit already has an apoapsis less than 400 km above Mars.\n');
+ return;
+end
+
+% Simulation loop to decrease apoapsis
+orbit_count = 0;
+ra = ra_initial;
+while ra > Apo_Radius
+ orbit_count = orbit_count + 1;
+
+ % Calculate velocity at periapsis
+ vp = sqrt(mu * (2/rp - 1/a));
+
+ % Calculate atmospheric density at periapsis
+ h = rp - Mars_radius; % Altitude above Mars' surface
+ rho = rho0 * exp(-(h - h0) / H);
+
+ % Drag force and velocity change
+ Fd = -0.5 * Cd * rho * vp^2 * A;
+ delta_v = Fd / m; % Change in velocity
+
+ % Update velocity and orbital elements
+ vp_new = vp + delta_v;
+ a = mu / (2 * mu / rp - vp_new^2); % New semi-major axis
+ e = 1 - rp / a; % New eccentricity
+
+ % Calculate new apoapsis radius
+ ra = a * (1 + e) - Mars_radius;
+end
+
+% Output results when apoapsis goes below 400 km
+fprintf('The radius of apoapsis became less than 400 km above Mars on orbit %d.\n', orbit_count);
+fprintf('Apoapsis Radius: %.4f km above Mars\n', ra);
+fprintf('Eccentricity at this orbit: %.4f\n', e);
+
+% Constants
+mu = 42828; % Gravitational parameter for Mars in km^3/s^2
+Mars_radius = 3396; % Radius of Mars in km
+
+% Orbit data
+Rp = Mars_radius + Altitude; % Radius of periapsis in km
+Ra = Mars_radius + Apo_Radius; % Radius of apoapsis in km
+
+% Semi-major axes
+a_initial = (Rp + Ra) / 2; % Semi-major axis of the elliptical orbit
+a_target = Ra; % Semi-major axis of the target circular orbit
+
+% Velocities
+vp_initial = sqrt(mu * (2 / Rp - 1 / a_initial)); % Velocity at periapsis in initial orbit
+va_initial = sqrt(mu * (2 / Ra - 1 / a_initial)); % Velocity at apoapsis in initial orbit
+vc_target = sqrt(mu / a_target); % Velocity in the target circular orbit
+
+% Delta-v calculations
+delta_v_circularize = vc_target - va_initial; % Delta-v required to circularize at apoapsis
+
+% Output results
+fprintf('Initial velocity at periapsis (vp): %.5f km/s\n', vp_initial);
+fprintf('Velocity at apoapsis before circularization (va): %.5f km/s\n', va_initial);
+fprintf('Target circular orbit velocity (vc_target): %.5f km/s\n', vc_target);
+fprintf('Delta-v required at apoapsis to circularize: %.5f km/s\n', delta_v_circularize);
+
+% Constants and Initial Conditions
+Mars_radius = 3396; % Radius of Mars in km
+rp = Mars_radius + Apo_Radius; % Periapsis distance from center of Mars, in km
+e = 0; % Eccentricity
+mu = 42828; % Mars' gravitational parameter in km^3/s^2
+
+% Orbital Elements Calculation
+a = calculate_semi_major_axis(rp, e); % Calculate semi-major axis
+
+% Inclination, RAAN, and Argument of Periapsis (converted to radians)
+i = deg2rad(167.90); % Inclination in radians
+Omega = deg2rad(276.56); % Right Ascension of the Ascending Node (RAAN) in radians
+omega = deg2rad(118.94 + 270); % Argument of periapsis in radians
+
+% Calculate the orbital period using Kepler's third law
+P = 2 * pi * sqrt(a^3 / mu); % Orbital period in seconds
+
+% Time vector spanning one orbital period
+t_vec = linspace(0, P, 1000);
+
+% True anomaly vector from 0 to 2*pi
+theta_vec = linspace(0, 2 * pi, length(t_vec));
+
+% Initialize matrix to store state vectors
+X_matrix = zeros(6, length(t_vec));
+
+% Calculate state vectors across the time vector using the updated COE2RV
+for k = 1:length(t_vec)
+ [RIJK, VIJK] = COE2RV_Elliptic([a, e, i, Omega, omega, theta_vec(k)], mu);
+ X_matrix(:, k) = [RIJK; VIJK]; % Storing position and velocity vectors
+end
+
+%% Plotting
+figure;
+hold on;
+
+% Plot Mars as a sphere with radius 6378 km
+[Xe, Ye, Ze] = sphere(50); % Create a 50x50 sphere
+MarsRadius = 3389.5; % Mars's radius in km
+surf(Xe * MarsRadius, Ye * MarsRadius, Ze * MarsRadius, 'EdgeColor', 'none', 'FaceColor', 'red', 'DisplayName', 'Mars');
+
+%% Plot Satellite's trajectory if X_matrix is defined
+plot3(X_matrix(1, :), X_matrix(2, :), X_matrix(3, :), 'g', 'LineWidth', 2, 'DisplayName', 'Satellite Trajectory');
+
+% Add legends and labels
+legend show;
+xlabel('X (km)');
+ylabel('Y (km)');
+zlabel('Z (km)');
+title('Circular Orbit Around Mars');
+grid on;
+axis equal;
+xlim([-1e4, 1e4]); % Set limits for the x-axis
+ylim([-1e4, 1e4]); % Set limits for the y-axis
+zlim([-1e4, 1e4]); % Set limits for the z-axis
+view(0, 0); % Set the view to 3D perspective
+
+% Legend to identify plot elements
+legend('show', 'Location', 'bestoutside');
+
+% Hold off to finalize the plot
+hold off;
\ No newline at end of file
--
GitLab