Skip to content
Snippets Groups Projects
Commit d81cf4fe authored by Thomas West's avatar Thomas West
Browse files

Trajectories of Earth and Mars around the Sun.

parent b6a061da
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 296 additions and 0 deletions
Earth_Mars.m 0 → 100644
+ 296
− 0
View file @ d81cf4fe
%% 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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment