diff --git a/Total_Delta_V_No_Lunar_Assist.m b/Total_Delta_V_No_Lunar_Assist.m
new file mode 100644
index 0000000000000000000000000000000000000000..90efc29bc7f6b14ac8d57e9ec0e2c0bd553f497f
--- /dev/null
+++ b/Total_Delta_V_No_Lunar_Assist.m
@@ -0,0 +1,526 @@
+%% Function Description (Trajectory Computation between Earth and Mars using Patch Conic Approximation)
+% This script computes the trajectory of a spacecraft traveling between Earth and Mars 
+% using the patch conic approximation. The patch conic approximation simplifies the 
+% problem by dividing the spacecraft's trajectory into separate conic sections, each governed 
+% by the gravitational influence of a single celestial body (Earth, Sun, or Mars). 
+% This script does not utilize a gravity assist; instead, it calculates a direct transfer trajectory 
+% from Earth to 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 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 relevant celestial bodies.
+% - Calculating state vectors (position and velocity) for Earth and the Moon over a specified observation period 
+%   using SPICE functions. These vectors help determine the spacecraft's trajectory relative to Earth and the Sun.
+% - Using Kepler's third law to compute the orbital period of the spacecraft's initial orbit around Earth and 
+%   propagating the orbit over one period to obtain the spacecraft's state vectors throughout the trajectory.
+% - Utilizing a Lambert solver to determine the initial and final velocity vectors required for the interplanetary 
+%   transfer from Earth to Mars. The Lambert solver accounts for the gravitational influence of the Sun and calculates 
+%   the optimal transfer trajectory.
+% - Iterating over a range of times of flight (TOFs) to find the optimal TOF that minimizes the delta-V (change in velocity) 
+%   required for the transfer. The script evaluates each TOF and selects the one that results in the lowest delta-V, 
+%   indicating the most efficient trajectory.
+% - Propagating the trajectory from Earth's parking orbit to the edge of Earth's Sphere of Influence (SOI) using numerical 
+%   integration (ode45). The output state vectors are used to plot the spacecraft's trajectory and visualize the transfer.
+% - Plotting the results, including the Earth's position, the spacecraft's trajectory, and other relevant elements, to 
+%   provide a visual representation of the transfer orbit from Earth to Mars.
+
+% 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 direct transfer trajectories for 
+% interplanetary missions without utilizing 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;
+
+%% Constants
+% Define the conversion factor from kilometers to Astronomical Units (AU)
+kmToAU = 1 / 149597870.7; % Inverse of the average distance from Earth 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 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(310.445);  % RAAN in radians 90deg + 51.82
+omega = 0;  % Argument of periapsis in radians
+meanAnomalyDeg = 0; % degrees
+mu = 3.986e5;  % Earth's gravitational parameter in km^3/s^2
+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.
+
+
+% Standard gravitational parameter of the Sun (km^3/s^2)
+mu_sun = 1.32712440018e11; % Sun's gravitational constant, affects orbital dynamics
+
+%% 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.)
+
+%% Earth moon State Vectors
+% Define start and end date for the observation period
+StartDate = '1 Nov 2026 20:04:03 (UTC)'; % Update this line with your desired start date
+EndDate = '5 Nov 2026 01:03:34 (UTC)'; % Update this line with your desired end date
+
+% Convert both dates to Ephemeris Time using SPICE's str2et function
+et0 = cspice_str2et(StartDate); % Converts the StartDate to ephemeris time (seconds past J2000)
+et_end = cspice_str2et(EndDate); % Converts the EndDate to ephemeris time (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
+
+% Retrieve state vectors for each time point using a loop
+for tt = 1:length(et)
+    % Get the state vector of moon from the SPICE toolkit
+    Xmoon(:, tt) = cspice_spkezr('301', et(tt), 'ECLIPJ2000', 'NONE', '399');
+    Xearth(:, tt) = cspice_spkezr('399', et(tt), 'ECLIPJ2000', 'NONE', '10');
+    % Get the state vector of Earth from the SPICE toolkit
+    % Convert from km to AU using conversion factor
+end
+
+% 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);
+
+% Print initial and final state vectors
+fprintf('Initial X:\n');
+disp(X_initial);
+
+% 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
+
+%% Desired final velocity at Mars
+v_mars_required = [-1.795014; 2.364287; 0.612314];
+v_mars_required_mag = norm(v_mars_required);  % Calculate magnitude of v_mars_required
+
+%% 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, idx90);  % Extract Earth's position vector at launch
+
+% Spacecraft final position vector at Mars
+rf = [-548056.087908; 720190.384622; 186378.606808];  % Extract Mars' position vector at the arrival time
+
+%% Lambert Solver
+% Define the range of TOFs
+tof_days = 3:1/86400:3.5;  % TOFs from 1 day to 8 days in 30-minute increments
+tof_seconds = tof_days * 86400;  % Convert days to seconds
+
+% Preallocate arrays to store initial and final velocity vectors for each TOF
+v0_arr = zeros(3, length(tof_seconds));
+vf_arr = zeros(3, length(tof_seconds));
+
+% Loop through each TOF
+for i = 1:length(tof_seconds)
+    % Call the Lambert solver for each TOF
+    [v0, vf] = lambert_solver_moon(r0, rf, tof_seconds(i), mu, DM);
+
+    % Store computed velocities
+    v0_arr(:, i) = v0;
+    vf_arr(:, i) = vf;
+
+    % Calculate delta v for each TOF
+    delta_v_arr(i) = norm(vf - v_mars_required);
+end
+
+% Find the minimum delta v and corresponding TOF
+[min_delta_v, min_index] = min(delta_v_arr);
+best_tof = tof_days(min_index);
+best_v0 = v0_arr(:, min_index) - X_matrix(4:6, 1);
+best_vf = vf_arr(:, min_index);
+best_v0_mag = norm(best_v0);
+best_vf_mag = norm(best_vf);
+
+% Display the best TOF and the minimum delta v
+fprintf('\nBest TOF: %.2f days with a delta v of %f km/s\n', best_tof, min_delta_v);
+fprintf('Corresponding initial velocity v0: [%f, %f, %f] with magnitude: %f km/s\n', ...
+        best_v0(1), best_v0(2), best_v0(3), best_v0_mag);
+fprintf('Corresponding final velocity vf: [%f, %f, %f] with magnitude: %f km/s\n', ...
+        best_vf(1), best_vf(2), best_vf(3), best_vf_mag);
+fprintf('Desired Mars velocity v_mars_required: [%f, %f, %f] with magnitude: %f km/s\n', ...
+        v_mars_required(1), v_mars_required(2), v_mars_required(3), v_mars_required_mag);
+
+v0 = v0_arr(:, min_index);
+
+% Define the time span of the simulation
+tspan = [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 = [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, state_out] = ode45(@(t, y) Orbit_Propagation(t, y, mu), tspan, initial_state, options);  % Numerically integrate the trajectory
+
+% Extract the position vectors from the output state for plotting
+positions = state_out(:, 1:3);  % Extract position data from the output state for subsequent use in plotting
+velocities = state_out(:, 4:6);  % Extract velocity data from the output state for subsequent use in plotting
+
+delta_v_final = 1.109200;
+
+% Display the delta-v
+disp('Delta-v for the inital transfer (km/s):');
+disp(best_v0_mag);  % Output the computed delta-v value
+
+% 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 + best_v0_mag);  % Output the computed delta-v value
+%% 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 Trajectory');
+    plot3(X_matrix(1, 1), X_matrix(2, 1), X_matrix(3, 1), 'g.', 'MarkerSize', 10, 'DisplayName', 'Satellite Initial Position');
+end
+
+% Plot transparent sphere of influence for the Moon at its final position
+[Xsoi, Ysoi, Zsoi] = sphere(50); % Generate sphere data
+
+% Plotting Earth as a reference point
+plot3(0, 0, 0, 'bo', 'MarkerFaceColor', 'blue', 'DisplayName', 'Earth');
+
+% Plot Moon's final position
+plot3(rf(1), rf(2), rf(3), 'o', 'MarkerFaceColor', [1, 0, 0], 'Color', [1, 0, 0.5], 'DisplayName', 'Spacecrafts Position');
+
+% Plotting the transfer trajectory
+plot3(positions(:, 1), positions(:, 2), positions(:, 3), 'm-', 'LineWidth', 1, 'DisplayName', 'Transfer Trajectory');
+
+% 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 the translated sphere of influence
+surf(Xsoi, Ysoi, Zsoi, 'EdgeColor', 'none', 'FaceColor', 'blue', 'FaceAlpha', 0.1, 'DisplayName', 'Moon Sphere of Influence');
+
+% Enable data tips and customize them
+dcm_obj = datacursormode(gcf);
+set(dcm_obj,'UpdateFcn',{@myupdatefcn, positions, velocities});
+
+% Add legends and labels
+legend show;
+xlabel('X (km)');
+ylabel('Y (km)');
+zlabel('Z (km)');
+title('Hyperbolic Transfer Orbit - From Earth Parking Orbit to SOI Boundary');
+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;
+
+V_Helio = [-19.108237304219; 26.8603108423787; 0.492875064188567];
+
+V_Earth_Helio = Xearth(4:6, end);
+
+% Convert the object's velocity to the geocentric frame
+object_geocentric_velocity = V_Helio - V_Earth_Helio;
+
+%% 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, idx90);  % Extract Earth's position vector at launch
+initial_velocity = X_matrix(4:6, idx90);  % Initial velocity vector of the orbit around Earth
+
+% Spacecraft final position vector at the Moon
+rf = [-534253.3462; 753776.5515 ; 13824.51521];  % Example final position vector
+
+%% Iterate over different times of flight
+days = linspace(2, 20, 100000);  % Time of flight from 1 day to 20 days with 100 steps for smoother plot
+total_delta_v_values = zeros(size(days));  % Array to store total delta_v values
+
+for idx = 1:length(days)
+    d = days(idx);
+    delta_t = d * 86400;  % Convert days to seconds
+    
+    % Call the Lambert solver function with the current time of flight
+    [v0, vf, exitflag] = lambert_solver(r0, rf, delta_t, mu, DM);
+    
+    if exitflag == 1
+        % Calculate the first delta-v (from parking orbit to hyperbolic trajectory)
+        delta_v1 = norm(v0 - initial_velocity);
+        
+        % Calculate the second delta-v (adjust velocity at the edge of SOI to match Mars trajectory)
+        delta_v2 = norm(vf - v_mars_required);
+        
+        % Sum the delta-v values to get the total delta-v
+        total_delta_v = delta_v1 + delta_v2;
+        
+        % Store the total delta-v value for the current TOF
+        total_delta_v_values(idx) = total_delta_v;
+    else
+        % If the solver didn't converge, set total_delta_v to Inf
+        total_delta_v_values(idx) = Inf;
+    end
+end
+
+%% 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 Trajectory');
+    plot3(X_matrix(1, 1), X_matrix(2, 1), X_matrix(3, 1), 'g.', 'MarkerSize', 10, 'DisplayName', 'Satellite Initial Position');
+end
+
+% Plot transparent sphere of influence for the Moon at its final position
+[Xsoi, Ysoi, Zsoi] = sphere(50); % Generate sphere data
+
+% Plotting Earth as a reference point
+plot3(0, 0, 0, 'bo', 'MarkerFaceColor', 'blue', 'DisplayName', 'Earth');
+% Spacecraft final position vector at Mars
+rf = [-548056.087908; 720190.384622; 186378.606808];  % Extract Mars' position vector at the arrival time
+
+% Plot Moon's final position
+plot3(rf(1), rf(2), rf(3), 'o', 'MarkerFaceColor', [1, 0, 0], 'Color', [1, 0, 0.5], 'DisplayName', 'Spacecrafts Position');
+
+% Plotting the transfer trajectory
+plot3(positions(:, 1), positions(:, 2), positions(:, 3), 'm-', 'LineWidth', 1, 'DisplayName', 'Transfer Trajectory');
+
+% Enable data tips and customize them
+dcm_obj = datacursormode(gcf);
+set(dcm_obj,'UpdateFcn',{@myupdatefcn, positions, velocities});
+
+% Add legends and labels
+legend show;
+xlabel('X (km)');
+ylabel('Y (km)');
+zlabel('Z (km)');
+title('Propagated Transfer Orbit from Earth to Moon');
+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;
+
+%% 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 Trajectory');
+    plot3(X_matrix(1, 1), X_matrix(2, 1), X_matrix(3, 1), 'g.', 'MarkerSize', 10, 'DisplayName', 'Satellite Initial Position');
+end
+
+% Plot transparent sphere of influence for the Moon at its final position
+[Xsoi, Ysoi, Zsoi] = sphere(50); % Generate sphere data
+
+% Plotting Earth as a reference point
+plot3(0, 0, 0, 'bo', 'MarkerFaceColor', 'blue', 'DisplayName', 'Earth');
+% Spacecraft final position vector at Mars
+rf = [-548056.087908; 720190.384622; 186378.606808];  % Extract Mars' position vector at the arrival time
+
+% Plot Moon's final position
+plot3(rf(1), rf(2), rf(3), 'o', 'MarkerFaceColor', [1, 0, 0], 'Color', [1, 0, 0.5], 'DisplayName', 'Spacecrafts Position');
+
+% Plotting the transfer trajectory
+plot3(positions(:, 1), positions(:, 2), positions(:, 3), 'm-', 'LineWidth', 1, 'DisplayName', 'Transfer Trajectory');
+
+% Enable data tips and customize them
+dcm_obj = datacursormode(gcf);
+set(dcm_obj,'UpdateFcn',{@myupdatefcn, positions, velocities});
+
+% Add legends and labels
+legend show;
+xlabel('X (km)');
+ylabel('Y (km)');
+zlabel('Z (km)');
+title('Hyperbolic Transfer Orbit - From Earth Parking Orbit to SOI Boundary');
+grid on;
+axis equal;
+xlim([-6e5, 1e5]);  % Set limits for the x-axis
+ylim([-1e4, 8e5]);  % Set limits for the y-axis
+zlim([-2e5, 2e5]); % 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;
+
+%% 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', 2, 'DisplayName', 'Satellite Trajectory');
+    plot3(X_matrix(1, 1), X_matrix(2, 1), X_matrix(3, 1), 'g.', 'MarkerSize', 30, 'DisplayName', 'Satellite Initial Position');
+end
+
+% Plot transparent sphere of influence for the Moon at its final position
+[Xsoi, Ysoi, Zsoi] = sphere(50); % Generate sphere data
+
+% Plotting Earth as a reference point
+plot3(0, 0, 0, 'bo', 'MarkerFaceColor', 'blue', 'DisplayName', 'Earth');
+% Spacecraft final position vector at Mars
+rf = [-548056.087908; 720190.384622; 186378.606808];  % Extract Mars' position vector at the arrival time
+
+% Plot Moon's final position
+plot3(rf(1), rf(2), rf(3), 'o', 'MarkerFaceColor', [1, 0, 0], 'Color', [1, 0, 0.5], 'DisplayName', 'Spacecrafts Position');
+
+% Plotting the transfer trajectory
+plot3(positions(:, 1), positions(:, 2), positions(:, 3), 'm-', 'LineWidth', 2, 'DisplayName', 'Transfer Trajectory');
+
+% Enable data tips and customize them
+dcm_obj = datacursormode(gcf);
+set(dcm_obj,'UpdateFcn',{@myupdatefcn, positions, velocities});
+
+% Add legends and labels
+legend show;
+xlabel('X (km)');
+ylabel('Y (km)');
+zlabel('Z (km)');
+title('Propagated Transfer Orbit from Earth to Moon');
+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(310.445 - 270, 0);  % 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
+
+function txt = myupdatefcn(~, event_obj, positions, velocities)
+    % Customizes text of data tips
+    pos = get(event_obj,'Position');  % Get the position of the data tip
+    idx = find(positions(:,1) == pos(1) & positions(:,2) == pos(2) & positions(:,3) == pos(3));
+    if isempty(idx)
+        txt = {'Position:', pos, 'Velocity: Not found'};
+    else
+        vel = velocities(idx, :);  % Extract corresponding velocity
+        txt = {['Position: (', num2str(pos(1)), ', ', num2str(pos(2)), ', ', num2str(pos(3)), ')'], ...
+               ['Velocity: (', num2str(vel(1)), ', ', num2str(vel(2)), ', ', num2str(vel(3)), ')']};
+    end
+end