Animation of the orbits of Earth and Mars

Animation_Earth_Mars.m

+%% Script Description (Earth Mars Animation)
+% This MATLAB script simulates and visualizes the orbital trajectories of Earth and Mars over a specified period. 
+% It uses the SPICE toolkit to retrieve precise ephemeris data, calculating and animating the planetary positions 
+% and velocities in three-dimensional space. The script highlights the dynamics of Earth and Mars as they orbit 
+% around the Sun, providing visual insights into their relative motions.
+%
+% Constants
+% Defines the necessary conversion factors for distances (kilometers to Astronomical Units) and angles (degrees to radians), 
+% and sets up essential orbital and gravitational parameters for the simulation.
+%
+% Calling SPICE Functions
+% Sets up paths for SPICE toolkits and loads ephemeris data which are crucial for retrieving accurate celestial mechanics 
+% data needed for the simulation.
+%
+% Earth Mars State Vectors
+% Calculates the state vectors (position and velocity) for Earth and Mars from a specified start date to an end date, converting
+% these vectors into Astronomical Units for visualization. This section forms the basis for plotting the trajectories.
+%
+% Orbital Velocity
+% Extracts the velocity components from the state vectors and calculates the magnitude of these velocities. This step
+% is vital for understanding the speeds at which the planets are moving in their orbits.
+%
+% Animation
+% Initializes a 3D plot and systematically updates it to create an animation that shows the orbits of Earth and Mars over the specified period. 
+% Features include:
+% - Dynamic updates of planetary positions with corresponding velocity vectors.
+% - A running date display that shows the simulation's current date.
+% - Plotting the Sun as a reference point in the center of the coordinate system.
+% - Use of different markers and colors to distinctly identify each celestial body and their paths.
+%
+% This script is useful for educational purposes, providing a visual demonstration of how Earth and Mars move relative to each other
+% and the Sun over time, specifically focused on the ideal observation or mission planning dates between 2024 and 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;
+
+%% 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; % 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 Mars State Vectors
+% Set the start and end dates for the simulation period
+StartDate = '15 Apr 2024 12:30:00 (UTC)';
+EndDate = '15 Dec 2028 12:30:00 (UTC)'; % Update this line with your desired end date
+
+% Convert both dates to Ephemeris Time using the SPICE toolkit
+et0 = cspice_str2et(StartDate);  % Convert start date to seconds past J2000
+et_end = cspice_str2et(EndDate);  % Convert end date to seconds past J2000
+
+% Calculate the total duration in seconds between start and end dates
+total_duration = et_end - et0;  % Time difference in seconds
+
+% Calculate the number of days between the start and end date
+num_days = total_duration / 86400; % Convert seconds to days by dividing by the number of seconds in a day
+
+% Round num_days to the nearest whole number to avoid fractional days
+num_days = round(num_days);  % Ensure the number of days is an integer
+
+% Generate a time vector from start to end date using the calculated number of days
+et = linspace(et0, et_end, num_days + 1); % Create linearly spaced vector including both endpoints
+
+% Retrieve state vectors for each time point
+for tt = 1:length(et)  % Loop through each time point
+    Xmars(:, tt) = cspice_spkezr('4', et(tt), 'ECLIPJ2000', 'NONE', '10');  % Get Mars (4) state vector at time tt
+    Xearth(:, tt) = cspice_spkezr('399', et(tt), 'ECLIPJ2000', 'NONE', '10');  % Get Earth (339) state vector at time tt
+    % Convert from km to AU
+    Xmars_AU(:, tt) = Xmars(:, tt) * kmToAU;  % Convert Mars' position from km to AU
+    Xearth_AU(:, tt) = Xearth(:, tt) * kmToAU;  % Convert Earth's position from km to AU
+end
+
+%% Orbital Velocity
+% Velocity components are typically the last three elements of the state vector
+v_earth_x = Xearth(4, :); % Extract the X velocity component of Earth from the state vector
+v_earth_y = Xearth(5, :); % Extract the Y velocity component of Earth from the state vector
+v_earth_z = Xearth(6, :); % Extract the Z velocity component of Earth from the state vector
+
+v_mars_x = Xmars(4, :); % Extract the X velocity component of Mars from the state vector
+v_mars_y = Xmars(5, :); % Extract the Y velocity component of Mars from the state vector
+v_mars_z = Xmars(6, :); % Extract the Z velocity component of Mars from the state vector
+
+% Calculating the magnitude of the velocity vector (orbital velocity)
+v_earth = sqrt(v_earth_x.^2 + v_earth_y.^2 + v_earth_z.^2); % Compute the magnitude of Earth's velocity vector
+v_mars = sqrt(v_mars_x.^2 + v_mars_y.^2 + v_mars_z.^2); % Compute the magnitude of Mars' velocity vector
+
+%% Animation
+% Initialize static title and figure
+figure;  % Create a new figure window
+xlabel('X Distance (AU)');  % Label the x-axis
+ylabel('Y Distance (AU)');  % Label the y-axis
+zlabel('Z Distance (AU)');  % Label the z-axis
+axis equal;  % Set equal scaling on all axes
+xlim([-2, 2]);  % Set x-axis limits
+ylim([-2, 2]);  % Set y-axis limits
+zlim([-2, 2]);  % Set z-axis limits
+grid on;  % Turn on the grid
+hold on;  % Retain plots so that new plots added to the figure do not delete existing plots
+title(sprintf('3D Orbit of Earth and Mars from %s to %s', StartDate, EndDate));  % Set title with formatted start and end dates
+plot3(0, 0, 0, 'yo', 'MarkerSize', 15, 'MarkerFaceColor', 'yellow'); % Plot the Sun at the origin
+
+% Create a placeholder for the date legend with an initial position
+dateLegendFontSize = 13; % Define font size for the date legend
+dateLegend = text(1, 1, 0, '', 'FontSize', dateLegendFontSize, 'HorizontalAlignment', 'right', 'VerticalAlignment', 'top', 'Units', 'normalized');
+
+% Days per increment variable (e.g., every day)
+daysPerIncrement = 5;  % Set the number of days per increment for the animation steps
+
+% Initialize empty arrays to store the trajectory points
+earthTrajectoryX = [];  % Initialize empty array for Earth's x-coordinates
+earthTrajectoryY = [];  % Initialize empty array for Earth's y-coordinates
+earthTrajectoryZ = [];  % Initialize empty array for Earth's z-coordinates
+marsTrajectoryX = [];   % Initialize empty array for Mars's x-coordinates
+marsTrajectoryY = [];   % Initialize empty array for Mars's y-coordinates
+marsTrajectoryZ = [];   % Initialize empty array for Mars's z-coordinates
+
+% Initialize markers for Earth and Mars at the starting position
+earthMarker = plot3(Xearth_AU(1, 1), Xearth_AU(2, 1), Xearth_AU(3, 1), 'bo', 'MarkerFaceColor', 'blue', 'MarkerSize', 5);  % Place an Earth marker on the plot
+marsMarker = plot3(Xmars_AU(1, 1), Xmars_AU(2, 1), Xmars_AU(3, 1), 'ro', 'MarkerFaceColor', 'red', 'MarkerSize', 5);  % Place a Mars marker on the plot
+
+% Initialize labels for Earth and Mars
+earthLabel = text(Xearth_AU(1, 1)+0.1, Xearth_AU(2, 1)+0.1, Xearth_AU(3, 1), 'Earth\n0 km/s', 'Color', 'blue', 'FontSize', 10, 'HorizontalAlignment', 'left');  % Label for Earth with initial velocity
+marsLabel = text(Xmars_AU(1, 1)+0.1, Xmars_AU(2, 1)+0.1, Xmars_AU(3, 1), 'Mars\n0 km/s', 'Color', 'red', 'FontSize', 10, 'HorizontalAlignment', 'left');  % Label for Mars with initial velocity
+
+% Animation loop for updated legend and labels
+for i = 1:daysPerIncrement:length(et)  % Loop through the ephemeris time array in increments defined by daysPerIncrement
+    % Append current position to trajectory arrays
+    earthTrajectoryX = [earthTrajectoryX Xearth_AU(1, i)];  % Append Earth's current X position
+    earthTrajectoryY = [earthTrajectoryY Xearth_AU(2, i)];  % Append Earth's current Y position
+    earthTrajectoryZ = [earthTrajectoryZ Xearth_AU(3, i)];  % Append Earth's current Z position
+    marsTrajectoryX = [marsTrajectoryX Xmars_AU(1, i)];  % Append Mars's current X position
+    marsTrajectoryY = [marsTrajectoryY Xmars_AU(2, i)];  % Append Mars's current Y position
+    marsTrajectoryZ = [marsTrajectoryZ Xmars_AU(3, i)];  % Append Mars's current Z position
+
+    % Update the plot for trajectories
+    plot3(earthTrajectoryX, earthTrajectoryY, earthTrajectoryZ, 'b', 'LineWidth', 1);  % Plot updated Earth trajectory
+    plot3(marsTrajectoryX, marsTrajectoryY, marsTrajectoryZ, 'r', 'LineWidth', 1);  % Plot updated Mars trajectory
+
+    % Update markers for Earth and Mars
+    set(earthMarker, 'XData', Xearth_AU(1, i), 'YData', Xearth_AU(2, i), 'ZData', Xearth_AU(3, i));  % Update Earth marker position
+    set(marsMarker, 'XData', Xmars_AU(1, i), 'YData', Xmars_AU(2, i), 'ZData', Xmars_AU(3, i));  % Update Mars marker position
+
+    % Calculate and update labels for Earth and Mars
+    v_earth = sqrt(Xearth(4, i)^2 + Xearth(5, i)^2 + Xearth(6, i)^2);  % Calculate Earth's velocity magnitude
+    v_mars = sqrt(Xmars(4, i)^2 + Xmars(5, i)^2 + Xmars(6, i)^2);  % Calculate Mars's velocity magnitude
+    set(earthLabel, 'Position', [Xearth_AU(1, i)+0.1, Xearth_AU(2, i)+0.1, Xearth_AU(3, i)], 'String', sprintf('Earth\n%.2f km/s', v_earth));  % Update Earth label with new velocity
+    set(marsLabel, 'Position', [Xmars_AU(1, i)+0.1, Xmars_AU(2, i)+0.1, Xmars_AU(3, i)], 'String', sprintf('Mars\n%.2f km/s', v_mars));  % Update Mars label with new velocity
+
+    % Convert the ephemeris time to a readable date
+    currentDate = cspice_et2utc(et(i), 'C', 0); % Convert ephemeris time to a UTC string
+
+    % Dynamically update the date legend's text
+    set(dateLegend, 'String', sprintf('Date: %s', currentDate)); % Update the date legend with the current date
+
+    % Adjust speed of animation
+    pause(0.1);  % Pause to slow down the animation for better visualization
+end
+
+%% Clean up by unloading SPICE kernels
+cspice_kclear;  % Unload all SPICE kernels and clear the SPICE subsystem
\ No newline at end of file