diff --git a/Animation_Earth_Mars.m b/Animation_Earth_Mars.m new file mode 100644 index 0000000000000000000000000000000000000000..6880e7edfc3c4678514f3af25af4cf280a386bb6 --- /dev/null +++ b/Animation_Earth_Mars.m @@ -0,0 +1,189 @@ +%% 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