%=========================================================================%
%=========================================================================%
% ------------------------------ Eas 4510 ------------------------------- %
% ---------------------------- Spring 2020 -------------------------------%
% -------------------------- Mini Project #4 -----------------------------%
% --------------------------- James Crouse -------------------------------%
clc; clear; close all;
%=========================================================================%
%=========================================================================%
% Question #1 : 30 Points
%=========================================================================%
%developing the two impulse Hohmann transfer
%-------------------------------------------------------------------------%
% earth parameters
mu = 398600; % earth's gravitational parameter
Re = 6371; % earth's radius (km)
M = 500; %number of desired points
%-------------------------------------------------------------------------%
% parameters initial orbit
e0 = 0; % eccentricity is zer b/c circular
alt0 = 500; %altitude of the orbit (km)
r0 = alt0 + Re; %semimajor axis
inc0 = 28.5*pi/180; % orbital inclination (rad)
Omega0 = 60*pi/180; % longitude of the ascending node (rad)
w0 = 0; % argument of the periapsis
nu0 = 0; % true anomaly
%-------------------------------------------------------------------------%
% parameters final orbit
ef = 0; % eccentricity is zer b/c circular
altf = 20200; %altitude of the orbit (km)
rf = altf + Re; %semimajor axis
incf = 57*pi/180; % orbital inclination (rad)
Omegaf = 120*pi/180; % longitude of the ascending node (rad)
wf = 0; % argument of the periapsis
nuf = 0; % true anomaly
[R,t,vv1delta,vv2delta] = twoImpulseHohmann_Crouse_James(r0,rf,inc0,incf,Omega0,Omegaf,mu);
%=========================================================================%
%=========================================================================%
% Question #2 and #3 : 70 Points
%=========================================================================%
N = [1:6]; %number of impulse pairs;
for jj = 1:length(N)
%3D plot
figure (jj);
earthSphere(50)
hold on
[R,t,vdeltp,vdelta,vdelt,tr] = twoNImpulseOrbitTransfer_Crouse_James(r0,rf,inc0,incf,Omega0,Omegaf,N(jj),mu);
%---------------------------------------------------------------------%
%saving variables for plotting
vdp(jj,1:N(jj)) = vdeltp;
vda(jj,1:N(jj)) = vdelta;
t_total(N(jj)) = tr;
%---------------------------------------------------------------------%
%plotting the transfer in 3D
str = sprintf('2xN Impulse Orbit Transfer, N = %g',jj);
t1 = title(str);
set(t1,'Fontweight','bold','FontSize',20,'Interpreter','LaTeX');
xl = xlabel('$x$ (km)');
yl = ylabel('$y$ (km)');
zl = zlabel('$z$ (km)');
set(xl,'FontSize',16,'Interpreter','LaTeX');
set(yl,'FontSize',16,'Interpreter','LaTeX');
set(zl,'FontSize',16,'Interpreter','LaTeX');
set(gca,'FontSize',16,'TickLabelInterpreter','LaTeX');
grid on;
axis equal;
%---------------------------------------------------------------------%
%plotting the mercator display
OmegaE = 2*pi/(24*60*60);
rvValuesECEF = eci2ecef(t,R,OmegaE);
[lonE,lat] = ecef2LonLat(rvValuesECEF);
figure(N(end)+jj);
clf
earth = imread('earth.jpg');
image('CData',earth,'XData',[-180 180],'YData',[90 -90])
hold on
plot(lonE*180/pi,lat*180/pi,'w.','Markersize',8);
str2 = sprintf('2xN Impulse Orbit Transfer, N = %g',jj);
title(str2);
hold off
end
%-------------------------------------------------------------------------%
%plotting impulse at apoapsis as a function of n
figure ((3*jj)+1);
for ii = 1:jj
plot(vda(ii,1:ii),'.-','LineWidth',2,...
'Markersize',12)
hold on
end
hold off
xlabel('n')
ylabel('delta v')
title('Total Impulse at Apoapsis vs n');
legend('N = 1','N = 2','N = 3','N = 4','N = 5','N = 6');
set(gca,'FontSize',16,'TickLabelInterpreter','LaTeX');
grid on;
%-------------------------------------------------------------------------%
%plotting impulse at periapsis as a function of n
figure ((3*jj)+2);
for ii = 1:jj
plot(vdp(ii,1:ii),'.-','LineWidth',2,...
'Markersize',12)
hold on
end
hold off
xlabel('n')
ylabel('delta v')
title('Total Impulse at Periapsis vs n');
legend('N = 1','N = 2','N = 3','N = 4','N = 5','N = 6');
set(gca,'FontSize',16,'TickLabelInterpreter','LaTeX');
grid on;
%-------------------------------------------------------------------------%
%plotting total impulse as a function of N
figure ((3*jj)+3);
for ii = 1:jj
vdt(ii) = sum(vdp(ii,:)) + sum(vdp(ii,:));
end
plot(vdt,'.-','LineWidth',2,...
'Markersize',12)
xlabel('n')
ylabel('delta v')
title('Total Impulse vs N');
set(gca,'FontSize',16,'TickLabelInterpreter','LaTeX');
grid on;
%-------------------------------------------------------------------------%
%plotting time required as a function of N
figure ((3*jj)+4);
plot(N(:),t_total(:),'.-','LineWidth',2,...
'Markersize',12)
xlabel('N')
ylabel('time [s]')
title('Time Required to Complete the Transfer vs N');
set(gca,'FontSize',16,'TickLabelInterpreter','LaTeX');
grid on;
%-------------------------------------------------------------------------%
%print statements
fprintf('The total impulse required to complete the transfer decreases\n');
fprintf('as a function of N. This is because the total energy change of\n');
fprintf('transfer stays constant, and change in energy = sum(1/2*deltav^2\n');
fprintf('+ 2*vinitial*deltav*cos(theta)).During an multiple impulse orbit\n');
fprintf('transfer this second term accounts for more of the total energy\n');
fprintf('change, so the deltav required becomes smaller.\n');
%=========================================================================%
%=========================================================================%