From d81cf4fea85845691496264dfdc9bcf7b38033ae Mon Sep 17 00:00:00 2001
From: Thomas West <tw00956@surrey.ac.uk>
Date: Tue, 3 Sep 2024 09:04:59 +0000
Subject: [PATCH] Trajectories of Earth and Mars around the Sun.

---
 Earth_Mars.m | 296 +++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 296 insertions(+)
 create mode 100644 Earth_Mars.m

diff --git a/Earth_Mars.m b/Earth_Mars.m
new file mode 100644
index 0000000..5f937d8
--- /dev/null
+++ b/Earth_Mars.m
@@ -0,0 +1,296 @@
+%% Script Description (Earth and Mars Positions)
+% This MATLAB script calculates and visualizes the trajectories of Earth and Mars relative to each other
+% and the Sun over a specified period using the SPICE toolkit. The script retrieves precise state vectors
+% (positions and velocities) of Earth and Mars from provided start and end dates, converts distances from
+% kilometers to Astronomical Units (AU), and displays the resulting trajectories in a 3D plot as well as
+% a top-down view. This visualization highlights the dynamic nature of planetary motion within the solar system.
+%
+% Constants
+% Establishes conversion factors for distances (km to AU) and angles (degrees to radians), and defines necessary
+% orbital parameters such as semi-major axes and eccentricities for Earth and Mars, along with the gravitational 
+% parameter of the Sun.
+%
+% Calling SPICE Functions
+% Configures paths and loads necessary SPICE kernels for retrieving planetary data, enabling precise calculations
+% of celestial positions and velocities.
+%
+% Earth Mars State Vectors
+% Calculates the ephemeris time for specified start and end dates and uses these times to retrieve and convert
+% state vectors of Earth and Mars for the entire period. This includes calculating the total duration, number of
+% days between the start and end dates, and generating a time vector to cover the interval.
+%
+% State Vector Analysis
+% Extracts and computes the initial and final position and velocity vectors for both planets. Calculates the
+% magnitude of these vectors to provide insights into the relative speeds and distances at the start and end
+% of the observation period.
+%
+% Free View Plot
+% This 3D plot provides a detailed visualization of the orbits of Earth and Mars, emphasizing changes in their 
+% eccentricity and inclination over the specified period. It marks the positions of the planets at the beginning
+% and end of the observation period, highlighting how these elements change with time.
+%
+% Top Down View
+% Provides a top-down view of the orbits, offering a different perspective on the alignment and spacing
+% between the celestial bodies over time. This view is useful for examining the orbital paths from a planar
+% perspective, highlighting conjunctions and oppositions.
+%
+% 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(groot, 'DefaultAxesFontSize', 24, 'DefaultTextFontSize', 24, '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; % 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
+% Define start and end date for the observation period
+StartDate = '1 Feb 2024 12:00:00 (UTC)'; % Update this line with your desired start date
+EndDate = '20 Mar 2028 12: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 / 86400; % 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 Mars from the SPICE toolkit
+    Xmars(:, tt) = cspice_spkezr('4', et(tt), 'ECLIPJ2000', 'NONE', '10');
+    % Get the state vector of Earth from the SPICE toolkit
+    Xearth(:, tt) = cspice_spkezr('399', et(tt), 'ECLIPJ2000', 'NONE', '10');
+    % Convert from km to AU using conversion factor
+    Xmars_AU(:, tt) = Xmars(:, tt) * kmToAU; % Converts Mars' position from km to AU
+    Xearth_AU(:, tt) = Xearth(:, tt) * kmToAU; % Converts Earth's position from km to AU
+end
+
+%% State Vector
+% Display initial and final state vectors
+% Extract initial position and velocity vectors for Earth and Mars
+initial_earth = Xearth(:, 1);  % Extracts the first column of Xearth, representing the initial state of Earth
+final_earth = Xearth(:, end);  % Extracts the last column of Xearth, representing the final state of Earth
+initial_mars = Xmars(:, 1);    % Extracts the first column of Xmars, representing the initial state of Mars
+final_mars = Xmars(:, end);    % Extracts the last column of Xmars, representing the final state of Mars
+
+% Calculate magnitudes of initial and final position and velocity vectors
+% Calculate the Euclidean norm of the first three elements (position vector)
+initial_earth_pos_mag = sqrt(sum(initial_earth(1:3).^2));
+final_earth_pos_mag = sqrt(sum(final_earth(1:3).^2));
+% Calculate the Euclidean norm of the last three elements (velocity vector)
+initial_earth_vel_mag = sqrt(sum(initial_earth(4:6).^2));
+final_earth_vel_mag = sqrt(sum(final_earth(4:6).^2));
+
+% Repeat for Mars
+initial_mars_pos_mag = sqrt(sum(initial_mars(1:3).^2));
+final_mars_pos_mag = sqrt(sum(final_mars(1:3).^2));
+initial_mars_vel_mag = sqrt(sum(initial_mars(4:6).^2));
+final_mars_vel_mag = sqrt(sum(final_mars(4:6).^2));
+
+% Display initial and final state vectors in a more formatted way using fprintf
+fprintf('Earth Initial State Vector (Position in km and Velocity in km/s):\n');
+fprintf('%12.3f km\n', initial_earth(1:3));  % Display initial position vector of Earth
+fprintf('%12.3f km/s\n', initial_earth(4:6));  % Display initial velocity vector of Earth
+fprintf('Magnitude of Position: %12.3f km\n', initial_earth_pos_mag);  % Display magnitude of initial position
+fprintf('Magnitude of Velocity: %12.3f km/s\n', initial_earth_vel_mag);  % Display magnitude of initial velocity
+fprintf('\n');
+
+fprintf('Earth Final State Vector (Position in km and Velocity in km/s):\n');
+fprintf('%12.3f km\n', final_earth(1:3));  % Display final position vector of Earth
+fprintf('%12.3f km/s\n', final_earth(4:6));  % Display final velocity vector of Earth
+fprintf('Magnitude of Position: %12.3f km\n', final_earth_pos_mag);  % Display magnitude of final position
+fprintf('Magnitude of Velocity: %12.3f km/s\n', final_earth_vel_mag);  % Display magnitude of final velocity
+fprintf('\n');
+
+fprintf('Mars Initial State Vector (Position in km and Velocity in km/s):\n');
+fprintf('%12.3f km\n', initial_mars(1:3));  % Display initial position vector of Mars
+fprintf('%12.3f km/s\n', initial_mars(4:6));  % Display initial velocity vector of Mars
+fprintf('Magnitude of Position: %12.3f km\n', initial_mars_pos_mag);  % Display magnitude of initial position
+fprintf('Magnitude of Velocity: %12.3f km/s\n', initial_mars_vel_mag);  % Display magnitude of initial velocity
+fprintf('\n');
+
+fprintf('Mars Final State Vector (Position in km and Velocity in km/s):\n');
+fprintf('%12.3f km\n', final_mars(1:3));  % Display final position vector of Mars
+fprintf('%12.3f km/s\n', final_mars(4:6));  % Display final velocity vector of Mars
+fprintf('Magnitude of Position: %12.3f km\n', final_mars_pos_mag);  % Display magnitude of final position
+fprintf('Magnitude of Velocity: %12.3f km/s\n', final_mars_vel_mag);  % Display magnitude of final velocity
+fprintf('\n');
+
+%% Free View Plot
+% Create a new figure for plotting
+figure();
+% Plot Earth's orbit in Astronomical Units (AU) as a blue line
+h_earthOrbit = plot3(Xearth_AU(1, :), Xearth_AU(2, :), Xearth_AU(3, :), 'b', 'DisplayName', 'Earth Orbit'); hold on;
+
+% Plot initial position of Earth on its orbit as a blue circle
+h_earthPosition = plot3(Xearth_AU(1, 1), Xearth_AU(2, 1), Xearth_AU(3, 1), 'bo', 'MarkerFaceColor', 'none', 'Color', 'blue', 'DisplayName', sprintf('Earth Initial Position on %s', StartDate));
+
+% Plot final position of Earth on its orbit as a filled blue circle
+h_earthFinal = plot3(Xearth_AU(1, end), Xearth_AU(2, end), Xearth_AU(3, end), 'bo', 'MarkerFaceColor', 'blue', 'DisplayName', sprintf('Earth Final Position on %s', EndDate));
+
+% Plot Mars' orbit in Astronomical Units (AU) as a red line
+h_marsOrbit = plot3(Xmars_AU(1, :), Xmars_AU(2, :), Xmars_AU(3, :), 'r', 'DisplayName', 'Mars Orbit');
+
+% Plot initial position of Mars on its orbit as a red circle
+h_marsPosition = plot3(Xmars_AU(1, 1), Xmars_AU(2, 1), Xmars_AU(3, 1), 'ro', 'MarkerFaceColor', 'none', 'Color', 'red', 'DisplayName', sprintf('Mars Initial Position on %s', StartDate));
+
+% Plot final position of Mars on its orbit as a filled red circle
+h_marsFinal = plot3(Xmars_AU(1, end), Xmars_AU(2, end), Xmars_AU(3, end), 'ro', 'MarkerFaceColor', 'red', 'DisplayName', sprintf('Mars Final Position on %s', EndDate));
+
+% Plot the Sun at the origin of the coordinate system as a yellow filled circle
+h_sun = plot3(0, 0, 0, 'yo', 'MarkerSize', 10, 'MarkerFaceColor', 'yellow', 'DisplayName', 'Sun');
+
+% Settings for velocity vector plots (arrows)
+scaleFactor = 0.0025; % Scale factor to adjust arrow sizes for visibility
+maxHeadSize = 10; % Maximum size of arrow heads
+
+% Plot velocity vectors of Earth and Mars at initial positions
+quiver3(Xearth_AU(1, 1), Xearth_AU(2, 1), Xearth_AU(3, 1), Xearth(4, 1) * scaleFactor, Xearth(5, 1) * scaleFactor, Xearth(6, 1) * scaleFactor, 'b', 'LineWidth', 1.5, 'MaxHeadSize', maxHeadSize, 'AutoScale', 'off');
+quiver3(Xmars_AU(1, 1), Xmars_AU(2, 1), Xmars_AU(3, 1), Xmars(4, 1) * scaleFactor, Xmars(5, 1) * scaleFactor, Xmars(6, 1) * scaleFactor, 'r', 'LineWidth', 1.5, 'MaxHeadSize', maxHeadSize, 'AutoScale', 'off');
+
+% Enhancements for 3D Visualization
+xlabel('X Distance (AU)');  % Label for the X-axis
+ylabel('Y Distance (AU)');  % Label for the Y-axis
+zlabel('Z Distance (AU)');  % Label for the Z-axis
+title(sprintf('3D Orbit Trajectories of Earth and Mars from %s to %s', StartDate, EndDate));  % Title of the plot
+legend([h_earthOrbit, h_earthPosition, h_earthFinal, h_marsOrbit, h_marsPosition, h_marsFinal, h_sun], 'Location', 'bestoutside');  % Legend to identify plot elements
+grid on;  % Enable grid for better visualization
+axis equal;  % Equal scaling for all axes
+xlim([-2, 2]);  % X-axis limits
+ylim([-2, 2]);  % Y-axis limits
+zlim([-2, 2]);  % Z-axis limits
+view(3);  % Set the view to 3D perspective
+
+%% Top Down View
+% Create a new figure for the top-down plot
+figure();
+% Plot Earth's orbit in Astronomical Units (AU) as a blue line
+h_earthOrbit = plot3(Xearth_AU(1, :), Xearth_AU(2, :), Xearth_AU(3, :), 'b', 'DisplayName', 'Earth Orbit'); hold on;
+
+% Plot initial position of Earth on its orbit as a blue circle
+h_earthPosition = plot3(Xearth_AU(1, 1), Xearth_AU(2, 1), Xearth_AU(3, 1), 'bo', 'MarkerFaceColor', 'none', 'Color', 'blue', 'DisplayName', sprintf('Earth Initial Position on %s', StartDate));
+
+% Plot final position of Earth on its orbit as a filled blue circle
+h_earthFinal = plot3(Xearth_AU(1, end), Xearth_AU(2, end), Xearth_AU(3, end), 'bo', 'MarkerFaceColor', 'blue', 'DisplayName', sprintf('Earth Final Position on %s', EndDate));
+
+% Plot Mars' orbit in Astronomical Units (AU) as a red line
+h_marsOrbit = plot3(Xmars_AU(1, :), Xmars_AU(2, :), Xmars_AU(3, :), 'r', 'DisplayName', 'Mars Orbit');
+
+% Plot initial position of Mars on its orbit as a red circle
+h_marsPosition = plot3(Xmars_AU(1, 1), Xmars_AU(2, 1), Xmars_AU(3, 1), 'ro', 'MarkerFaceColor', 'none', 'Color', 'red', 'DisplayName', sprintf('Mars Initial Position on %s', StartDate));
+
+% Plot final position of Mars on its orbit as a filled red circle
+h_marsFinal = plot3(Xmars_AU(1, end), Xmars_AU(2, end), Xmars_AU(3, end), 'ro', 'MarkerFaceColor', 'red', 'DisplayName', sprintf('Mars Final Position on %s', EndDate));
+
+% Plot the Sun at the origin of the coordinate system as a yellow filled circle
+h_sun = plot3(0, 0, 0, 'yo', 'MarkerSize', 10, 'MarkerFaceColor', 'yellow', 'DisplayName', 'Sun');
+
+% Quiver settings for plotting velocity vectors as arrows showing orbit direction
+scaleFactor = 0.0025; % Scale factor to adjust arrow sizes for better visibility
+maxHeadSize = 10; % Maximum size of arrow heads for visibility
+
+% Plot velocity vectors of Earth and Mars at initial positions
+quiver3(Xearth_AU(1, 1), Xearth_AU(2, 1), Xearth_AU(3, 1), Xearth(4, 1) * scaleFactor, Xearth(5, 1) * scaleFactor, Xearth(6, 1) * scaleFactor, 'b', 'LineWidth', 1.5, 'MaxHeadSize', maxHeadSize, 'AutoScale', 'off');
+quiver3(Xmars_AU(1, 1), Xmars_AU(2, 1), Xmars_AU(3, 1), Xmars(4, 1) * scaleFactor, Xmars(5, 1) * scaleFactor, Xmars(6, 1) * scaleFactor, 'r', 'LineWidth', 1.5, 'MaxHeadSize', maxHeadSize, 'AutoScale', 'off');
+
+% Enhancements for 3D Visualization in a 2D projection
+xlabel('X Distance (AU)');  % Label for the X-axis
+ylabel('Y Distance (AU)');  % Label for the Y-axis
+zlabel('Z Distance (AU)');  % Label for the Z-axis (though not visible in top-down view)
+title(sprintf('3D Orbit Trajectories of Earth and Mars from %s to %s', StartDate, EndDate)); % Title of the plot
+legend([h_earthOrbit, h_earthPosition, h_earthFinal, h_marsOrbit, h_marsPosition, h_marsFinal, h_sun], 'Location', 'bestoutside'); % Legend to identify plot elements
+grid on;  % Enable grid for better visualization
+axis equal;  % Equal scaling for all axes
+xlim([-2, 2]);  % X-axis limits
+ylim([-2, 2]);  % Y-axis limits
+zlim([-2, 2]);  % Z-axis limits (mostly for maintaining aspect ratio)
+view(0, 90);  % Set the view to top-down (along the z-axis)
+
+%% Top Down View
+% Create a new figure for the top-down plot
+figure();
+% Plot Earth's orbit in Astronomical Units (AU) as a blue line
+h_earthOrbit = plot3(Xearth_AU(1, :), Xearth_AU(2, :), Xearth_AU(3, :), 'b', 'DisplayName', 'Earth Orbit'); hold on;
+
+% Plot initial position of Earth on its orbit as a blue circle
+h_earthPosition = plot3(Xearth_AU(1, 1), Xearth_AU(2, 1), Xearth_AU(3, 1), 'bo', 'MarkerFaceColor', 'none', 'Color', 'blue', 'DisplayName', sprintf('Earth Initial Position on %s', StartDate));
+
+% Plot final position of Earth on its orbit as a filled blue circle
+h_earthFinal = plot3(Xearth_AU(1, end), Xearth_AU(2, end), Xearth_AU(3, end), 'bo', 'MarkerFaceColor', 'blue', 'DisplayName', sprintf('Earth Final Position on %s', EndDate));
+
+% Plot Mars' orbit in Astronomical Units (AU) as a red line
+h_marsOrbit = plot3(Xmars_AU(1, :), Xmars_AU(2, :), Xmars_AU(3, :), 'r', 'DisplayName', 'Mars Orbit');
+
+% Plot initial position of Mars on its orbit as a red circle
+h_marsPosition = plot3(Xmars_AU(1, 1), Xmars_AU(2, 1), Xmars_AU(3, 1), 'ro', 'MarkerFaceColor', 'none', 'Color', 'red', 'DisplayName', sprintf('Mars Initial Position on %s', StartDate));
+
+% Plot final position of Mars on its orbit as a filled red circle
+h_marsFinal = plot3(Xmars_AU(1, end), Xmars_AU(2, end), Xmars_AU(3, end), 'ro', 'MarkerFaceColor', 'red', 'DisplayName', sprintf('Mars Final Position on %s', EndDate));
+
+% Plot the Sun at the origin of the coordinate system as a yellow filled circle
+h_sun = plot3(0, 0, 0, 'yo', 'MarkerSize', 10, 'MarkerFaceColor', 'yellow', 'DisplayName', 'Sun');
+
+% Quiver settings for plotting velocity vectors as arrows showing orbit direction
+scaleFactor = 0.0025; % Scale factor to adjust arrow sizes for better visibility
+maxHeadSize = 10; % Maximum size of arrow heads for visibility
+
+% Plot velocity vectors of Earth and Mars at initial positions
+quiver3(Xearth_AU(1, 1), Xearth_AU(2, 1), Xearth_AU(3, 1), Xearth(4, 1) * scaleFactor, Xearth(5, 1) * scaleFactor, Xearth(6, 1) * scaleFactor, 'b', 'LineWidth', 1.5, 'MaxHeadSize', maxHeadSize, 'AutoScale', 'off');
+quiver3(Xmars_AU(1, 1), Xmars_AU(2, 1), Xmars_AU(3, 1), Xmars(4, 1) * scaleFactor, Xmars(5, 1) * scaleFactor, Xmars(6, 1) * scaleFactor, 'r', 'LineWidth', 1.5, 'MaxHeadSize', maxHeadSize, 'AutoScale', 'off');
+
+% Enhancements for 3D Visualization in a 2D projection
+xlabel('X Distance (AU)');  % Label for the X-axis
+ylabel('Y Distance (AU)');  % Label for the Y-axis
+zlabel('Z Distance (AU)');  % Label for the Z-axis (though not visible in top-down view)
+title(sprintf('3D Orbit Trajectories of Earth and Mars from %s to %s', StartDate, EndDate)); % Title of the plot
+legend([h_earthOrbit, h_earthPosition, h_earthFinal, h_marsOrbit, h_marsPosition, h_marsFinal, h_sun], 'Location', 'bestoutside'); % Legend to identify plot elements
+grid on;  % Enable grid for better visualization
+axis equal;  % Equal scaling for all axes
+xlim([-2, 2]);  % X-axis limits
+ylim([-2, 2]);  % Y-axis limits
+zlim([-2, 2]);  % Z-axis limits (mostly for maintaining aspect ratio)
+view(0, 0);  % Set the view to top-down (along the z-axis)
+
+%% Clean up by unloading SPICE kernels
+cspice_kclear;  % Unload all SPICE kernels and clear the SPICE subsystem
\ No newline at end of file
-- 
GitLab