diff --git a/Mars_Capture.m b/Mars_Capture.m new file mode 100644 index 0000000000000000000000000000000000000000..6bb5bc9874fc134eb818b4c2099a3e84654455ca --- /dev/null +++ b/Mars_Capture.m @@ -0,0 +1,362 @@ +%% Function Description (Mars Capture Orbit using Patch Conic Approximation) +% This script computes the capture orbit of a spacecraft around Mars, focusing on the Mars +% section of the patch conics approximation transfer. The patch conic approximation divides +% the interplanetary transfer into segments, each governed by the gravitational influence of +% a single celestial body. In this case, the script handles the transition from a hyperbolic +% approach trajectory to an elliptical capture orbit around Mars. + +% The script begins by clearing the workspace, closing figures, and setting the output format +% to ensure a fresh environment. It then defines several constants needed for the simulation, +% including conversion factors, gravitational parameters, and the initial conditions for the +% spacecraft's orbit upon entering Mars' Sphere of Influence (SOI). + +% Key steps and calculations in this script are: +% - Defining the initial and final conditions for the spacecraft's transfer from Earth to Mars, +% including launch and arrival dates, positions, and velocities. The script converts these dates +% to Ephemeris Time using SPICE functions to ensure accurate orbital calculations. +% - Loading SPICE kernels to access planetary ephemerides and gravitational models, +% providing accurate state vectors for Mars and other celestial bodies. +% - Using a Lambert solver to determine the initial and final velocity vectors required for +% the interplanetary transfer. The Lambert solver calculates the optimal trajectory using +% the patched conic approximation, which simplifies the trajectory into separate segments. +% - Computing the spacecraft's hyperbolic trajectory as it approaches Mars and transitions into +% a capture orbit. The script calculates orbital elements such as inclination, right ascension +% of the ascending node (RAAN), argument of periapsis, and eccentricity to characterize the +% capture orbit. +% - Iterating over a range of times of flight (TOFs) and transfer modes to find the optimal +% configuration that minimizes the delta-V required for the Mars capture maneuver. The script +% evaluates each configuration and selects the one that results in the lowest delta-V, +% indicating the most efficient trajectory. +% - Propagating the trajectory using numerical integration (ode45) to simulate the spacecraft's +% path around Mars. The output state vectors are used to plot the spacecraft's trajectory and +% visualize the transition from hyperbolic to elliptical orbit. +% - Plotting the results, including Mars' position, the spacecraft's trajectory, and other relevant +% elements, to provide a visual representation of the capture orbit around Mars. + +% The script is designed for flexibility, allowing for adjustments to various parameters, such as +% initial orbital elements, observation periods, and transfer modes. It provides a comprehensive +% approach to computing capture orbits for interplanetary missions using the patch conic approximation. + +% 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; + +%% Constants +% Define the conversion factor from kilometers to Astronomical Units (AU) +kmToAU = 1 / 149597870.7; % Inverse of the average distance from Mars to the Sun in km +% Define the conversion factor from degrees to radians +deg2Rad = pi / 180; % Pi radians equals 180 degrees + +% Semi-major axes of Mars and Mars orbits in km (average distances) +a_Mars = 149597870.7; % 1 AU in km, average distance from the Mars to the Sun +a_mars = 227939200; % Approximately 1.524 AU in km, average distance from Mars to the Sun + +mu_mars = 4.28284e4; % Gravitational parameter of the sun, needed for orbital mechanics calculations +%% Calling Spice Functions +% Add paths to required SPICE toolboxes for planetary data access +addpath(genpath('C:/Users/TomWe/OneDrive/Desktop/University [MSc]/Major Project/MATLAB Code/SPICE/mice')); +addpath(genpath('C:/Users/TomWe/OneDrive/Desktop/University [MSc]/Major Project/MATLAB Code/SPICE/generic_kernels')); +addpath(genpath('C:/Users/TomWe/OneDrive/Desktop/University [MSc]/Major Project/MATLAB Code/Functions')); + +% Load SPICE kernels for planetary and lunar ephemerides and gravitational models +cspice_furnsh(which('de430.bsp')); % Binary SPK file containing planetary ephemeris +cspice_furnsh(which('naif0012.tls')); % Leapseconds kernel, for time conversion +cspice_furnsh(which('gm_de431.tpc')); % Text PCK file containing gravitational constants +cspice_furnsh(which('pck00010.tpc')); % Planetary constants kernel (sizes, shapes, etc.) + +%% Launch Details and Configuration +LaunchDate = '2 Sep 2027 12:30:00 (UTC)'; % Define the launch date in UTC +ArrivalDate = '5 Sep 2027 00:40:00 (UTC)'; % Define the arrival date in UTC + +%% Data Collection +% Convert both dates to Ephemeris Time using SPICE functions +et0 = cspice_str2et(LaunchDate); % Convert the launch date to seconds past J2000 +et_end = cspice_str2et(ArrivalDate); % Convert the arrival date to seconds past J2000 + +% Calculate the total duration in seconds between start and end dates +total_duration = et_end - et0; % Compute the time difference in seconds + +% Define a finer resolution for the time vector +time_step = 3600; % Define time step in seconds (e.g., 60 seconds for 1-minute intervals) + +% Generate a time vector from start to end date using the finer resolution +et = et0:time_step:et_end; % Generate evenly spaced time points including both endpoints + +% Retrieve state vectors for each time point +num_points = length(et); % Calculate the number of points in the time vector +Xmars = zeros(6, num_points); % Preallocate array for Mars state vectors +XMars = zeros(6, num_points); % Preallocate array for Mars state vectors + +for tt = 1:num_points % Loop over each time point + Xmars(:, tt) = cspice_spkezr('4', et(tt), 'ECLIPJ2000', 'NONE', '10'); % Retrieve Mars' state vector at time tt + XMars(:, tt) = cspice_spkezr('399', et(tt), 'ECLIPJ2000', 'NONE', '10'); % Retrieve Mars's state vector at time tt +end + +radius = 150 + 3396; +r = [521709.873992; -223450.000055; -105659.613183]; +v = [-2.329180; 0.997457; 0.471722]; + +mu = 42828; % Gravitational parameter of Mars (km^3/s^2) + +% Calculate orbital elements +[inclination, RAAN, argument_of_periapsis, eccentricity] = geometry_calculations(r, v, mu); + +% Display the calculated values +fprintf('Inclination: %.2f degrees\n', inclination); +fprintf('RAAN: %.2f degrees\n', RAAN); +fprintf('Argument of Periapsis: %.2f degrees\n', argument_of_periapsis); +fprintf('Eccentricity: %.4f\n', eccentricity); + +% Constants and Initial Conditions +Mars_radius = 3396; % Radius of Mars in km +rp = Mars_radius + 100; % Periapsis distance from center of Mars, in km +e = 0.9; % 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 + +% Display the semi-major axis +fprintf('The semi-major axis for periapsis %d km and eccentricity %.2f is %.2f km.\n', rp, e, a); + +% 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 + +% Spacecraft's initial position vector (example values) +r0 = [521709.873992; -223450.000055; -105659.613183]; + +rf = X_matrix(1:3, 300); + + +%% Launch and Arrival Details Configuration +ArrivalDate = '5 Sep 2027 00:40:00 (UTC)'; % Fixed arrival date in UTC + +% Convert the fixed arrival date to Ephemeris Time using SPICE functions +et_end = cspice_str2et(ArrivalDate); % Convert the arrival date to seconds past J2000 + +%% Simulation Configuration +DM = 1; % Transfer mode setting (+1 for Type I trajectory) + +% Time of Flight (TOF) range and step size +min_TOF_days = 2; % Minimum TOF in days +max_TOF_days = 3; % Maximum TOF in days +time_step_minutes = 0.1; % Increment step in minutes +time_step = time_step_minutes * 60; % Convert minutes to seconds + +% Generate TOFs from 1 day to 4 days in 20 minute increments +TOFs = min_TOF_days*86400:time_step:max_TOF_days*86400; + +% Preallocate arrays for initial and final velocities +v0s = zeros(3, length(TOFs)); +vfs = zeros(3, length(TOFs)); + +%% Lambert Solver Simulation Loop +for i = 1:length(TOFs) + % Calculate the corresponding launch date for each TOF + et0 = et_end - TOFs(i); + LaunchDate = cspice_et2utc(et0, 'C', 0); % Convert back to readable date, if needed for display + + % Compute the transfer time in seconds (delta_t is TOF) + delta_t = TOFs(i); + + try + % Call the Lambert solver function with the current TOF + [v0, vf] = lambert_solver(r0, rf, delta_t, mu, DM); + + % Store results in preallocated arrays + v0s(:, i) = v0; + vfs(:, i) = vf; + catch ME + % Handle errors by logging or setting to NaN + fprintf('Failed to solve Lambert problem for TOF %d seconds: %s\n', delta_t, ME.message); + v0s(:, i) = NaN(3,1); + vfs(:, i) = NaN(3,1); + end +end + +% Given target initial velocity vector +v_enter_value = [-2.329180; 0.997457; 0.471722]; + +% Assuming v0s is the matrix containing all calculated v0 vectors for different TOFs +% Each column in v0s corresponds to a different TOF +% TOFs is the array containing all the TOFs used in the simulation + +% Calculate the Euclidean distances between each calculated v0 and the given v_enter_value +distances = sqrt(sum((v0s - v_enter_value).^2, 1)); + +% Find the index of the minimum distance +[min_distance, min_index] = min(distances); + +% Retrieve the TOF corresponding to the minimum distance +closest_TOF = TOFs(min_index); +closest_v0 = v0s(:, min_index); + +% Fixed Arrival Date +ArrivalDate = '5 Sep 2027 00:40:00 (UTC)'; + +% Convert the arrival date to Ephemeris Time using SPICE functions +et_end = cspice_str2et(ArrivalDate); + +% Calculate the departure date's Ephemeris Time +et_departure = et_end - closest_TOF; + +% Convert the departure date's Ephemeris Time back to a UTC string +DepartureDate = cspice_et2utc(et_departure, 'C', 0); + +% Display results +fprintf('The closest TOF is %d seconds, corresponding to %d days and %d hours.\n', closest_TOF, floor(closest_TOF/86400), mod(floor(closest_TOF/3600), 24)); +fprintf('The closest initial velocity v0 is: [%f; %f; %f]\n', closest_v0); +fprintf('The distance to the target velocity is %f.\n', min_distance); +% Display the calculated departure date +fprintf('The departure date corresponding to the closest TOF is: %s\n', DepartureDate); + +% Define TOF, r0, mu, and DM_values +TOF = closest_TOF; % Closest TOF previously determined +DM_values = [1, -1]; % DM values to loop over + +% Preallocate structures to store results +min_delta_v = inf; % Initialize with a large number +best_result = struct('DM', 0, 'VectorNumber', 0, 'DeltaV', inf, 'rf', zeros(3,1), 'vf', zeros(3,1)); + +% Loop over the two DM values +for dm_index = 1:length(DM_values) + DM = DM_values(dm_index); % Current DM value + + % Loop through each column of X_matrix (up to the 400th column) + for i = 1:400 + % Extract the current rf and corresponding velocity vector from the matrix + rf = X_matrix(1:3, i); + velocity_vector = X_matrix(4:6, i); + + try + % Call the Lambert solver function + [~, vf] = lambert_solver(r0, rf, TOF, mu, DM); + + % Calculate the delta v between the computed vf and the velocity vector from X_matrix + delta_v = norm(vf - velocity_vector); + + % Check if this delta v is the smallest found so far + if delta_v < min_delta_v + % Update minimum delta v and store the corresponding data + min_delta_v = delta_v; + best_result.DM = DM; + best_result.VectorNumber = i; + best_result.DeltaV = delta_v; + best_result.rf = rf; + best_result.vf = vf; + end + catch ME + % Log errors if the Lambert solver fails (keep this simple) + fprintf('Error for rf vector %d with DM %d: %s\n', i, DM, ME.message); + end + end +end + +% Display the best results +fprintf('Best Match Details:\n'); +fprintf('DM Direction: %d\n', best_result.DM); +fprintf('Vector Number: %d\n', best_result.VectorNumber); +fprintf('Delta V: %f\n', best_result.DeltaV); +fprintf('Position Vector rf: [%f; %f; %f]\n', best_result.rf); +fprintf('Velocity Vector vf: [%f; %f; %f]\n', best_result.vf); + +DM = best_result.DM; +idx = best_result.VectorNumber; +rf = X_matrix(1:3, idx); +delta_t = closest_TOF; + +% Call the Lambert solver function with the gathered data +[v0, vf] = lambert_solver(r0, rf, delta_t, mu, DM); % Compute initial and final velocity vectors using the Lambert solver + +v_enter_value = [-2.329180; 0.997457; 0.471722]; % Given value +v_orbit_value = X_matrix(4:6, idx); + +% Delta-v Calculation for Initial Burn (Entering Mars SOI to v0) +delta_v0 = norm(v_enter_value - v0); +delta_vf = norm(v_orbit_value - vf); +delta_v_total = delta_v0 + delta_vf; + +% Display the results +fprintf('Delta-v for Initial Burn (Entering Mars SOI to v0): %.6f km/s\n', delta_v0); +fprintf('Delta-v for Final Burn (Mars Orbit Insertion): %.6f km/s\n', delta_vf); +fprintf('Total Delta-v for the Mars Orbit Insertion: %.6f km/s\n', delta_v_total); + +% Initial state vector (position and velocity) +initial_state = [r0; v0]; % Combine the position vector and initial velocity vector into a single state vector + +%% Propagate Orbit +% Define the time span of the simulation with more points +num_points = 100000; % Define the number of points for finer resolution +tspan = linspace(0, delta_t, num_points); % Generate a finely spaced time vector for the simulation + +% Propagate the trajectory using ode45 +options = odeset('RelTol', 1e-13, 'AbsTol', 1e-13); % Set options for numerical solver to high precision +[t_out, state_out] = ode45(@(t, y) Orbit_Propagation(t, y, mu), tspan, initial_state, options); % Numerically integrate the trajectory + +% Extract position and velocity vectors +positions = state_out(:, 1:3); % Position data +velocities = state_out(:, 4:6); % Velocity data + + +%% 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, :), 'm', 'LineWidth', 0.5, 'DisplayName', 'Satellite Trajectory'); + +% Plotting the transfer trajectory +plot3(positions(:, 1), positions(:, 2), positions(:, 3), 'm-', 'LineWidth', 1, 'DisplayName', 'Transfer Trajectory'); + +% Add legends and labels +legend show; +xlabel('X (km)'); +ylabel('Y (km)'); +zlabel('Z (km)'); +title('Mars Capture Orbit (Hyperbolic to Elliptical)'); +grid on; +axis equal; +xlim([-1e5, 1e5]); % Set limits for the x-axis +ylim([-1e5, 1e5]); % Set limits for the y-axis +zlim([-1e5, 1e5]); % Set limits for the z-axis +view(3); % 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