diff --git a/Earth_Mars.m b/Earth_Mars.m new file mode 100644 index 0000000000000000000000000000000000000000..5f937d833e3798705976b3664d17675a562944c1 --- /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