diff --git a/Trajectory_Earth_Mars.m b/Trajectory_Earth_Mars.m
new file mode 100644
index 0000000000000000000000000000000000000000..3014e4b808d95b87d8e81d523da845f1ba7513ef
--- /dev/null
+++ b/Trajectory_Earth_Mars.m
@@ -0,0 +1,330 @@
+%% Script Description (Earth to Mars Trajectory)
+% This MATLAB script calculates and visualizes the interplanetary trajectory of a spacecraft transferring from Earth to Mars.
+% It leverages the SPICE toolkit to handle precise celestial mechanics computations, including state vectors of Earth and Mars 
+% for a given departure and arrival date. The script uses these vectors to solve Lambert's problem to determine the required 
+% velocities for departure from Earth and arrival at Mars, integrating these velocities to plot the spacecraft's trajectory.
+%
+% Key components of the script include:
+% - Conversion factors for distances (km to AU) and angles (degrees to radians).
+% - Definitions of orbital parameters for Earth and Mars.
+% - SPICE toolkit setup for retrieving planetary data.
+% - Implementation of the Lambert solver to compute initial and final velocity vectors needed for the transfer.
+% - Numerical integration (using `ode45`) of the spacecraft's trajectory from Earth to Mars.
+% - Visualization of the orbits of Earth and Mars, as well as the spacecraft trajectory in a 3D plot.
+%
+% Constants
+% Includes standard gravitational parameters and semi-major axes for Earth and Mars. These constants are essential for
+% the calculations and conversions throughout the script.
+%
+% Calling SPICE and Functions
+% Sets up paths and loads necessary SPICE kernels for planetary data access, enabling the extraction of precise 
+% celestial positions and velocities.
+%
+% Lambert Solver
+% Calculates the necessary velocities at departure and arrival using the Lambert problem solution, which is critical 
+% for determining the feasible path the spacecraft will follow.
+%
+% Trajectory Propagation
+% Uses `ode45` to simulate the spacecraft's orbit based on the initial conditions derived from the Lambert solver.
+% This step is vital for understanding how the spacecraft will move through space under the influence of celestial 
+% gravities.
+%
+% Author: Thomas West
+% Date: April 18, 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 Mars orbits in km (average distances)
+a_earth = 149597870.7; % 1 AU in km, average distance from the Earth to the Sun
+a_mars = 227939200; % Approximately 1.524 AU in km, average distance from Mars to the Sun
+
+% Standard gravitational parameter of the Sun (km^3/s^2)
+mu = 1.32712440018e11; % 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 = '1 Nov 2026 12:30:00 (UTC)';  % Define the launch date in UTC
+ArrivalDate = '7 Sep 2027 12:30:00 (UTC)'; % Define the arrival date in UTC
+
+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;  % Compute the time difference in seconds
+
+% Calculate the number of days between the start and end date
+num_days = total_duration; % Convert the duration from seconds to days
+
+% Round num_days to the nearest whole number to avoid fractional days
+num_days = round(num_days);  % Round the number of days to the nearest integer
+
+% Generate a time vector from start to end date using the calculated number of days
+et = linspace(et0, et_end, num_days + 1); % Generate evenly spaced time points including both endpoints
+
+% Retrieve state vectors for each time point
+for tt = 1:length(et)  % Loop over each time point
+    Xmars(:, tt) = cspice_spkezr('4', et(tt), 'ECLIPJ2000', 'NONE', '10');  % Retrieve Mars' state vector at time tt
+    Xearth(:, tt) = cspice_spkezr('399', et(tt), 'ECLIPJ2000', 'NONE', '10');  % Retrieve Earth's state vector at time tt
+    % Convert from km to AU
+    Xmars_AU(:, tt) = Xmars(:, tt) * kmToAU;  % Convert Mars' position from km to Astronomical Units
+    Xearth_AU(:, tt) = Xearth(:, tt) * kmToAU; % Convert Earth's position from km to Astronomical Units
+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 = Xearth(1:3, 1);  % Extract Earth's position vector at launch
+
+% Spacecraft final position vector at Mars
+rf = Xmars(1:3, end);  % Extract Mars' position vector at the arrival time
+
+% 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
+
+% Output the initial and final velocities in km/s
+disp('Initial velocity vector (km/s):');  % Display message for initial velocity
+disp(v0);  % Display the computed initial velocity vector
+disp('Final velocity vector (km/s):');  % Display message for final velocity
+disp(vf);  % Display the computed final velocity vector
+
+% Define the time span of the simulation with more points
+tspan = linspace(0, delta_t, num_days); % Generate a finely spaced time vector 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);
+
+%% Calculate Distances and Find Exit from Earth's SOI
+SOI_distance_km = 924000;  % Define the distance for SOI exit
+exit_index = 0;  % Initialize exit index
+
+% Loop through all time points to calculate distances
+for tt = 1:length(t_out)
+    % Calculate the distance between Earth and the spacecraft at time tt
+    distance = norm(positions(tt, :) - Xearth(1:3, tt)');  % Make sure Xearth is interpolated if necessary
+    
+    % Check if the distance is greater than the SOI distance
+    if distance > SOI_distance_km
+        exit_index = tt;  % Save the index when spacecraft exits the SOI
+        break;  % Exit the loop as we found the first exit point
+    end
+end
+
+if exit_index > 0
+    % Earth's state vector at exit time (assuming Earth data is interpolated to match spacecraft data points)
+    Earth_position = Xearth(1:3, exit_index);  % Earth's position vector in km
+    Earth_velocity = Xearth(4:6, exit_index);  % Earth's velocity vector in km/s
+    
+    % Spacecraft's heliocentric state vector at exit time
+    Spacecraft_position = positions(exit_index, :);  % Spacecraft's position vector in km
+    Spacecraft_velocity = velocities(exit_index, :);  % Spacecraft's velocity vector in km/s
+    
+    % Convert to geocentric frame
+    Geocentric_position = Spacecraft_position - Earth_position';  % Convert position from heliocentric to geocentric
+    Geocentric_velocity = Spacecraft_velocity - Earth_velocity';  % Convert velocity from heliocentric to geocentric
+
+    % Time when the spacecraft exits the SOI in seconds and hours
+    exit_time_seconds = t_out(exit_index);  % Time in seconds since the start of the simulation
+    exit_time_hours = exit_time_seconds / 3600;  % Convert time to hours
+    exit_time_days = exit_time_seconds / 86400;  % Convert time to days
+
+    % Display the results
+    fprintf('Time of exit from Earths SOI: %f seconds (%f hours, %f days)\n', exit_time_seconds, exit_time_hours, exit_time_days);
+    fprintf('Earths position vector at exit time (km): [%f, %f, %f]\n', Earth_position);
+    fprintf('Earths velocity vector at exit time (km): [%f, %f, %f]\n', Earth_velocity);
+    fprintf('Spacecrafts heliocentric position at exit time (km): [%f, %f, %f]\n', Spacecraft_position);
+    fprintf('Spacecrafts heliocentric velocity at exit time (km/s): [%f, %f, %f]\n', Spacecraft_velocity);
+    fprintf('Geocentric position vector of the spacecraft at exit time (km): [%f, %f, %f]\n', Geocentric_position);
+    fprintf('Geocentric velocity vector of the spacecraft at exit time (km/s): [%f, %f, %f]\n', Geocentric_velocity);
+else
+    fprintf('Spacecraft does not exit Earths SOI within the simulation period.\n');
+end
+
+%$ Mars proximity threshold in kilometers
+Mars_proximity_km = 577300;  % Define the proximity distance for Mars approach
+approach_index = 0;  % Initialize approach index
+
+% Loop through all time points to calculate distances
+for tt = 1:length(t_out)
+    % Calculate the distance between Mars and the spacecraft at time tt
+    distance = norm(positions(tt, :) - Xmars(1:3, tt)');  % Ensure Xmars is interpolated
+    
+    % Check if the distance is less than the Mars proximity distance
+    if distance < Mars_proximity_km
+        approach_index = tt;  % Save the index when spacecraft approaches Mars
+        break;  % Exit the loop as we found the first approach point
+    end
+end
+
+% Output results
+if approach_index > 0
+    approach_time_seconds = t_out(approach_index);  % Time in seconds since the start of the simulation
+    approach_time_hours = approach_time_seconds / 3600;  % Convert time to hours
+    approach_time_days = approach_time_seconds / 86400;  % Convert time to days
+    
+    % Extract Mars and spacecraft state vectors at approach time
+    Mars_position = Xmars(1:3, approach_index);  % Mars' position vector in km
+    Mars_velocity = Xmars(4:6, approach_index);  % Mars' velocity vector in km/s
+    Spacecraft_position = positions(approach_index, :);  % Spacecraft's position vector in km
+    Spacecraft_velocity = velocities(approach_index, :);  % Spacecraft's velocity vector in km/s
+    
+    % Convert to areocentric frame
+    Areocentric_position = Spacecraft_position - Mars_position';  % Convert position from heliocentric to areocentric
+    Areocentric_velocity = Spacecraft_velocity - Mars_velocity';  % Convert velocity from heliocentric to areocentric
+    
+    fprintf('Spacecraft approaches Mars within %f km at t = %f seconds (%f hours, %f days)\n', Mars_proximity_km, approach_time_seconds, approach_time_hours, approach_time_days);
+    fprintf('Mars position vector at entry time (km): [%f, %f, %f]\n', Mars_position);
+    fprintf('Mars velocity vector at entry time (km): [%f, %f, %f]\n', Mars_velocity);
+    fprintf('Spacecrafts heliocentric position at entry time (km): [%f, %f, %f]\n', Spacecraft_position);
+    fprintf('Spacecrafts heliocentric velocity at entry time (km/s): [%f, %f, %f]\n', Spacecraft_velocity);
+    fprintf('Areocentric position vector of the spacecraft at approach time (km): [%f, %f, %f]\n', Areocentric_position);
+    fprintf('Areocentric velocity vector of the spacecraft at approach time (km/s): [%f, %f, %f]\n', Areocentric_velocity);
+else
+    fprintf('Spacecraft does not approach Mars within %f km during the simulation period.\n', Mars_proximity_km);
+end
+
+% Assuming Xearth and v0 are defined
+% Extract the vector from Xearth
+vec1 = Xearth(4:6, 1);
+
+% Assume v0 is already defined as a vector
+vec2 = v0;
+
+% Calculate the dot product of vec1 and vec2
+dot_product = dot(vec1, vec2);
+
+% Calculate the magnitudes of vec1 and vec2
+mag_vec1 = norm(vec1);
+mag_vec2 = norm(vec2);
+
+% Calculate the cosine of the angle
+cos_theta = dot_product / (mag_vec1 * mag_vec2);
+
+% Calculate the angle in degrees
+angle = acosd(cos_theta);
+
+% Display the angle
+disp(['The angle between the vectors is ', num2str(angle), ' degrees.']);
+
+v0
+Xearth(4:6, 1)
+
+%% Free View Plot
+% Use start and end dates from earlier in the script for labeling
+StartDate = LaunchDate;
+EndDate = ArrivalDate;
+
+% Get the initial and final positions for Earth and Mars in Astronomical Units
+earthInitialPos = Xearth_AU(1:3, 1);  % Extract the initial position of Earth
+earthFinalPos = Xearth_AU(1:3, end);  % Extract the final position of Earth
+marsInitialPos = Xmars_AU(1:3, 1);    % Extract the initial position of Mars
+marsFinalPos = Xmars_AU(1:3, end);    % Extract the final position of Mars
+
+% Create a new figure for plotting
+figure();
+% Plot Earth's orbit using blue lines
+h_earthOrbit = plot3(Xearth_AU(1, :), Xearth_AU(2, :), Xearth_AU(3, :), 'b', 'DisplayName', 'Earth Orbit'); hold on;
+% Mark Earth's final position with a blue circle
+h_earthPosition = plot3(Xearth_AU(1, 1), Xearth_AU(2, 1), Xearth_AU(3, 1), 'bo', 'MarkerFaceColor', 'blue', 'DisplayName', sprintf('Earth'));
+
+% Plot Mars's orbit using red lines
+h_marsOrbit = plot3(Xmars_AU(1, :), Xmars_AU(2, :), Xmars_AU(3, :), 'r', 'DisplayName', 'Mars Orbit');
+% Mark Mars's final position with a red circle
+h_marsPosition = plot3(Xmars_AU(1, 1), Xmars_AU(2, 1), Xmars_AU(3, 1), 'ro', 'MarkerFaceColor', 'red', 'DisplayName', sprintf('Mars'));
+
+% Define an offset value to position the text above the marker
+offset = 0.05; % adjust the offset value as needed
+
+% Mark Earth's initial position with a blue circle and add text annotation
+h_earthInitial = plot3(earthInitialPos(1), earthInitialPos(2), earthInitialPos(3), 'bo', 'MarkerFaceColor', 'blue', 'DisplayName', 'Earth Initial Position');
+text(earthInitialPos(1), earthInitialPos(2), earthInitialPos(3) + offset, 'Earth Initial', 'VerticalAlignment', 'cap', 'HorizontalAlignment', 'center', 'Color', 'blue');
+
+% Earth Final Position Marker and Text
+h_earthFinal = plot3(earthFinalPos(1), earthFinalPos(2), earthFinalPos(3), 'bo', 'MarkerFaceColor', 'blue', 'DisplayName', 'Earth Final Position');
+text(earthFinalPos(1) + offset, earthFinalPos(2) - offset, earthFinalPos(3), 'Earth Final', 'VerticalAlignment', 'top', 'HorizontalAlignment', 'left', 'Color', 'blue');
+
+% Mars Initial Position Marker and Text
+h_marsInitial = plot3(marsInitialPos(1), marsInitialPos(2), marsInitialPos(3), 'ro', 'MarkerFaceColor', 'red', 'DisplayName', 'Mars Initial Position');
+text(marsInitialPos(1) + offset, marsInitialPos(2) + offset, marsInitialPos(3), 'Mars Initial', 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'left', 'Color', 'red');
+
+% Mars Final Position Marker and Text
+h_marsFinal = plot3(marsFinalPos(1), marsFinalPos(2), marsFinalPos(3), 'ro', 'MarkerFaceColor', 'red', 'DisplayName', 'Mars Final Position');
+text(marsFinalPos(1), marsFinalPos(2) - offset, marsFinalPos(3) + offset, 'Mars Final', 'VerticalAlignment', 'top', 'HorizontalAlignment', 'center', 'Color', 'red');
+
+% Plot the spacecraft's transfer trajectory
+% Convert positions from km to AU for plotting
+positions_AU = positions * kmToAU;
+h_spacecraftTrajectory = plot3(positions_AU(:, 1), positions_AU(:, 2), positions_AU(:, 3), 'm-', 'LineWidth', 1, 'DisplayName', 'Spacecraft Trajectory');
+
+% Plot the Sun at the origin with a legend handle
+h_sun = plot3(0, 0, 0, 'yo', 'MarkerSize', 10, 'MarkerFaceColor', 'yellow', 'DisplayName', 'Sun');
+
+% Determine the title string based on DM
+if DM == 1  % Check if DM equals 1 (short way)
+    dmTitle = 'DM = +1 (Short Way Transfer)';  % Assign the title for short way transfer
+else  % If DM does not equal 1, it must be -1 (long way)
+    dmTitle = 'DM = -1 (Long Way Transfer)';  % Assign the title for long way transfer
+end
+
+% Enhancements for 3D Visualization
+xlabel('X Distance (AU)');  % Label the x-axis of the plot
+ylabel('Y Distance (AU)');  % Label the y-axis of the plot
+zlabel('Z Distance (AU)');  % Label the z-axis of the plot
+title(sprintf('Spacecraft Transfer Trajectory from Earth to Mars (%s to %s) (%s)', LaunchDate, ArrivalDate, dmTitle));  % Title the plot with the mission timeframe and direction of motion
+legend([h_earthOrbit, h_marsOrbit, h_earthPosition, h_marsPosition, h_spacecraftTrajectory, h_sun], 'Location', 'bestoutside');  % Create a legend and place it outside the plot
+grid on;  % Enable grid for better visibility of plot scale and distances
+axis equal;  % Set all axis scales to be equal
+xlim([-2, 2]);  % Set limits for the x-axis
+ylim([-2, 2]);  % Set limits for the y-axis
+zlim([-2, 2]); % Set limits for the z-axis
+view(3); % Set the view angle to 3D for a better perspective of the plot
+
+% Ensure the plot holds for additional plotting
+hold off; % Release the hold on the current plot
+
+%% Clean up by unloading SPICE kernels
+cspice_kclear;  % Unload all SPICE kernels and clear the SPICE subsystem