diff --git a/V_Inf_Globe_Iterative.m b/V_Inf_Globe_Iterative.m
new file mode 100644
index 0000000000000000000000000000000000000000..b02a00f1e4fa0d1b4a1d28d129b5d0f3916e35dd
--- /dev/null
+++ b/V_Inf_Globe_Iterative.m
@@ -0,0 +1,167 @@
+%% Function Description (Optimal Gravity Assist Calculation using V-Infinity Globe)
+% This script calculates the optimal input variables to optimize a gravity
+% assist maneuver around the Moon to reach a desired final orbit or destination.
+% The optimization process involves minimizing the required delta-V (change in velocity)
+% for a spacecraft to transition from one orbit to another, using the Moon's gravity
+% to assist in the maneuver.
+%
+% The script utilizes the SPICE toolkit to load planetary and lunar ephemerides
+% and gravitational models, providing accurate state vectors for celestial bodies.
+% It then iterates over a range of possible arrival and departure dates to find
+% the optimal conditions for the gravity assist.
+%
+% Key steps and calculations include:
+% - Loading SPICE kernels for accurate celestial mechanics data, including
+% planetary positions, velocities, and gravitational parameters.
+% - Defining the range and resolution for the angular variables alpha and kappa,
+% which describe the direction of the incoming and outgoing hyperbolic excess
+% velocity vectors (v-infinity) relative to the Moon.
+% - Using a Lambert solver to compute the initial and final velocity vectors
+% required to transfer between Earth orbit and a specified point near the Moon,
+% considering Earth's gravitational parameter and the desired transfer mode
+% (e.g., Type I transfer for shorter paths).
+% - Calculating the new v-infinity vector after the gravity assist, taking into
+% account the Moon's gravity and the spacecraft's initial approach velocity.
+% - Computing the turning angle of the spacecraft's trajectory during the gravity
+% assist and evaluating specific orbital energy to ensure the spacecraft does
+% not crash into the Moon or fail to escape Earth's sphere of influence.
+% - Iteratively searching over a range of potential time of flights (TOFs) and
+% arrival dates to identify the optimal set of conditions that minimize the
+% required delta-V for the desired gravity assist maneuver.
+%
+% The script outputs the optimal departure date, arrival date, minimum delta-V
+% required, and the optimal time of flight (TOF) to achieve the target gravity
+% assist maneuver.
+%
+% 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;
+
+%% 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.)
+
+mu = 3.986e5; % Earth's gravitational parameter in km^3/s^2
+
+%% Constants and Configuration
+mu_Earth = 398600.4418; % Earth's gravitational parameter in km^3/s^2
+mu_Moon = 4902.8; % Moon's gravitational parameter in km^3/s^2
+critical_rp = 1737.4 + 50; % Critical radius of periapsis in km
+
+%% Initialize Variables
+DM = 1; % Transfer mode (+1 for Type I, shorter path)
+minDeltaV = inf; % Initialize minimum delta-v to a large number
+optimalDate = ''; % To store the date corresponding to the minimum delta-v
+optimalTOF = 0; % To store the optimal time of flight
+
+%% Define the range and resolution of alpha and kappa
+alpha_range = -180:1:180;
+kappa_range = -90:1:90;
+v_mars_required = [-1.973070; 2.268442; 0.492901]; % Required Mars velocity vector
+
+%% Date Range for the Simulation
+startDate = datetime(2026, 10, 26, 00, 0, 0);
+endDate = datetime(2026, 10, 26, 23, 0, 0);
+
+%% Prepare to Store Results Hourly
+results = struct('Date', {}, 'MinDeltaV', {}, 'TurningAngle', {});
+
+%% Loop over each day within the range
+for currentArrivalDate = startDate:hours(1):endDate
+ % Loop over TOF from 1 day to 2.5 days in 1-hour increments
+ for tofHours = 43:1:44 % From 24 hours to 60 hours
+ currentLaunchDate = currentArrivalDate - hours(tofHours);
+
+ % Convert the current dates to string in the required format
+ LaunchDate = datestr(currentLaunchDate, 'dd mmm yyyy HH:MM:SS (UTC)');
+ ArrivalDate = datestr(currentArrivalDate, 'dd mmm yyyy HH:MM:SS (UTC)');
+
+ % Convert launch and arrival dates to Ephemeris Time
+ et0 = cspice_str2et(LaunchDate);
+ et_end = cspice_str2et(ArrivalDate);
+ delta_t = et_end - et0; % Transfer time in seconds
+
+ % Get Moon's state vector at arrival
+ Xmoon = cspice_spkezr('301', et_end, 'ECLIPJ2000', 'NONE', '399');
+ r_enc_vector = Xmoon(1:3); % Position vector of the Moon
+ v_ga_vector = Xmoon(4:6); % Gravity assist velocity vector
+
+ % Normalize vectors for q1, q2, q3
+ q2 = v_ga_vector / norm(v_ga_vector);
+ q3 = cross(r_enc_vector, v_ga_vector);
+ q3 = q3 / norm(q3); % Normalize to make it a unit vector
+ q1 = cross(q2, q3); % This is already unit length if q2 and q3 are unit vectors
+
+ % Lambert solver to find initial and final velocity vectors
+ [v0, vf] = lambert_solver_moon([0; 0; 6800], r_enc_vector, delta_t, mu_Earth, DM);
+
+ v_inf_vector = vf - v_ga_vector; % Calculate the departure velocity vector
+ v_inf = norm(v_inf_vector); % Magnitude of departure velocity vector
+
+ % Initialize DeltaVGrid for current date and TOF
+ DeltaVGrid = zeros(length(kappa_range), length(alpha_range));
+
+ % Loop over all alpha and kappa combinations
+ for i = 1:length(alpha_range)
+ alpha = deg2rad(alpha_range(i));
+ for j = 1:length(kappa_range)
+ kappa = deg2rad(kappa_range(j));
+
+ % Calculate new v_infinity after gravity assist using q1, q2, q3
+ v_infinity_new = v_inf * (sin(alpha) * cos(kappa) * q1 + ...
+ cos(alpha) * q2 + ...
+ sin(alpha) * sin(kappa) * q3);
+
+ % Calculate geocentric velocity
+ 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
+
+ % Check for crashing into the Moon or not escaping Earth's SOI
+ 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));
+
+ 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_mars_required);
+ end
+ end
+ end
+
+ % Find the minimum delta-V for this TOF
+ [currentMinDeltaV, ~] = min(DeltaVGrid(:), [], 'omitnan');
+ if currentMinDeltaV < minDeltaV
+ minDeltaV = currentMinDeltaV;
+ optimalDate = ArrivalDate;
+ optimalDepartureDate = LaunchDate; % Store the departure date corresponding to the optimal TOF
+ optimalTOF = tofHours; % Store the optimal TOF in hours
+ end
+ end
+end
+
+% Display the result
+fprintf('Optimal Departure Date: %s\nOptimal Arrival Date: %s\nMinimum Delta-V: %.5f km/s\nOptimal TOF: %d hours\n', optimalDepartureDate, optimalDate, minDeltaV, optimalTOF);
\ No newline at end of file