LaunchDate='1 Nov 2026 12:30:00 (UTC)';% Define the launch date in UTC
ArrivalDate='7 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;% 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 time point
Xmars(:,tt)=cspice_spkezr('4',et(tt),'ECLIPJ2000','NONE','10');% Retrieve Mars' state vector at time tt
Xearth(:,tt)=cspice_spkezr('399',et(tt),'ECLIPJ2000','NONE','10');% Retrieve Earth's state vector at time tt
% Convert from km to AU
Xmars_AU(:,tt)=Xmars(:,tt)*kmToAU;% Convert Mars' position from km to Astronomical Units
Xearth_AU(:,tt)=Xearth(:,tt)*kmToAU;% Convert Earth's position from km to Astronomical Units
end
%% Lambert Solver
delta_t=(et_end-et0);% Calculate the transfer time in seconds
epsilon=1e-6;% Set a small number for numerical precision in calculations
% Spacecraft initial position vector
r0=Xearth(1:3,1);% Extract Earth's position vector at launch
% Spacecraft final position vector at Mars
rf=Xmars(1:3,end);% Extract Mars' position vector at the arrival time
% Call the Lambert solver function with the gathered data
[v0,vf]=lambert_solver(r0,rf,delta_t,mu,DM);% Compute initial and final velocity vectors using the Lambert solver
% Output the initial and final velocities in km/s
disp('Initial velocity vector (km/s):');% Display message for initial velocity
disp(v0);% Display the computed initial velocity vector
disp('Final velocity vector (km/s):');% Display message for final velocity
disp(vf);% Display the computed final velocity vector
% Define the time span of the simulation with more points
tspan=linspace(0,delta_t,num_days);% Generate a finely spaced time vector for the simulation
% Initial state vector (position and velocity)
initial_state=[r0;v0];% Combine the position vector and initial velocity vector into a single state vector
%% Propagate Orbit
% Propagate the trajectory using ode45
options=odeset('RelTol',1e-13,'AbsTol',1e-13);% Set options for numerical solver to high precision
[t_out,state_out]=ode45(@(t,y)Orbit_Propagation(t,y,mu),tspan,initial_state,options);% Numerically integrate the trajectory
% Extract the position vectors from the output state for plotting
positions=state_out(:,1:3);% Extract position data from the output state for subsequent use in plotting
velocities=state_out(:,4:6);
%% Calculate Distances and Find Exit from Earth's SOI
SOI_distance_km=924000;% Define the distance for SOI exit
exit_index=0;% Initialize exit index
% Loop through all time points to calculate distances
fortt=1:length(t_out)
% Calculate the distance between Earth and the spacecraft at time tt
distance=norm(positions(tt,:)-Xearth(1:3,tt)');% Make sure Xearth is interpolated if necessary
% Check if the distance is greater than the SOI distance
ifdistance>SOI_distance_km
exit_index=tt;% Save the index when spacecraft exits the SOI
break;% Exit the loop as we found the first exit point
end
end
ifexit_index>0
% Earth's state vector at exit time (assuming Earth data is interpolated to match spacecraft data points)
Earth_position=Xearth(1:3,exit_index);% Earth's position vector in km
Earth_velocity=Xearth(4:6,exit_index);% Earth's velocity vector in km/s
% Spacecraft's heliocentric state vector at exit time
Spacecraft_position=positions(exit_index,:);% Spacecraft's position vector in km
Spacecraft_velocity=velocities(exit_index,:);% Spacecraft's velocity vector in km/s
% Convert to geocentric frame
Geocentric_position=Spacecraft_position-Earth_position';% Convert position from heliocentric to geocentric
Geocentric_velocity=Spacecraft_velocity-Earth_velocity';% Convert velocity from heliocentric to geocentric
% Time when the spacecraft exits the SOI in seconds and hours
exit_time_seconds=t_out(exit_index);% Time in seconds since the start of the simulation
exit_time_hours=exit_time_seconds/3600;% Convert time to hours
exit_time_days=exit_time_seconds/86400;% Convert time to days
% Display the results
fprintf('Time of exit from Earths SOI: %f seconds (%f hours, %f days)\n',exit_time_seconds,exit_time_hours,exit_time_days);
fprintf('Earths position vector at exit time (km): [%f, %f, %f]\n',Earth_position);
fprintf('Earths velocity vector at exit time (km): [%f, %f, %f]\n',Earth_velocity);
fprintf('Spacecrafts heliocentric position at exit time (km): [%f, %f, %f]\n',Spacecraft_position);
fprintf('Spacecrafts heliocentric velocity at exit time (km/s): [%f, %f, %f]\n',Spacecraft_velocity);
fprintf('Geocentric position vector of the spacecraft at exit time (km): [%f, %f, %f]\n',Geocentric_position);
fprintf('Geocentric velocity vector of the spacecraft at exit time (km/s): [%f, %f, %f]\n',Geocentric_velocity);
else
fprintf('Spacecraft does not exit Earths SOI within the simulation period.\n');
end
%$ Mars proximity threshold in kilometers
Mars_proximity_km=577300;% Define the proximity distance for Mars approach
approach_index=0;% Initialize approach index
% Loop through all time points to calculate distances
fortt=1:length(t_out)
% Calculate the distance between Mars and the spacecraft at time tt
distance=norm(positions(tt,:)-Xmars(1:3,tt)');% Ensure Xmars is interpolated
% Check if the distance is less than the Mars proximity distance
ifdistance<Mars_proximity_km
approach_index=tt;% Save the index when spacecraft approaches Mars
break;% Exit the loop as we found the first approach point
end
end
% Output results
ifapproach_index>0
approach_time_seconds=t_out(approach_index);% Time in seconds since the start of the simulation
approach_time_hours=approach_time_seconds/3600;% Convert time to hours
approach_time_days=approach_time_seconds/86400;% Convert time to days
% Extract Mars and spacecraft state vectors at approach time
Mars_position=Xmars(1:3,approach_index);% Mars' position vector in km
Mars_velocity=Xmars(4:6,approach_index);% Mars' velocity vector in km/s
Spacecraft_position=positions(approach_index,:);% Spacecraft's position vector in km
Spacecraft_velocity=velocities(approach_index,:);% Spacecraft's velocity vector in km/s
% Convert to areocentric frame
Areocentric_position=Spacecraft_position-Mars_position';% Convert position from heliocentric to areocentric
Areocentric_velocity=Spacecraft_velocity-Mars_velocity';% Convert velocity from heliocentric to areocentric
fprintf('Spacecraft approaches Mars within %f km at t = %f seconds (%f hours, %f days)\n',Mars_proximity_km,approach_time_seconds,approach_time_hours,approach_time_days);
fprintf('Mars position vector at entry time (km): [%f, %f, %f]\n',Mars_position);
fprintf('Mars velocity vector at entry time (km): [%f, %f, %f]\n',Mars_velocity);
fprintf('Spacecrafts heliocentric position at entry time (km): [%f, %f, %f]\n',Spacecraft_position);
fprintf('Spacecrafts heliocentric velocity at entry time (km/s): [%f, %f, %f]\n',Spacecraft_velocity);
fprintf('Areocentric position vector of the spacecraft at approach time (km): [%f, %f, %f]\n',Areocentric_position);
fprintf('Areocentric velocity vector of the spacecraft at approach time (km/s): [%f, %f, %f]\n',Areocentric_velocity);
else
fprintf('Spacecraft does not approach Mars within %f km during the simulation period.\n',Mars_proximity_km);
end
% Assuming Xearth and v0 are defined
% Extract the vector from Xearth
vec1=Xearth(4:6,1);
% Assume v0 is already defined as a vector
vec2=v0;
% Calculate the dot product of vec1 and vec2
dot_product=dot(vec1,vec2);
% Calculate the magnitudes of vec1 and vec2
mag_vec1=norm(vec1);
mag_vec2=norm(vec2);
% Calculate the cosine of the angle
cos_theta=dot_product/(mag_vec1*mag_vec2);
% Calculate the angle in degrees
angle=acosd(cos_theta);
% Display the angle
disp(['The angle between the vectors is ',num2str(angle),' degrees.']);
v0
Xearth(4:6,1)
%% Free View Plot
% Use start and end dates from earlier in the script for labeling
StartDate=LaunchDate;
EndDate=ArrivalDate;
% Get the initial and final positions for Earth and Mars in Astronomical Units
earthInitialPos=Xearth_AU(1:3,1);% Extract the initial position of Earth
earthFinalPos=Xearth_AU(1:3,end);% Extract the final position of Earth
marsInitialPos=Xmars_AU(1:3,1);% Extract the initial position of Mars
marsFinalPos=Xmars_AU(1:3,end);% Extract the final position of Mars
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
% Enhancements for 3D Visualization
xlabel('X Distance (AU)');% Label the x-axis of the plot
ylabel('Y Distance (AU)');% Label the y-axis of the plot
zlabel('Z Distance (AU)');% Label the z-axis of the plot
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
legend([h_earthOrbit,h_marsOrbit,h_earthPosition,h_marsPosition,h_spacecraftTrajectory,h_sun],'Location','bestoutside');% Create a legend and place it outside the plot
gridon;% Enable grid for better visibility of plot scale and distances
axisequal;% Set all axis scales to be equal
xlim([-2,2]);% Set limits for the x-axis
ylim([-2,2]);% Set limits for the y-axis
zlim([-2,2]);% Set limits for the z-axis
view(3);% Set the view angle to 3D for a better perspective of the plot
% Ensure the plot holds for additional plotting
holdoff;% Release the hold on the current plot
%% Clean up by unloading SPICE kernels
cspice_kclear;% Unload all SPICE kernels and clear the SPICE subsystem
