Mini-Project 3

%=========================================================================%

%---------------------------| EAS4510 |---------------------------%

%---------------------------| Mini Project #3 |---------------------------%

%---------------------------| Spring 2020 |---------------------------%

%---------------------------| James Crouse |---------------------------%

%=========================================================================%

%=========================================================================%

close all; clc; clear;

%=========================================================================%

%=========================================================================%

%Given Parameters

mu = 398600; %gravitational parameter

tau = 24*60*60; %period in seconds

e = 0.25; %eccentricity

inc = 63.4*pi/180; %orbital inclination in radians

w = [90 270].*pi/180; %argument of periapsis in radians

Omega = [0:45:180].*pi/180; %Longitude of the ascending node in radians

%setting up the time array

t0 = 0; %epoch time

deltat = 5*60; %timesteps are 5 minutes

N = tau/deltat; %number of time steps

t = linspace(t0,tau,N); %creates the time array

%calculationg rPCI0 and vPCI0

a = ((tau/(2*pi))^2*mu)^(1/3);

nu0 = 0;%because the initial time is = epoch time

for kk = 1:length(w) %this loop cycles through the different values of w

for jj =1:length(Omega) %this loop cycles through the different values of Omega

oe = [a e Omega(jj) inc w(kk) nu0]';

[rPCI0,vPCI0] = oe2rv_Crouse_James(oe,mu);

rprop(1,:) = rPCI0';

vprop(1,:) = vPCI0';

for ii = 2:N

t0 = t(ii-1);

tf = t(ii);

rPCI0 = rprop(ii-1,:);

vPCI0 = vprop(ii-1,:);

[r,v] = KeplerSolver_Crouse_James(rPCI0,vPCI0,mu,t0,tf);

rprop(ii,:) = r';

vprop(ii,:) = v';

end

%task 2-----------------------------------------------------------%

xprop = rprop(:,1); %position arrays

yprop = rprop(:,2);

zprop = rprop(:,3);

vxprop = vprop(:,1); %velocity arrays

vyprop = vprop(:,2);

vzprop = vprop(:,3);

column = (5*(kk-1))+jj;%sets up a column to be able to store variables to an array

rvValuesECI = [xprop yprop zprop];

tValues = t';

OmegaE = 2*pi/(24*60*60); %angular velocity of the earth

rvValuesECEF = eci2ecef(tValues,rvValuesECI,OmegaE);

[lonE(:,column),lat(:,column)] = ecef2LonLat(rvValuesECEF);

xECI(:,column) = xprop;

yECI(:,column) = yprop;

zECI(:,column) = zprop;

end

end

%finds the index corresponding to the apoapsis for each case

for ii = 1:(length(w)*length(Omega))

radius(:,ii) = sqrt(xECI(:,ii).^2+xECI(:,ii).^2+xECI(:,ii).^2);

[maxrad(ii),index(ii)] = max(radius(:,ii));

end

%plotting

for ii = 1:length(w)

if ii == 1

column = 1:5;

num = 1;

else

column = 6:10;

num = 3;

end

Markers = {'o','s','d','v','>','o','s','d','v','>'};%want to label each plot in

%a figure with a different marker for clarity.

%plotting the position

figure (num);

for kk = column(1):column(5)

plot3(xECI(:,kk),yECI(:,kk),zECI(:,kk),Markers{kk})

hold on

end

axis equal

str = sprintf('Position for w = %g [deg]',w(ii)*180/pi);

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;

hold on;

earthSphere(50)%plots earthSphere function

ll = legend('Omega = 0 [deg]','Omega = 45 [deg]','Omega = 90 [deg]','Omega = 135 [deg]','Omega = 180 [deg]','Location','northeast');

set(ll,'FontSize',16,'Interpreter','LaTeX');

hold off

%plotting the groundtrack

figure(num+1);

clf

earth = imread('earth.jpg');

image('CData',earth,'XData',[-180 180],'YData',[90 -90])

hold on

for kk = column(1):column(5)

plot(lonE(:,kk)*180/pi,lat(:,kk)*180/pi,Markers{kk});%plots the groundtrack for each case

hold on

end

% plot(lonE(index(column),column)*180/pi,lat(index(column),column)*180/pi,'wp');%plots the apoapsis for each case

hold off

str = sprintf('Groundtrack for w = %g [deg]',w(ii)*180/pi);

t1 = title(str);

set(t1,'Fontweight','bold','FontSize',20,'Interpreter','LaTeX');

ll = legend('Omega = 0 [deg]','Omega = 45 [deg]','Omega = 90 [deg]','Omega = 135 [deg]','Omega = 180 [deg]','Location','northeast');

set(ll,'FontSize',16,'Interpreter','LaTeX');

end

fprintf('<strong>Task 2</strong>\n');

fprintf('====================================================================================\n')

fprintf('the two cases, w = 90 and 270 [deg], are mirrored about lattitude = 0\n');

fprintf('with respect to the longitudinal distance of the apoapsis to the equator\n');

fprintf('====================================================================================\n')

%finding the time spent in each hemisphere for one of the orbits this can

%be done because the cases where w = 90 and 270 are mirror image;

%therefore, the time spent in the northern hemisphere for case 1 = time

%spent in the southern hemisphere for case 2.

lattitude = lat(:,1);

for ii = 2:length(lattitude)

if lattitude(ii-1) >= 0 && lattitude(ii)<= 0

index1 = ii;

t1 = t(ii);

elseif lattitude(ii) >= 0 && lattitude(ii-1)<= 0

index2 = ii;

t2 = t(ii);

else

end

end

tnorth = abs(t1-t2)/3600;

tsouth = (tau/3600) - tnorth;

fprintf('<strong>Task 3</strong>\n');

fprintf('====================================================================================\n')

fprintf('<strong>a:</strong>\n');

fprintf('w = 90\n')

fprintf('most of the orbit occurs in the Northern hemisphere\n');

fprintf('the spacecraft spends %g [hours] from 90 to 270 degrees\n',tnorth);

fprintf('------------------------------------------------------------------------------------\n')

fprintf('w = 270\n')

fprintf('most of the orbit occurs in the Southern hemisphere\n');

fprintf('the spacecraft spends %g [hours] from 90 to 270 degrees\n',tsouth);

fprintf('====================================================================================\n')

fprintf('<strong>b:</strong>\n');

fprintf('The orbit when w = 90 [deg] and Omega = 90 [deg] has an apoapsis over Australia\n')

fprintf('The orbit when w = 90 [deg] and Omega = 180 [deg] has an apoapsis over South America\n')

fprintf('====================================================================================\n')

fprintf('<strong>c:</strong>\n');

fprintf('The orbit when w = 270 [deg] and Omega = 0 [deg] has an apoapsis over North America\n')

fprintf('The orbit when w = 270 [deg] and Omega = 180 [deg] has an apoapsis over Russia\n')

fprintf('====================================================================================\n')

%=========================================================================%

%=========================================================================%