LaunchDate='01 Nov 2026 12:30:00 (UTC)';% Define the launch date in UTC
ArrivalDate='07 Sep 2027 12:30:00 (UTC)';% Define the arrival date in UTC
DM=-1;% Transfer mode setting:
% -1 for a Type II trajectory, which takes the longer path, more than 180 degrees around the central body.
% This path may be chosen for specific mission requirements like thermal management or avoiding solar conjunctions.
% +1 for a Type I trajectory, representing the shorter path, less than 180 degrees around the central body.
% This is generally more energy-efficient and quicker, suitable for missions prioritizing speed and fuel efficiency.
%% Data Collection
% Convert both dates to Ephemeris Time using SPICE functions
et0=cspice_str2et(LaunchDate);% Convert the launch date to seconds past J2000
et_end=cspice_str2et(ArrivalDate);% Convert the arrival date to seconds past J2000
% Calculate the total duration in seconds between start and end dates
total_duration=et_end-et0;% Compute the time difference in seconds
% Calculate the number of days between the start and end date
num_days=total_duration/86400;% Convert the duration from seconds to days
% Round num_days to the nearest whole number to avoid fractional days
num_days=round(num_days);% Round the number of days to the nearest integer
% Generate a time vector from start to end date using the calculated number of days
et=linspace(et0,et_end,num_days+1);% Generate evenly spaced time points including both endpoints
% Retrieve state vectors for each time point
fortt=1:length(et)% Loop over each ephemeris time in the vector
Xmars(:,tt)=cspice_spkezr('4',et(tt),'ECLIPJ2000','NONE','10');% Retrieve the state vector of Mars at time tt
Xearth(:,tt)=cspice_spkezr('399',et(tt),'ECLIPJ2000','NONE','10');% Retrieve the state vector of Earth at time tt
% Convert from km to AU
Xmars_AU(:,tt)=Xmars(:,tt)*kmToAU;% Convert Mars' distance from km to AU
Xearth_AU(:,tt)=Xearth(:,tt)*kmToAU;% Convert Earth's distance from km to AU
end
%% Lambert Solver
delta_t=(et_end-et0);% Calculate the time difference between the start and end ephemeris times
epsilon=1e-6;% Set a small epsilon value for numerical stability or small error margin in calculations
% Spacecraft initial position vector
r0=Xearth(1:3,1);% Extract the initial position vector of the spacecraft relative to Earth
% Spacecraft final position vector at Mars
rf=Xmars(1:3,end);% Extract the final position vector of the spacecraft relative to Mars
% Call the Lambert solver function with the gathered data
[v0,vf]=lambert_solver(r0,rf,delta_t,mu,DM);
% Output the initial and final velocities in km/s
disp('Initial velocity vector (km/s):');% Display message indicating that the next line will show the initial velocity vector in km/s.
disp(v0);% Output the initial velocity vector of the spacecraft.
disp('Final velocity vector (km/s):');% Display message indicating that the next line will show the final velocity vector in km/s.
disp(vf);% Output the final velocity vector of the spacecraft.
%% Propagate Orbit
% Define the time span of the simulation
tspan=[0delta_t];% Set the time span for the orbit propagation from 0 to the total duration in seconds.
% Initial state vector (position and velocity)
initial_state=[r0;v0];% Create an initial state vector that includes both position and velocity of the spacecraft.
% Propagate the trajectory using ode45
options=odeset('RelTol',1e-33,'AbsTol',1e-33);% Set options for the ODE solver to use very small tolerance values for high accuracy.
[~,~]=ode45(@(t,y)Orbit_Propagation(t,y,mu),tspan,initial_state,options);% Use the ODE45 solver to simulate the orbit, ignoring the usual output arguments.
% Calculate spacecraft trajectory using ode45
tspan=[0total_duration];% Define the time span from 0 to total mission duration in seconds.
initial_state=[r0;v0];% Define the initial state vector combining the position and velocity vectors.
options=odeset('RelTol',1e-13,'AbsTol',1e-13);% Set the solver options for precision.
[t_out,state_out]=ode45(@(t,y)Orbit_Propagation(t,y,mu),tspan,initial_state,options);% Run the ODE solver with the defined function, time span, initial state, and options.
positions=state_out(:,1:3);% Extract the position vectors from the solution.
positions_AU=positions*kmToAU;% Convert the position vectors from kilometers to Astronomical Units.
% Normalize spacecraft trajectory time to match et
spacecraft_time_normalized=linspace(0,1,length(t_out));% Normalize the time output of the spacecraft trajectory.
et_normalized=linspace(0,1,length(et));% Normalize the ephemeris time vector to the same scale as the spacecraft time.
%% Animation
maxDataPoint=max(abs(positions_AU),[],'all');% Find the maximum absolute value from all data points in the positions matrix for normalization purposes
buffer=0.37*maxDataPoint;% Calculate a buffer (37% of the maximum data point) to ensure all data points are comfortably within the plot boundaries
% Initialize figure and static elements
figure;% Create a new figure for plotting
holdon;% Retain all plots in the figure
gridon;% Enable grid for better visualization of space
axisequal;% Ensure that the scale of the x, y, and z axes are the same
% Set dynamic limits based on the data
xlim([-maxDataPoint-buffer,maxDataPoint+buffer]);% Set the x-axis limits based on the maximum data point plus buffer
ylim([-maxDataPoint-buffer,maxDataPoint+buffer]);% Set the y-axis limits similarly
zlim([-maxDataPoint-buffer,maxDataPoint+buffer]);% Set the z-axis limits similarly
% Determine the title string based on DM
ifDM==1% Check if DM equals 1 (short way)
dmTitle='DM = +1 (Short Way Transfer)';% Assign the title for short way transfer
else% If DM does not equal 1, it must be -1 (long way)
dmTitle='DM = -1 (Long Way Transfer)';% Assign the title for long way transfer
end
xlabel('X Distance (AU)');% Label the x-axis
ylabel('Y Distance (AU)');% Label the y-axis
zlabel('Z Distance (AU)');% Label the z-axis
title(sprintf('Spacecraft Transfer Trajectory from Earth to Mars (%s to %s) (%s)',LaunchDate,ArrivalDate,dmTitle));% Title the plot with the mission timeframe and direction of motion
% Plot Sun
sun=plot3(0,0,0,'yo','MarkerSize',15,'MarkerFaceColor','yellow');% Plot the sun as a yellow circle at the origin
% Initialize plot objects for Earth, Mars, and the Spacecraft trajectory
h_earth=plot3(NaN,NaN,NaN,'b-','LineWidth',1,'MarkerFaceColor','blue');% Initialize the Earth plot with no data
h_mars=plot3(NaN,NaN,NaN,'r-','LineWidth',1,'MarkerFaceColor','red');% Initialize the Mars plot with no data
h_spacecraft=plot3(NaN,NaN,NaN,'m-','LineWidth',1);% Initialize the spacecraft trajectory plot with no data
% Dynamic markers that move along the orbits
earth_marker=plot3(NaN,NaN,NaN,'bo','MarkerSize',5,'MarkerFaceColor','blue');% Marker for Earth
mars_marker=plot3(NaN,NaN,NaN,'ro','MarkerSize',5,'MarkerFaceColor','red');% Marker for Mars
spacecraft_marker=plot3(NaN,NaN,NaN,'m^','MarkerSize',5,'MarkerFaceColor','magenta');% Marker for the spacecraft
% Add legends for the markers after initializing all plot objects
