diff --git a/Trajectory_High_Inclination_To_SOI.m b/Trajectory_High_Inclination_To_SOI.m new file mode 100644 index 0000000000000000000000000000000000000000..c8bd210461221765890a452d0e7be977fc42752a --- /dev/null +++ b/Trajectory_High_Inclination_To_SOI.m @@ -0,0 +1,234 @@ +%% Function Description (Orbit Propagation from Earth Parking Orbit to Leaving Earth's SOI) +% This script propagates the orbit of a spacecraft from an initial Earth parking orbit +% to the point of leaving Earth's Sphere of Influence (SOI). The goal is to optimize the +% spacecraft's trajectory to minimize the required delta-V (change in velocity) while +% achieving a successful escape from Earth's gravitational influence. The script considers +% orbital mechanics, spacecraft dynamics, and gravitational interactions with Earth and the Moon. +% +% The script starts by clearing the workspace, closing figures, and setting the format for +% numerical output to ensure a fresh and clean environment. It then defines constants, including +% conversion factors and gravitational parameters, to set up the simulation. +% +% The key steps and calculations in this script are: +% - Defining the orbital elements of the initial Earth parking orbit, such as semi-major axis, +% eccentricity, inclination, right ascension of the ascending node (RAAN), and argument of periapsis. +% - Loading SPICE kernels to retrieve accurate planetary and lunar ephemerides, gravitational +% constants, and other necessary data for precise orbit calculations. +% - Calculating state vectors (position and velocity) for Earth and the Moon over the defined +% observation period using SPICE functions. These vectors are used to determine the spacecraft's +% trajectory relative to the Earth and Moon. +% - Propagating the spacecraft's orbit over one orbital period using Kepler's third law and calculating +% the state vectors across the time vector. The state vectors provide the position and velocity of +% the spacecraft at different points along its orbit. +% - Utilizing a Lambert solver to calculate the initial and final velocity vectors required for the +% optimal transfer from Earth to the desired position near Mars. The Lambert solver accounts for the +% gravitational influence of the Sun and the specific orbital characteristics of the spacecraft. +% - Iterating over a range of possible times of flight (TOFs) to find the optimal TOF that minimizes the +% delta-V required to achieve the desired transfer. The script evaluates each TOF and selects the one +% that results in the lowest delta-V, indicating the most efficient trajectory. +% - Outputting the optimal initial and final velocities, the best TOF, the minimum delta-V, and the precise +% arrival date at the destination. These outputs are essential for planning the spacecraft's maneuvers +% and ensuring a successful mission trajectory. +% +% The script is designed to be flexible and allows for adjustments to various parameters, such as orbital +% elements, gravitational constants, and observation periods. It provides a comprehensive approach to +% optimizing spacecraft trajectories for missions involving Earth escape and interplanetary transfers. +% +% 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(313.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 = '3 Nov 2026 16:00:00 (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 + +moon_pos = Xmoon(1:3, :); + +Earth_helio_pos = Xearth(1:3, 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.795; 2.364; 0.612]; +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.0879; 720190.384; 186378.606]; % Extract Mars' position vector at the arrival time + +%% Lambert Solver +% Define the range of TOFs +tof_days = 3:1/86400:4; % 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); + +% Define the start date and the best TOF found +StartDate = '1 Nov 2026 20:04:03 (UTC)'; % Update this line with your desired start date +best_tof_days = best_tof; % Example best TOF in days (replace with the actual best TOF found) + +% Remove the time zone for parsing +StartDateClean = extractBefore(StartDate, ' (UTC)'); + +% Convert the cleaned start date to datetime +start_date_dt = datetime(StartDateClean, 'InputFormat', 'd MMM yyyy HH:mm:ss', 'TimeZone', 'UTC'); + +% Calculate the arrival date by adding the best TOF (in days) to the start date +arrival_date_dt = start_date_dt + days(best_tof_days); + +% Display the calculated arrival date +fprintf('The precise arrival date is: %s\n', datestr(arrival_date_dt, 'dd-mmm-yyyy HH:MM:SS (UTC)')); + +% Optional: Convert the arrival date to ephemeris time using SPICE (if needed) +arrival_et = cspice_str2et(datestr(arrival_date_dt, 'dd-mmm-yyyy HH:MM:SS (UTC)')); + +% Display the ephemeris time (if needed) +fprintf('Ephemeris time of arrival: %.2f seconds past J2000\n', arrival_et); \ No newline at end of file