From 02906c76ff72e69c09ed49bd766b24f967b4b0f4 Mon Sep 17 00:00:00 2001 From: Thomas West <tw00956@surrey.ac.uk> Date: Tue, 3 Sep 2024 15:11:56 +0000 Subject: [PATCH] Animation for the transfer trajectory between Earth and Mars --- Animation_Trajectory_Earth_Mars.m | 238 ++++++++++++++++++++++++++++++ 1 file changed, 238 insertions(+) create mode 100644 Animation_Trajectory_Earth_Mars.m diff --git a/Animation_Trajectory_Earth_Mars.m b/Animation_Trajectory_Earth_Mars.m new file mode 100644 index 0000000..d8023cc --- /dev/null +++ b/Animation_Trajectory_Earth_Mars.m @@ -0,0 +1,238 @@ +%% Script Description (Earth To Mars Trajectory Animation) +% This MATLAB script simulates the trajectory of a spacecraft during an interplanetary transfer from Earth to Mars. +% The script utilizes the SPICE toolkit to access precise planetary ephemeris data, computes necessary orbital +% parameters using the Lambert solver, and visualizes the resulting transfer trajectory in a dynamic 3D animation. +% +% Constants +% Establishes conversion factors for units and defines necessary orbital elements and gravitational parameters. +% These constants are crucial for the calculations that determine the trajectory. +% +% Calling SPICE Functions +% Sets up the SPICE toolkit paths and loads planetary data kernels. These functions allow the script to retrieve +% accurate positional data for Earth and Mars for specific time points, which are essential for trajectory calculations. +% +% Date Calculations +% Computes the precise launch and arrival dates based on a specified mission duration. This section is vital for planning +% the mission timeline and ensuring the calculations are based on accurate temporal data. +% +% Lambert Solver +% Implements Lambert's problem to find the optimal velocities at departure and arrival. This part is crucial for +% determining the most efficient path (in terms of delta-v) for the spacecraft to follow. +% +% Propagate Orbit +% Uses numerical methods (ode45) to propagate the spacecraft's trajectory from Earth to Mars. This part simulates +% the spacecraft’s path based on the initial conditions derived from the Lambert solution. +% +% Animation and Visualization +% Dynamically plots the orbits of Earth and Mars, and the spacecraft's trajectory in 3D space. This visualization includes: +% - Real-time updating of Earth and Mars positions. +% - Spacecraft trajectory shown with progressing animation to reflect its journey. +% - Annotations and markers to clearly identify the celestial bodies and spacecraft at any given simulation time. +% This section helps in visual assessment and presentation of the transfer trajectory, enhancing understanding of the +% mission dynamics. +% +% 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 = '01 Nov 2026 12:30:00 (UTC)'; % Define the launch date in UTC +ArrivalDate = '07 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 / 86400; % 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 ephemeris time in the vector + Xmars(:, tt) = cspice_spkezr('4', et(tt), 'ECLIPJ2000', 'NONE', '10'); % Retrieve the state vector of Mars at time tt + Xearth(:, tt) = cspice_spkezr('399', et(tt), 'ECLIPJ2000', 'NONE', '10'); % Retrieve the state vector of Earth at time tt + % Convert from km to AU + Xmars_AU(:, tt) = Xmars(:, tt) * kmToAU; % Convert Mars' distance from km to AU + Xearth_AU(:, tt) = Xearth(:, tt) * kmToAU; % Convert Earth's distance from km to AU +end + +%% Lambert Solver +delta_t = (et_end - et0); % Calculate the time difference between the start and end ephemeris times +epsilon = 1e-6; % Set a small epsilon value for numerical stability or small error margin in calculations + +% Spacecraft initial position vector +r0 = Xearth(1:3, 1); % Extract the initial position vector of the spacecraft relative to Earth + +% Spacecraft final position vector at Mars +rf = Xmars(1:3, end); % Extract the final position vector of the spacecraft relative to Mars + +% Call the Lambert solver function with the gathered data +[v0, vf] = lambert_solver(r0, rf, delta_t, mu, DM); + +% Output the initial and final velocities in km/s +disp('Initial velocity vector (km/s):'); % Display message indicating that the next line will show the initial velocity vector in km/s. +disp(v0); % Output the initial velocity vector of the spacecraft. +disp('Final velocity vector (km/s):'); % Display message indicating that the next line will show the final velocity vector in km/s. +disp(vf); % Output the final velocity vector of the spacecraft. + +%% Propagate Orbit +% Define the time span of the simulation +tspan = [0 delta_t]; % Set the time span for the orbit propagation from 0 to the total duration in seconds. + +% Initial state vector (position and velocity) +initial_state = [r0; v0]; % Create an initial state vector that includes both position and velocity of the spacecraft. + +% Propagate the trajectory using ode45 +options = odeset('RelTol', 1e-33, 'AbsTol', 1e-33); % Set options for the ODE solver to use very small tolerance values for high accuracy. +[~, ~] = ode45(@(t, y) Orbit_Propagation(t, y, mu), tspan, initial_state, options); % Use the ODE45 solver to simulate the orbit, ignoring the usual output arguments. + +% Calculate spacecraft trajectory using ode45 +tspan = [0 total_duration]; % Define the time span from 0 to total mission duration in seconds. +initial_state = [r0; v0]; % Define the initial state vector combining the position and velocity vectors. +options = odeset('RelTol', 1e-13, 'AbsTol', 1e-13); % Set the solver options for precision. +[t_out, state_out] = ode45(@(t, y) Orbit_Propagation(t, y, mu), tspan, initial_state, options); % Run the ODE solver with the defined function, time span, initial state, and options. +positions = state_out(:, 1:3); % Extract the position vectors from the solution. +positions_AU = positions * kmToAU; % Convert the position vectors from kilometers to Astronomical Units. + +% Normalize spacecraft trajectory time to match et +spacecraft_time_normalized = linspace(0, 1, length(t_out)); % Normalize the time output of the spacecraft trajectory. +et_normalized = linspace(0, 1, length(et)); % Normalize the ephemeris time vector to the same scale as the spacecraft time. + +%% Animation +maxDataPoint = max(abs(positions_AU), [], 'all'); % Find the maximum absolute value from all data points in the positions matrix for normalization purposes +buffer = 0.37 * maxDataPoint; % Calculate a buffer (37% of the maximum data point) to ensure all data points are comfortably within the plot boundaries + +% Initialize figure and static elements +figure; % Create a new figure for plotting +hold on; % Retain all plots in the figure +grid on; % Enable grid for better visualization of space +axis equal; % Ensure that the scale of the x, y, and z axes are the same + +% Set dynamic limits based on the data +xlim([-maxDataPoint-buffer, maxDataPoint+buffer]); % Set the x-axis limits based on the maximum data point plus buffer +ylim([-maxDataPoint-buffer, maxDataPoint+buffer]); % Set the y-axis limits similarly +zlim([-maxDataPoint-buffer, maxDataPoint+buffer]); % Set the z-axis limits similarly + +% 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 + +xlabel('X Distance (AU)'); % Label the x-axis +ylabel('Y Distance (AU)'); % Label the y-axis +zlabel('Z Distance (AU)'); % Label the z-axis +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 + +% Plot Sun +sun = plot3(0, 0, 0, 'yo', 'MarkerSize', 15, 'MarkerFaceColor', 'yellow'); % Plot the sun as a yellow circle at the origin + +% Initialize plot objects for Earth, Mars, and the Spacecraft trajectory +h_earth = plot3(NaN, NaN, NaN, 'b-', 'LineWidth', 1, 'MarkerFaceColor', 'blue'); % Initialize the Earth plot with no data +h_mars = plot3(NaN, NaN, NaN, 'r-', 'LineWidth', 1, 'MarkerFaceColor', 'red'); % Initialize the Mars plot with no data +h_spacecraft = plot3(NaN, NaN, NaN, 'm-', 'LineWidth', 1); % Initialize the spacecraft trajectory plot with no data + +% Dynamic markers that move along the orbits +earth_marker = plot3(NaN, NaN, NaN, 'bo', 'MarkerSize', 5, 'MarkerFaceColor', 'blue'); % Marker for Earth +mars_marker = plot3(NaN, NaN, NaN, 'ro', 'MarkerSize', 5, 'MarkerFaceColor', 'red'); % Marker for Mars +spacecraft_marker = plot3(NaN, NaN, NaN, 'm^', 'MarkerSize', 5, 'MarkerFaceColor', 'magenta'); % Marker for the spacecraft + +% Add legends for the markers after initializing all plot objects +legend([h_earth, earth_marker, h_mars, mars_marker, h_spacecraft, spacecraft_marker, sun], ... + {'Earth Orbit', 'Earth', 'Mars Orbit', 'Mars', 'Spacecraft Trajectory', 'Spacecraft', 'Sun'}, ... + 'Location', 'bestoutside'); % Create a legend that clearly identifies each plot element + +% Set desired font size for the date legend +dateLegendFontSize = 13; + +% Create a text object for the date legend at the top right of the plot +dateLegend = text(0.99, 0.99, 'Date: ', 'FontSize', dateLegendFontSize, ... + 'HorizontalAlignment', 'right', 'VerticalAlignment', 'top', ... + 'Units', 'normalized', 'Color', 'k', 'FontWeight', 'bold', 'BackgroundColor', 'white'); % Place a date legend in the top right corner + +% Animation loop +for tt = 1:length(et) + % Update Earth and Mars positions + Xearth = cspice_spkezr('399', et(tt), 'ECLIPJ2000', 'NONE', '10'); % Retrieve Earth's state vector at time tt + Xmars = cspice_spkezr('4', et(tt), 'ECLIPJ2000', 'NONE', '10'); % Retrieve Mars' state vector at time tt + Xearth_AU = Xearth(1:3) * kmToAU; % Convert Earth's position from km to AU + Xmars_AU = Xmars(1:3) * kmToAU; % Convert Mars' position from km to AU + + % Append new positions to Earth and Mars orbit traces + set(h_earth, 'XData', [get(h_earth, 'XData') Xearth_AU(1)], 'YData', [get(h_earth, 'YData') Xearth_AU(2)], 'ZData', [get(h_earth, 'ZData') Xearth_AU(3)]); + set(h_mars, 'XData', [get(h_mars, 'XData') Xmars_AU(1)], 'YData', [get(h_mars, 'YData') Xmars_AU(2)], 'ZData', [get(h_mars, 'ZData') Xmars_AU(3)]); + + % Update marker positions to the latest Earth and Mars positions + set(earth_marker, 'XData', Xearth_AU(1), 'YData', Xearth_AU(2), 'ZData', Xearth_AU(3)); + set(mars_marker, 'XData', Xmars_AU(1), 'YData', Xmars_AU(2), 'ZData', Xmars_AU(3)); + + % Spacecraft trajectory data + [~, closestIndex] = min(abs(spacecraft_time_normalized - et_normalized(tt))); % Find the index of the closest normalized time for the spacecraft's position + if closestIndex <= length(positions_AU) + set(h_spacecraft, 'XData', positions_AU(1:closestIndex, 1), 'YData', positions_AU(1:closestIndex, 2), 'ZData', positions_AU(1:closestIndex, 3)); + set(spacecraft_marker, 'XData', positions_AU(closestIndex, 1), 'YData', positions_AU(closestIndex, 2), 'ZData', positions_AU(closestIndex, 3)); + end + + % Update the date display + current_date = cspice_et2utc(et(tt), 'C', 0); % Convert the current ephemeris time to a human-readable UTC date + set(dateLegend, 'String', sprintf('Date: %s', current_date)); % Update the date legend with the current date + + drawnow; % Update figure window with the latest changes + pause(0.0001); % Short pause to slow down the animation slightly +end + +hold off; + +%% Clean up by unloading SPICE kernels +cspice_kclear; % Unload all SPICE kernels and clear the SPICE subsystem \ No newline at end of file -- GitLab