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