diff --git a/Total_Delta_V_Lunar_Assist.m b/Total_Delta_V_Lunar_Assist.m
deleted file mode 100644
index 3292e676cb754f6bbb5a4efbd0780914f4e886d9..0000000000000000000000000000000000000000
--- a/Total_Delta_V_Lunar_Assist.m
+++ /dev/null
@@ -1,602 +0,0 @@
-%% Function Description (Trajectory Computation between Earth and Mars with Lunar Gravity Assist)
-% This script computes the trajectory of a spacecraft traveling from Earth to Mars using
-% a lunar gravity assist. The lunar gravity assist is a maneuver that uses the Moon's gravitational
-% field to alter the spacecraft's trajectory and velocity, enabling a more efficient transfer to Mars.
-% This approach reduces the total delta-V required for the mission, saving fuel and optimizing the transfer path.
-
-% The script begins by clearing the workspace, closing figures, and setting the output format
-% for a fresh environment. It also sets text size and formatting options for consistency in plots
-% and visualizations. It then defines several constants needed for the simulation, including conversion
-% factors, gravitational parameters, and the initial orbital elements for the spacecraft's Earth parking orbit.
-
-% Key steps and calculations in this script are:
-% - Defining the initial conditions for the spacecraft's orbit around Earth, including the semi-major axis,
-% eccentricity, inclination, right ascension of the ascending node (RAAN), and argument of periapsis.
-% - Loading SPICE kernels to access planetary and lunar ephemerides and gravitational models,
-% providing accurate state vectors for the Earth, Moon, and other celestial bodies involved in the transfer.
-% - Calculating state vectors (position and velocity) for Earth and the Moon over a specified observation period
-% using SPICE functions. These vectors are crucial for determining the spacecraft's trajectory relative to
-% Earth, the Moon, and the Sun.
-% - Utilizing a Lambert solver to determine the initial and final velocity vectors required for the interplanetary
-% transfer from Earth to the Moon, and subsequently from the Moon to Mars. The Lambert solver accounts for
-% gravitational influences and computes the optimal transfer trajectory using the patched conic approximation.
-% - Implementing a lunar gravity assist by computing the change in the spacecraft's velocity vector (v-infinity)
-% due to the Moon's gravitational field. This change is used to adjust the spacecraft's trajectory towards Mars.
-% - Iterating over a range of possible pump (alpha) and crank (kappa) angles to find the optimal configuration
-% that minimizes the delta-V required for the transfer. 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 through
-% Earth's Sphere of Influence (SOI), its encounter with the Moon, and its eventual exit from Earth's SOI.
-% - Plotting the results, including the Earth's position, the Moon's trajectory, the spacecraft's path, and the
-% V-infinity globe showing the delta-V contours for different lunar gravity assist configurations.
-
-% The script is designed for flexibility, allowing for adjustments to various parameters, such as the initial
-% orbital elements, observation periods, and transfer modes. It provides a comprehensive approach to computing
-% efficient transfer trajectories for interplanetary missions using gravity assists.
-
-% 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;
-% Set the text size format
-set(groot, 'DefaultAxesFontSize', 18, 'DefaultTextFontSize', 18, 'DefaultAxesFontName', 'Times New Roman', 'DefaultTextFontName', 'Times New Roman');
-
-%% 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.)
-
-% Semi-major axes of Earth and moon orbits in km (average distances)
-a_earth = 149597870.7; % 1 AU in km, average distance from the Earth to the Sun
-
-a = 6878; % Semi-major axis in km
-e = 0.0; % Eccentricity
-i = deg2rad(75); % Inclination in radians
-Omega = deg2rad(270 + 309.1746);
-
-omega = 0; % Argument of periapsis in radians
-meanAnomalyDeg = 0; % degrees
-mu = 3.986e5; % Earth's gravitational parameter in km^3/s^2
-
-%% Launch Details and Configuration
-LaunchDate = '24 Oct 2026 16:32:00 (UTC)'; % Update this line with your desired start date 17, 3.77096302004596 24
-ArrivalDate = '26 Oct 2026 13:16:00 (UTC)'; % Update this line with your desired end date
-
-Step = 0.1;
-
-DM = 1; % Transfer mode setting:
-% -1 for a Type II trajectory, which takes the longer path, more than 180 degrees around the central body.
-% This path may be chosen for specific mission requirements like thermal management or avoiding solar conjunctions.
-
-% +1 for a Type I trajectory, representing the shorter path, less than 180 degrees around the central body.
-% This is generally more energy-efficient and quicker, suitable for missions prioritizing speed and fuel efficiency.
-
-%% 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; % Difference in ephemeris time gives duration in seconds
-
-% Calculate the number of days between the start and end date
-num_days = total_duration / 60; % Convert seconds to days using the number of seconds in a day
-
-% Round num_days to the nearest whole number to avoid fractional days
-num_days = round(num_days); % Rounding to ensure an integer number of days for vector generation
-
-% Generate a time vector from start to end date using the calculated number of days
-et = linspace(et0, et_end, num_days + 1); % Generates evenly spaced time points including both endpoints
-
-for tt = 1:length(et) % Loop over each time point
- Xmoon(:, tt) = cspice_spkezr('301', et(tt), 'ECLIPJ2000', 'NONE', '399');
-end
-
-% Constants
-mu_Earth = 398600.4418; % Gravitational parameter of Earth in km^3/s^2
-
-% Define gravitational parameter for Moon (adjust as needed)
-mu_moon = 4902.8; % Moon's gravitational parameter in km^3/s^2
-
-% Define the critical radius of periapsis
-critical_rp = 1737.4 + 100; % Radius that indicates a potential crash
-
-% Calculate the orbital period (s) of the satellite using Kepler's third law
-P = 2 * pi * sqrt(a^3 / mu);
-
-% Create a time vector 't_vec' with 1000 points spanning one orbital period
-t_vec = linspace(0, P, 100000);
-
-% Calculate the Mean motion (rad/s)
-n = sqrt(mu / a^3);
-
-% Define time vector and orbital elements
-theta_vec = linspace(0, 2*pi, length(t_vec)); % True anomaly from 0 to 2*pi
-
-% Initialize matrix to store state vectors
-X_matrix = zeros(6, length(t_vec));
-
-% Calculate state vectors across the time vector
-for k = 1:length(t_vec)
- X_matrix(:, k) = COE2RV([a, e, i, Omega, omega, theta_vec(k)], mu);
-end
-
-% Extract initial and final state vectors
-X_initial = X_matrix(:, 1);
-
-% Plot Satellite's trajectory if X_matrix is defined
-if exist('X_matrix', 'var')
-
- % Assuming a circular orbit for simplicity where mean anomaly = true anomaly
- meanAnomalyRad = deg2rad(meanAnomalyDeg); % Convert degrees to radians
-
- % Find the index where the true anomaly (theta_vec) is closest to 90 degrees
- [~, idx90] = min(abs(theta_vec - meanAnomalyRad));
-
-end
-
-%% Lambert Solver
-delta_t = (et_end - et0); % Calculate the transfer time in seconds
-epsilon = 1e-6; % Set a small number for numerical precision in calculations
-
-% Spacecraft initial position vector
-r0 = X_matrix(1:3, 1); % Extract Earth's position vector at launch
-
-% Spacecraft final position vector at Moon
-rf = Xmoon(1:3, end); % Extract Moon' position vector at the arrival time
-
-% Call the Lambert solver function with the gathered data
-[v0, vf] = lambert_solver_moon(r0, rf, delta_t, mu, DM); % Compute initial and final velocity vectors using the Lambert solver
-
-% Vectors
-r_enc_vector = Xmoon(1:3, end);
-v_sc_vector = vf;
-v_ga_vector = Xmoon(4:6, end);
-
-% Assume v_sc is the same as the escape velocity at this point
-v_inf_vector = vf - v_ga_vector;
-
-% Calculate magnitudes of vectors
-r_enc = norm(r_enc_vector);
-v_inf = norm(v_inf_vector);
-v_ga = norm(v_ga_vector);
-v_sc = norm(v_sc_vector);
-
-% Display the results
-fprintf('V-Infinity Value = %.2f Km/s.\n', v_inf);
-
-% Calculate vector q3 as the normalized version of v_ga
-q2 = v_ga_vector / norm(v_ga_vector);
-
-% Calculate vector q2 as the normalized cross product of v_ga and r_enc
-q3 = cross(r_enc_vector, v_ga_vector);
-q3 = q3 / norm(q3); % Normalize to make it a unit vector
-
-% Calculate vector q1 as the cross product of q2 and q3
-q1 = cross(q2, q3); % This is already unit length if q2 and q3 are unit vectors
-
-% Initialize a figure for plotting
-figure;
-hold on; % Hold the plot for multiple data series
-xlabel('Alpha (pump angle) [degrees]');
-ylabel('Kappa (crank angle) [degrees]');
-title('V-Infinity Globe Delta-V Contour Plot for Lunar Gravity Assist');
-grid on;
-
-% Define the range and resolution of alpha and kappa
-alpha_range = -180:Step:180;
-kappa_range = -90:Step:90;
-[Alpha, Kappa] = meshgrid(alpha_range, kappa_range);
-DeltaVGrid = zeros(size(Alpha)); % Initialize delta-V grid to zeros
-
-% Define v_Moon_required
-v_Moon_required = [-1.973070; 2.268442; 0.492901];
-
-% Loop over alpha and kappa angles
-for i = 1:length(alpha_range)
- for j = 1:length(kappa_range)
- alpha_deg = alpha_range(i);
- alpha = deg2rad(alpha_deg);
- kappa_deg = kappa_range(j);
- kappa = deg2rad(kappa_deg);
-
- % Calculate new v_infinity for current alpha and kappa
- v_infinity_new = v_inf * sin(alpha) * cos(kappa) * q1 + ...
- v_inf * cos(alpha) * q2 + ...
- v_inf * sin(alpha) * sin(kappa) * q3;
-
- % Calculate geocentric velocity
- v_geocentric = v_infinity_new + v_ga_vector;
- v_magnitude = norm(v_geocentric);
-
- % Calculate delta-V
- DeltaVGrid(j, i) = norm(v_geocentric - v_Moon_required);
-
- % Specific orbital energy
- epsilon = (v_magnitude^2) / 2 - mu_Earth / norm(r_enc_vector); % Calculate specific orbital energy
-
- % Turning angle and eccentricity
- dot_product = dot(v_infinity_new, v_inf_vector);
- angle_magnitude = norm(v_infinity_new) * norm(v_inf_vector);
- turning_angle = acos(dot_product / angle_magnitude);
- eccentricity = -1 / cos(turning_angle);
-
- rp = mu_moon / (norm(v_infinity_new)^2 * (eccentricity - 1));
-
- % Check for crashing into the Moon or not escaping Earth's SOI
- if rp < critical_rp || epsilon < 0
- DeltaVGrid(j, i) = NaN; % Set to NaN to avoid plotting
- else
- % Calculate delta-V
- DeltaVGrid(j, i) = norm(v_geocentric - v_Moon_required);
- end
- end
-end
-
-% Define the delta-V range and step size
-deltaV_min = 0; % Minimum delta-V in km/s
-deltaV_max = 4; % Maximum delta-V in km/s
-deltaV_step = 0.1; % Step size for each color change in km/s
-
-% Calculate the number of levels in the colormap based on the range and step size
-nLevels = floor((deltaV_max - deltaV_min) / deltaV_step);
-
-% Find the minimum delta-V in the grid excluding NaN values
-[minDeltaV, idx] = min(DeltaVGrid(:));
-[iRow, iCol] = ind2sub(size(DeltaVGrid), idx); % Convert linear index to subscript
-
-% Extract the corresponding alpha and kappa values
-minAlpha = alpha_range(iCol); % Corresponding alpha value
-minKappa = kappa_range(iRow); % Corresponding kappa value
-
-% Display the results
-fprintf('The minimum delta-V is %.5f km/s, found at alpha (Pump) = %.2f degrees and kappa (Crank) = %.2f degrees.\n', minDeltaV, minAlpha, minKappa);
-
-% Plot the contour map
-contourf(Alpha, Kappa, DeltaVGrid, linspace(deltaV_min, deltaV_max, nLevels));
-colormap(jet(nLevels)); % Color map setup
-hColorbar = colorbar; % Display a color bar indicating the delta-V values
-caxis([deltaV_min deltaV_max]); % Set color axis scaling
-
-% Define the tick marks for the color bar
-ticks = deltaV_min:0.2:deltaV_max;
-set(hColorbar, 'Ticks', ticks, 'TickLabels', arrayfun(@num2str, ticks, 'UniformOutput', false));
-
-% Optional: Customize tick label formatting if necessary
-set(hColorbar, 'TickLabelInterpreter', 'latex', 'FontSize', 10);
-
-% Mark the minimum delta-V point with a hollow magenta circle
-hold on;
-plot(minAlpha, minKappa, 'mx', 'MarkerSize', 10, 'MarkerFaceColor', 'none', 'MarkerEdgeColor', 'm', 'LineWidth', 2);
-
-% Set axes labels and title
-xlabel('Pump Angle [degrees]');
-ylabel('Crank Angle [degrees]');
-title('V-Infinity Globe Delta-V Contour Plot for Lunar Gravity Assist');
-
-% Add legend entry for the minimum delta-V point with more precise formatting
-legend('Minimum Delta-V', sprintf('Min Δv: %.6f km/s at α (Pump Angle) = %.2f°, κ (Crank Angle) = %.2f°', minDeltaV, minAlpha, minKappa));
-
-% Adjust the legend position and box
-hLegend = legend('show');
-set(hLegend, 'Location', 'bestoutside', 'Box', 'on');
-
-hold off;
-
-% Legend and colormap settings
-legend('show');
-
-v0_delta = norm(v0 - X_matrix(4:6, idx90));
-v0
-X_matrix(4:6, idx90)
-v0 - X_matrix(4:6, idx90)
-% Display the delta-v
-disp('Delta-v for the inital transfer (km/s):');
-disp(v0_delta); % Output the computed delta-v value
-
-delta_v_final = 1.109200;
-
-% Display the delta-v
-disp('Delta-v for the final transfer (km/s):');
-disp(delta_v_final); % Output the computed delta-v value
-
-% Display the delta-v
-disp('Total Delta-v for the transfer (km/s):');
-disp(delta_v_final + v0_delta); % Output the computed delta-v value
-
-Pump = deg2rad(minAlpha);
-Crank = deg2rad(minKappa);
-
-% Calculate new v_infinity for current alpha and kappa
-v_infinity_new = v_inf * sin(Pump) * cos(Crank) * q1 + v_inf * cos(Pump) * q2 + v_inf * sin(Pump) * sin(Crank) * q3;
-
-v_geocentric = v_infinity_new + v_ga_vector;
-
-v_magnitude = norm(v_geocentric);
-
-% Specific orbital energy
-epsilon = (v_magnitude^2) / 2 - mu_Earth / norm(r_enc_vector); % Calculate specific orbital energy
-
-% Turning angle and eccentricity
-dot_product = dot(v_infinity_new, v_inf_vector);
-angle_magnitude = norm(v_infinity_new) * norm(v_inf_vector);
-turning_angle = acos(dot_product / angle_magnitude);
-turning_angle_deg = rad2deg(turning_angle);
-eccentricity = -1 / cos(turning_angle);
-
-rp = mu_moon / (norm(v_infinity_new)^2 * (eccentricity - 1));
-
-Velocity_Increase = norm(v_geocentric) - norm(vf);
-
-% Display the results
-% Display vector information along with magnitudes
-fprintf('Velocity Vector of Spacecraft Before Lunar Assist (Geocentric) (km/s): [%f, %f, %f], Magnitude: %f km/s\n', vf, norm(vf));
-fprintf('Velocity Vector of Spacecraft After Lunar Assist (Geocentric) (km/s): [%f, %f, %f], Magnitude: %f km/s\n', v_geocentric, norm(v_geocentric));
-fprintf('Velocity Increase in Geocentric Frame from Lunar Gravity Assist: %f km/s\n', Velocity_Increase); % Assuming vf is before and v_geocentric is after for increase calculation
-fprintf('Required Moon Velocity (v_Moon_required) (km/s): [%f, %f, %f], Magnitude: %f km/s\n', v_Moon_required, norm(v_Moon_required));
-fprintf('Velocity Vector of Spacecraft Before Lunar Assist (Selenocentric) (km/s): [%f, %f, %f], Magnitude: %f km/s\n', v_inf_vector, norm(v_inf_vector));
-fprintf('Velocity Vector of Spacecraft After Lunar Assist (Selenocentric) (km/s): [%f, %f, %f], Magnitude: %f km/s\n', v_infinity_new, norm(v_infinity_new));
-
-fprintf('Specific Orbital Energy (epsilon): %f km^2/s^2\n', epsilon);
-fprintf('Turning Angle: %f degrees\n', turning_angle_deg);
-fprintf('Eccentricity: %f\n', eccentricity);
-fprintf('Periapsis Distance (rp): %f km\n', rp);
-
-% Define the time span of the simulation
-tspan_1 = [0 delta_t]; % Time span from t = 0 to the transfer time in seconds for the simulation
-
-% Initial state vector (position and velocity)
-initial_state_1 = [r0; v0]; % Combine the position vector and initial velocity vector into a single state vector
-
-%% Propagate Orbit
-% Propagate the trajectory using ode45
-options = odeset('RelTol', 1e-13, 'AbsTol', 1e-13); % Set options for numerical solver to high precision
-[t_out_1, state_out_1] = ode45(@(t, y) Orbit_Propagation(t, y, mu), tspan_1, initial_state_1, options); % Numerically integrate the trajectory
-
-% Extract the position vectors from the output state for plotting
-positions_1 = state_out_1(:, 1:3); % Extract position data from the output state for subsequent use in plotting
-velocities_1 = state_out_1(:, 4:6); % Extract velocity data from the output state for subsequent use in plotting
-
-%% Propagate Orbit
-TOF_Leave_SOI = delta_t * 1.935807;
-
-tspan_2 = [0 TOF_Leave_SOI];
-
-% Initial state vector (position and velocity)
-initial_state_2 = [r_enc_vector; v_geocentric]; % Combine the position vector and initial velocity vector into a single state vector
-
-% Propagate the trajectory using ode45
-[t_out_2, state_out_2] = ode45(@(t, y) Orbit_Propagation(t, y, mu), tspan_2, initial_state_2, options); % Numerically integrate the trajectory
-
-v_leaving_Moon = state_out_2(1, 4:6);
-disp('Velocity vector after Lunar Assist (km/s):');
-disp(v_leaving_Moon); % Output the computed delta-v value
-
-v_leaving_SOI = state_out_2(end, 4:6);
-disp('Velocity vector leaving Earths SOI (km/s):');
-disp(v_leaving_SOI); % Output the computed delta-v value
-
-disp('Velocity vector required to transfer to Mars at SOI (km/s):');
-disp(v_Moon_required'); % Output the computed delta-v value
-
-diff = v_Moon_required' - v_leaving_SOI;
-norm_diff = norm(diff);
-
-disp('Delta V vector (km/s):');
-disp(diff); % Output the computed delta-v value
-disp('Delta V magnitude (km/s):');
-disp(norm_diff); % Output the computed delta-v value
-
-Total_Delta_V = v0_delta + norm_diff;
-disp('Total velocity required for inital transfer (Lunar Assist) (km/s):');
-disp(Total_Delta_V); % Output the computed delta-v value
-
-disp('Total velocity required for inital transfer (No Lunar Assist) (km/s):');
-disp(4.0357); % Output the computed delta-v value
-
-V_Saved = 4.03578041578003 - Total_Delta_V;
-disp('Total velocity saved using a Lunar gravity assist (km/s):');
-disp(V_Saved); % Output the computed delta-v value
-
-%% Propagate Orbit
-% 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_1] = ode45(@(t, y) Orbit_Propagation(t, y, mu), tspan_1, initial_state_1, options); % Numerically integrate the trajectory
-
-% Extract the position vectors from the output state for plotting
-positions_1 = state_out_1(:, 1:3); % Extract position data from the output state for subsequent use in plotting
-velocities_1 = state_out_1(:, 4:6); % Extract velocity data from the output state for subsequent use in plotting
-
-%% Find the index where the distance to Xmoon(1:3, end) is < 66100 km on approach
-Xmoon_end = Xmoon(1:3, end); % Assuming Xmoon is defined with your data
-approach_idx = find(arrayfun(@(idx) norm(positions_1(idx, :) - Xmoon_end') < 66100, 1:size(positions_1, 1)), 1, 'first');
-
-%% Propagate Orbit
-TOF_Leave_SOI = delta_t * 1.935807
-
-tspan_2 = [0 TOF_Leave_SOI];
-
-Exit_SOI = [245873; 273749; 41212.8];
-
-% Initial state vector (position and velocity)
-initial_state_2 = [r_enc_vector; v_geocentric]; % Combine the position vector and initial velocity vector into a single state vector
-
-% Propagate the trajectory using ode45
-[t_out, state_out_2] = ode45(@(t, y) Orbit_Propagation(t, y, mu), tspan_2, initial_state_2, options); % Numerically integrate the trajectory
-
-% Extract the position vectors from the output state for plotting
-positions_2 = state_out_2(:, 1:3); % Extract position data from the output state for subsequent use in plotting
-velocities_2 = state_out_2(:, 4:6); % Extract velocity data from the output state for subsequent use in plotting
-
-%% Find the index where the distance to Xmoon(1:3, end) is > 66100 km on departure
-departure_idx = find(arrayfun(@(idx) norm(positions_2(idx, :) - Xmoon_end') > 66100, 1:size(positions_2, 1)), 1, 'first');
-
-%% Display results
-fprintf('Approach Index: %d, Time: %f seconds\n', approach_idx, t_out_1(approach_idx));
-fprintf('Departure Index: %d, Time: %f seconds\n', departure_idx, t_out_2(departure_idx));
-
-%% Display Geocentric Position Vectors at Critical Indices
-% Approach Index Geocentric Position
-approach_pos_geocentric = positions_1(approach_idx, :);
-fprintf('Approach Position Geocentric (km): [%f, %f, %f]\n', approach_pos_geocentric);
-
-% Departure Index Geocentric Position
-departure_pos_geocentric = positions_2(departure_idx, :);
-fprintf('Departure Position Geocentric (km): [%f, %f, %f]\n', departure_pos_geocentric);
-
-%% Calculate and Display Selenocentric Position Vectors
-% Approach Position Selenocentric
-approach_pos_selenocentric = approach_pos_geocentric - Xmoon_end';
-approach_pos_selenocentric_mag = norm(approach_pos_selenocentric); % Calculate magnitude
-fprintf('Approach Position Selenocentric (km): [%f, %f, %f], Magnitude: %f km\n', approach_pos_selenocentric, approach_pos_selenocentric_mag);
-
-% Departure Position Selenocentric
-departure_pos_selenocentric = departure_pos_geocentric - Xmoon_end';
-departure_pos_selenocentric_mag = norm(departure_pos_selenocentric); % Calculate magnitude
-fprintf('Departure Position Selenocentric (km): [%f, %f, %f], Magnitude: %f km\n', departure_pos_selenocentric, departure_pos_selenocentric_mag);
-
-%% Plotting
-figure;
-hold on;
-
-% Plot Earth as a sphere with radius 6378 km
-[Xe, Ye, Ze] = sphere(50); % Create a 50x50 sphere
-earthRadius = 6378; % Earth's radius in km
-surf(Xe * earthRadius, Ye * earthRadius, Ze * earthRadius, 'EdgeColor', 'none', 'FaceColor', 'blue', 'DisplayName', 'Earth');
-
-%% Plot Satellite's trajectory if X_matrix is defined
-if exist('X_matrix', 'var')
- plot3(X_matrix(1, :), X_matrix(2, :), X_matrix(3, :), 'm--', 'LineWidth', 0.5, 'DisplayName', 'Satellite Parking Orbit Trajectory');
- plot3(X_matrix(1, 1), X_matrix(2, 1), X_matrix(3, 1), 'g.', 'MarkerSize', 10, 'DisplayName', 'Satellite Initial Position');
-end
-
-% Plot Moon trajectory
-plot3(Xmoon(1, :), Xmoon(2, :), Xmoon(3, :), 'Color', [0.5, 0.5, 0.5], 'LineWidth', 2, 'DisplayName', 'Moon Trajectory');
-
-% Annotations for initial and final positions of the Moon
-plot3(Xmoon(1, 1), Xmoon(2, 1), Xmoon(3, 1), 'o', 'MarkerFaceColor', 'none', 'Color', [0.5, 0.5, 0.5], 'DisplayName', 'Moon Initial Position');
-plot3(Xmoon(1, end), Xmoon(2, end), Xmoon(3, end), 'o', 'MarkerFaceColor', [0.5, 0.5, 0.5], 'Color', [0.5, 0.5, 0.5], 'DisplayName', 'Moon Final Position');
-
-% Define the radius of the Moon's Sphere of Influence
-moonSOIRadius = 66100; % Modify this value as needed
-
-% Plot transparent sphere of influence for the Moon at its final position
-[Xsoi, Ysoi, Zsoi] = sphere(50); % Generate sphere data
-
-% Translate sphere data to the final position of the Moon
-Xsoi = Xsoi * moonSOIRadius + Xmoon(1, end);
-Ysoi = Ysoi * moonSOIRadius + Xmoon(2, end);
-Zsoi = Zsoi * moonSOIRadius + Xmoon(3, end);
-
-
-% Plotting Earth as a reference point
-plot3(0, 0, 0, 'bo', 'MarkerFaceColor', 'blue', 'DisplayName', 'Earth');
-
-
-% Plotting the transfer trajectory
-plot3(positions_1(:, 1), positions_1(:, 2), positions_1(:, 3), 'm-', 'LineWidth', 1, 'DisplayName', 'Transfer Trajectory (Earth To Moon)');
-
-% Plotting the transfer trajectory
-plot3(positions_2(:, 1), positions_2(:, 2), positions_2(:, 3), 'm-', 'LineWidth', 1, 'DisplayName', 'Transfer Trajectory (Moon To SOI Exit)');
-
-% Enable data tips and customize them
-dcm_obj = datacursormode(gcf);
-set(dcm_obj,'UpdateFcn',{@myupdatefcn, positions_1, velocities_1});
-
-% Add legends and labels
-legend show;
-xlabel('X (km)');
-ylabel('Y (km)');
-zlabel('Z (km)');
-title('Propagated Transfer Orbit from Earth to SOI Edge Via Lunar Gravity Assist');
-grid on;
-axis equal;
-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;
-
-fprintf('Specific Orbital Energy (epsilon): %f km^2/s^2\n', epsilon);
-fprintf('Turning Angle: %f degrees\n', turning_angle_deg);
-fprintf('Eccentricity: %f\n', eccentricity);
-fprintf('Periapsis Distance (rp): %f km\n', rp);
-
-mu = 4902.8; % Moon's gravitational parameter in km^3/s^2
-a = -mu / 2 * epsilon; % Semi-major axis in km
-e = eccentricity; % Eccentricity
-i = deg2rad(0); % Inclination in radians
-Omega = deg2rad(0);
-
-omega = 0; % Argument of periapsis in radians
-meanAnomalyDeg = 359.9644; % degrees
-
-% Calculate the orbital period (s) of the satellite using Kepler's third law
-P = 2 * pi * sqrt(a^3 / mu);
-
-% Create a time vector 't_vec' with 1000 points spanning one orbital period
-t_vec = linspace(0, P, 100000);
-
-% Calculate the Mean motion (rad/s)
-n = sqrt(mu / a^3);
-
-% Define time vector and orbital elements
-theta_vec = linspace(0, 2*pi, length(t_vec)); % True anomaly from 0 to 2*pi
-
-% Initialize matrix to store state vectors
-X_matrix = zeros(6, length(t_vec));
-
-% Calculate state vectors across the time vector
-for k = 1:length(t_vec)
- X_matrix(:, k) = COE2RV([a, e, i, Omega, omega, theta_vec(k)], mu);
-end
-
-% Extract initial and final state vectors
-X_initial = X_matrix(:, 1);
-
-%% Plotting
-figure;
-hold on;
-
-%% Plot Satellite's trajectory if X_matrix is defined
-if exist('X_matrix', 'var')
- plot3(X_matrix(1, :), X_matrix(2, :), X_matrix(3, :), 'b--', 'LineWidth', 0.5, 'DisplayName', 'Satellite Trajectory');
-end
-
-% Add legends and labels
-legend show;
-xlabel('X (km)');
-ylabel('Y (km)');
-zlabel('Z (km)');
-title('2D Gravity Assist Layout');
-grid on;
-axis equal;
-xlim([-10e4, 10e4]); % Set limits for the x-axis
-ylim([-10e4, 10e4]); % Set limits for the y-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;
-
-%% Clean up by unloading SPICE kernels
-cspice_kclear; % Unload all SPICE kernels and clear the SPICE subsystem
\ No newline at end of file