From 0d6f8ce135307e94cad746f24d37e19031b5200e Mon Sep 17 00:00:00 2001
From: Thomas West <tw00956@surrey.ac.uk>
Date: Tue, 3 Sep 2024 16:24:42 +0000
Subject: [PATCH] Porkchop plots for the 2024-2028 launch window.
---
Porkchop_Plot_Earth_Mars.m | 293 +++++++++++++++++++++++++++++++++++++
1 file changed, 293 insertions(+)
create mode 100644 Porkchop_Plot_Earth_Mars.m
diff --git a/Porkchop_Plot_Earth_Mars.m b/Porkchop_Plot_Earth_Mars.m
new file mode 100644
index 0000000..eeb39d5
--- /dev/null
+++ b/Porkchop_Plot_Earth_Mars.m
@@ -0,0 +1,293 @@
+%% Script Description (Porkchop Plot)
+% This MATLAB script is designed to plan interplanetary missions between Earth and Mars by calculating
+% potential launch windows and optimal transfer trajectories using the porkchop plot method. The script
+% leverages NASA's SPICE toolkit to access precise planetary ephemeris data and computes the delta-v
+% requirements for various combinations of launch dates and times of flight (TOF).
+%
+% Constants
+% Defines necessary constants such as unit conversions, orbital elements, and the gravitational parameter
+% of the Sun, which are crucial for the trajectory calculations.
+%
+% Calling SPICE Functions
+% Sets up necessary directories for the SPICE toolkit and loads essential planetary data kernels. This
+% functionality is vital for retrieving accurate position and velocity data for Earth and Mars.
+%
+% Launch and Porkchop Plot Details
+% Specifies the ideal launch windows identified from previous analyses and sets the range for exploring
+% potential launch dates and TOF. This section is fundamental in determining the scope of the mission planning.
+%
+% Data Collection
+% Collects state vectors for Earth and Mars over a range of dates and TOFs. This data collection is crucial
+% for calculating the delta-v requirements for different mission scenarios.
+%
+% Orbit Calculations
+% Processes the collected data to calculate orbital parameters and prepares for delta-v calculations. This
+% step includes initializing arrays for storing results and setting up loops for data processing.
+%
+% Lambert Solver for Short and Long Way
+% Implements the Lambert solver to calculate the necessary velocities for both the short and long transfer
+% paths. This calculation determines the initial and final velocities required for the spacecraft to
+% transition between Earth and Mars orbits.
+%
+% Calculate and Filter Delta-V for Both Short and Long Ways
+% Uses the results from the Lambert solver to calculate the total delta-v for each combination of departure
+% date and TOF. Filters out combinations that exceed a specified maximum delta-v limit, focusing on feasible
+% mission profiles.
+%
+% Plotting the Pork Chop Plot
+% Visualizes the results using a contour plot, commonly known as a porkchop plot. This plot illustrates the
+% delta-v requirements across different launch dates and TOFs, helping identify the most efficient windows
+% for mission launch. Highlights include:
+% - Display of delta-v values using color contours.
+% - Identification and annotation of the minimum delta-v points for both short and long transfer paths.
+% - Interactive features that allow users to explore different data points on the plot.
+%
+% This code is specifically focused on identifying ideal launch dates for Earth-to-Mars missions within the years 2024-2028.
+%
+% 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;
+% Set the text size format
+set(groot, 'DefaultAxesFontSize', 18, 'DefaultTextFontSize', 18, 'DefaultAxesFontName', 'Times New Roman', 'DefaultTextFontName', 'Times New Roman');
+
+%% 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_sun = 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 and Porkchop Plot Details
+% For 2024-2028 Launch Window there are two ideal launch dates:
+% 1. Closest approach date and time: 2025 JAN 12 13:37:09
+% 2. Closest approach date and time: 2027 FEB 20 00:12:53
+Ideal_Departure_1 = '12 Jan 2025 12:37:09 (UTC)'; % Define the first ideal departure date and time
+Ideal_Departure_2 = '20 Feb 2027 00:12:53 (UTC)'; % Define the second ideal departure date and time
+
+% Select your ideal departure date:
+Ideal_Departure = Ideal_Departure_2; % Choose the first ideal departure date for further calculations
+
+% First and Last possible departure dates:
+Before_ID = -335; % Days before the ideal departure date to start the departure window
+After_ID = 396; % Days after the ideal departure date to end the departure window
+dLD = 1; % Day increments for each step in the departure window
+
+% Time of flight (TOF)
+TOF_min = 50; % Minimum time of flight in days
+TOF_max = 900; % Maximum time of flight in days
+dTOF = 1; % Increment in days for time of flight
+
+Max_v = 20; % Maximum allowable velocity in km/s for spacecraft
+
+% Convert ideal departure date to ephemeris time (ET)
+Ideal_Departure_et = cspice_str2et(Ideal_Departure); % Convert the selected ideal departure date to ephemeris time
+et0 = Ideal_Departure_et + Before_ID * 86400; % Calculate the start of the launch window in ephemeris time
+etf = Ideal_Departure_et + After_ID * 86400; % Calculate the end of the launch window in ephemeris time
+num_departure_days = (After_ID - Before_ID) / dLD + 1; % Calculate the total number of days in the departure window
+
+%% Data Collection
+num_TOF_days = TOF_max - TOF_min + 1; % Calculate the total number of TOF days to analyze
+
+% Initialize matrices to store the state vectors of Mars and Earth
+Xmars = zeros(6, num_departure_days, num_TOF_days); % Create a 3D array to store Mars' state vector for each TOF and each launch date
+XEarth = zeros(6, num_departure_days); % Create a 2D array to store Earth's state vector for each launch date
+
+% Loop over each departure day
+for d = 1:num_departure_days
+ current_et = et0 + (d - 1) * 86400; % Calculate the current ephemeris time for the departure day, incrementing by one day in seconds
+ XEarth(:, d) = cspice_spkezr('399', current_et, 'ECLIPJ2000', 'NONE', '10'); % Retrieve and store Earth's state vector from SPICE at current_et
+
+ % Nested loop for each TOF from minimum to maximum TOF
+ for tof = TOF_min:dTOF:TOF_max
+ mars_et = current_et + tof * 86400; % Calculate the ephemeris time for Mars' position at the end of the current TOF
+ Xmars(:, d, tof - TOF_min + 1) = cspice_spkezr('4', mars_et, 'ECLIPJ2000', 'NONE', '10'); % Retrieve and store Mars' state vector at mars_et
+ end
+end
+
+% Display or analyze the collected data
+disp('Data collection complete.'); % Notify the user that data collection is finished
+% Further processing can include calculations for the pork chop plot using the collected data.
+
+%% Programming Calculations
+total_iterations = num_departure_days * num_TOF_days; % Calculate the total number of iterations needed by multiplying the number of departure days by the number of TOF days
+iteration_count = 0; % Initialize a counter to keep track of iterations
+
+C3_shorts = inf(num_departure_days, num_TOF_days); % Create a matrix to store characteristic energy (C3) values, initialized with 'inf' (infinity)
+v_inf_shorts = inf(num_departure_days, num_TOF_days); % Create a matrix to store excess velocity (v-infinity) values, initialized with 'inf' (infinity)
+
+% Display the Number of Combinations
+num_departure_days = floor((etf - et0) / 86400) + 1; % Recalculate the number of departure days from the range of ephemeris time (ET) by dividing the total seconds by the number of seconds per day and adding one to include both endpoints
+num_TOF_days = TOF_max - TOF_min + 1; % Recalculate the total number of time-of-flight days by subtracting the minimum TOF from the maximum TOF and adding one to include both endpoints
+total_combinations = num_departure_days * num_TOF_days; % Calculate the total number of combinations possible between departure days and TOF days
+
+fprintf('Total combinations of departure and arrival days: %d\n', total_combinations); % Print the total number of combinations to the console
+
+%% Lambert Solver for Short and Long Way
+% Initialize four matrices to store initial and final velocities for both short and long transfer trajectories.
+V0_short = zeros(3, num_departure_days, num_TOF_days); % Initial velocity for short-way transfers
+Vf_short = zeros(3, num_departure_days, num_TOF_days); % Final velocity for short-way transfers
+V0_long = zeros(3, num_departure_days, num_TOF_days); % Initial velocity for long-way transfers
+Vf_long = zeros(3, num_departure_days, num_TOF_days); % Final velocity for long-way transfers
+
+% Loop over each departure day
+for d = 1:num_departure_days
+ % Nested loop for each time of flight (TOF) index
+ for tof_index = 1:num_TOF_days
+ r0 = XEarth(1:3, d); % Extract the initial position vector of Earth for the current departure day
+ rf = Xmars(1:3, d, tof_index); % Extract the final position vector of Mars for the current TOF
+ Dt = (TOF_min + tof_index - 1) * 86400; % Calculate the duration in seconds for the current TOF
+
+ % Calculate short-way trajectory using Lambert's solver
+ [v0_short, vf_short, exitflag_short] = lambert_solver(r0, rf, Dt, mu_sun, +1);
+ if isempty(v0_short) || isempty(vf_short) % Check if the Lambert solver failed to find a solution
+ fprintf('Short way solution not found for d = %d, TOF = %d\n', d, Dt/86400); % Log failure message for short-way
+ continue; % Skip to the next iteration if no solution is found
+ end
+
+ % Calculate long-way trajectory using Lambert's solver
+ [v0_long, vf_long, exitflag_long] = lambert_solver(r0, rf, Dt, mu_sun, -1);
+ if isempty(v0_long) || isempty(vf_long) % Check if the Lambert solver failed to find a solution
+ fprintf('Long way solution not found for d = %d, TOF = %d\n', d, Dt/86400); % Log failure message for long-way
+ continue; % Skip to the next iteration if no solution is found
+ end
+
+ % Store the calculated velocities in their respective matrices
+ V0_short(:, d, tof_index) = v0_short; % Store the initial velocity for the short-way
+ Vf_short(:, d, tof_index) = vf_short; % Store the final velocity for the short-way
+ V0_long(:, d, tof_index) = v0_long; % Store the initial velocity for the long-way
+ Vf_long(:, d, tof_index) = vf_long; % Store the final velocity for the long-way
+ end
+end
+
+%% Calculate and Filter Delta-V for Both Short and Long Ways
+delta_v_short = calculate_delta_v(V0_short, Vf_short, XEarth, Xmars, Max_v, num_departure_days, num_TOF_days, TOF_min, dTOF); % Calculate delta-v values for the short-way transfer using the 'calculate_delta_v' function
+delta_v_long = calculate_delta_v(V0_long, Vf_long, XEarth, Xmars, Max_v, num_departure_days, num_TOF_days, TOF_min, dTOF); % Calculate delta-v values for the long-way transfer using the 'calculate_delta_v' function
+
+%% Date Formatting
+% Base date for J2000
+j2000 = datetime(2000, 1, 1, 12, 0, 0); % Define the exact J2000 epoch, January 1, 2000, at noon
+
+% Conversion from ET to datetime, considering et0 and etf as SPICE ET
+departure_dates_vec = j2000 + seconds(et0:86400*dLD:etf); % Create a vector of datetime objects representing each departure day by adding seconds to the J2000 base date
+
+% Convert datetime to datenum for contourf plotting
+departure_dates_num = datenum(departure_dates_vec); % Convert the vector of datetime objects into MATLAB datenum format for use in plotting functions like contourf
+
+%% Plotting the Pork Chop Plot
+figure; % Create a new figure window for the plot
+hold on; % Retain plots so that new plots added to the figure do not delete existing plots
+
+% Define TOF range and number of levels for contour plotting
+TOF_range = TOF_min:dTOF:(TOF_min + num_TOF_days - 1); % Create a range of Time of Flight (TOF) values from minimum to maximum
+% Define levels from 5 to Max_v, each corresponding to an integer
+levels = 5:0.5:Max_v; % 1 km/s steps from 5 to Max_v
+Num_Levels = numel(levels) - 1; % Number of transitions
+
+% Plot data with specified levels
+[~, hShort] = contourf(departure_dates_num, TOF_range, delta_v_short', [levels inf], 'LineColor', 'none');
+[C, hLong] = contourf(departure_dates_num, TOF_range, delta_v_long', [levels inf], 'LineColor', 'none');
+
+% Set the colormap to 'jet' with one color per level
+colormap(jet(Num_Levels));
+
+% Customize colorbar
+hColorbar = colorbar; % Add a colorbar to the figure
+caxis([5 Max_v]); % Set the color axis scaling to range from 5 to Max_v
+set(hColorbar, 'Ticks', 5:1:Max_v); % Set ticks at each level from 5 to Max_v
+set(hColorbar, 'TickLabels', arrayfun(@(x) sprintf('%.1f km/s', x), 5:1:Max_v, 'UniformOutput', false)); % Label each tick with velocity
+
+% Ensure colors map directly to ticks by using 'edges' of each color
+set(hColorbar, 'Limits', [5 Max_v]); % Ensures limits are exactly from 5 to Max_v
+% Find minimum delta-V points for both short and long ways
+[minValShort, idxShort] = min(delta_v_short(:)); % Find the minimum delta-V value and index for the short way
+[minValLong, idxLong] = min(delta_v_long(:)); % Find the minimum delta-V value and index for the long way
+[rowShort, colShort] = ind2sub(size(delta_v_short), idxShort); % Convert linear index to subscript indices for short way
+[rowLong, colLong] = ind2sub(size(delta_v_long), idxLong); % Convert linear index to subscript indices for long way
+
+minDateShort = departure_dates_vec(rowShort); % Determine the departure date for the minimum delta-V short way
+minTOFShort = TOF_range(colShort); % Determine the TOF for the minimum delta-V short way
+minDateLong = departure_dates_vec(rowLong); % Determine the departure date for the minimum delta-V long way
+minTOFLong = TOF_range(colLong); % Determine the TOF for the minimum delta-V long way
+arrivalDateShort = minDateShort + days(minTOFShort); % Calculate the arrival date based on the departure date and TOF for short way
+arrivalDateLong = minDateLong + days(minTOFLong); % Calculate the arrival date based on the departure date and TOF for long way
+
+% Plot minimum points
+hPointShort = plot(datenum(minDateShort), minTOFShort, 'm^', 'MarkerSize', 4, 'MarkerFaceColor', 'm'); % Mark the minimum delta-V point for the short way with a magenta triangle
+hPointLong = plot(datenum(minDateLong), minTOFLong, 'mv', 'MarkerSize', 4, 'MarkerFaceColor', 'm'); % Mark the minimum delta-V point for the long way with a magenta inverted triangle
+
+% Convert ideal departure date to ephemeris time (ET)
+Ideal_Departure_et = cspice_str2et(Ideal_Departure); % Convert the selected ideal departure date to ephemeris time
+et0 = Ideal_Departure_et + Before_ID * 86400; % Adjust the start of the time window by converting days to seconds
+etf = Ideal_Departure_et + After_ID * 86400; % Adjust the end of the time window by converting days to seconds
+num_departure_days = (After_ID - Before_ID) / dLD + 1; % Calculate the total number of departure days within the window
+
+% Convert ET dates back to datetime for year extraction
+departure_date_0 = datetime(cspice_et2utc(et0, 'C', 0), 'InputFormat', 'yyyy MMM dd HH:mm:ss', 'TimeZone', 'UTC'); % Convert start ET to datetime
+departure_date_f = datetime(cspice_et2utc(etf, 'C', 0), 'InputFormat', 'yyyy MMM dd HH:mm:ss', 'TimeZone', 'UTC'); % Convert end ET to datetime
+
+% Extract years from the calculated datetime
+min_year = year(departure_date_0); % Extract the year from the start date
+max_year = year(departure_date_f); % Extract the year from the end date
+
+% Add label to the colorbar
+ylabel(hColorbar, '∆V (km/s)'); % Label the colorbar to indicate it represents delta-V in km/s
+
+% Labels, title, and legend
+xlabel('Launch Date'); % Label the x-axis with 'Launch Date'
+ylabel('Time of Flight (days)'); % Label the y-axis with 'Time of Flight (days)'
+title_str = sprintf('%d-%d Earth-to-Mars Transfer Porkchop Plot with ∆V Contours', min_year, max_year); % Create a title string that includes the year range
+title(title_str); % Set the title of the plot
+
+% Add a legend to explain plot elements
+legend([hShort, hPointShort, hPointLong], ...
+ 'Total ∆V', ...
+ sprintf('Min ∆V Short Way: Departure: %s, Arrival: %s, Time of Flight: %g days, ∆V: %.5g km/s', datestr(minDateShort, 'dd-mmm-yyyy'), datestr(arrivalDateShort, 'dd-mmm-yyyy'), minTOFShort, minValShort), ... % Detailed legend entry for the minimum delta-V short way
+ sprintf('Min ∆V Long Way: Departure: %s, Arrival: %s, Time of Flight: %g days, ∆V: %.5g km/s', datestr(minDateLong, 'dd-mmm-yyyy'), datestr(arrivalDateLong, 'dd-mmm-yyyy'), minTOFLong, minValLong), ... % Detailed legend entry for the minimum delta-V long way
+ 'Location', 'northwest'); % Specify the location of the legend
+
+% Date formatting and axis settings
+datetick('x','dd-mmm-yyyy') % Set the x-axis ticks to display dates in 'day-month-year' format
+axis tight % Tighten the axis limits to the data
+grid on; % Enable the grid for better readability
+ax = gca; % Get the current axes object
+set(ax, 'Layer', 'top', 'GridLineStyle', '-', 'GridAlpha', 0.3, 'GridColor', 'k', 'LineWidth', 0.2); % Set grid properties and layer order
+
+% Set ticks on the x-axis
+ax.XTick = linspace(departure_dates_num(1), departure_dates_num(end), 12); % Distribute x-axis ticks evenly across the displayed date range
+datetick('x', 'dd-mmm-yyyy', 'keeplimits', 'keepticks'); % Maintain the date ticks formatting even when resizing the plot
+
+% Attach the data cursor update function
+dcm_obj = datacursormode(gcf); % Enable data cursor mode on the current figure
+set(dcm_obj, 'UpdateFcn', @Cursor_Function); % Set the function that updates the text when the data cursor is used
+
+hold off; % Release the plot hold
+
+%% Clean up by unloading SPICE kernels
+cspice_kclear; % Unload all SPICE kernels and clear the SPICE subsystem
\ No newline at end of file
--
GitLab