List of Functions

function TI2E = dcmI2E(t,OmegaE)

% --------------------------------------------------------------%

% -------------------------- DCMI2E.M --------------------------%

% --------------------------------------------------------------%

% Compute the transformation matrix from Earth-centered inertial%

% Cartesian coordinates to Earth-centered Earth-fixed Cartesian %

% coordinates. The transformation matrix is of size 3 by 3 and %

% represents the transformation matrix at the time t %

% Inputs: %

% t: the time at which the transformation is desired %

% OmegaE: the rotation rate of the Earth %

% Outputs: %

% TI2E: 3 by 3 matrix that transforms components of a %

% vector expressed in Earth-centered inertial %

% Cartesian coordinates to Earth-centered %

% Earth-fixed Cartesian coordinates %

% --------------------------------------------------------------%

% NOTE: the units of the inputs T and OMEGAE must be consistent %

% and the angular rate must be in RADIANS PER TIME UNIT %

% --------------------------------------------------------------%

TI2E = [cos(OmegaE*t) -sin(OmegaE*t) 0; sin(OmegaE*t) cos(OmegaE*t) 0; 0 0 1]';

end

function [xx,yy,zz] = earthSphere(varargin)

%EARTH_SPHERE Generate an earth-sized sphere.

% [X,Y,Z] = EARTH_SPHERE(N) generates three (N+1)-by-(N+1)

% matrices so that SURFACE(X,Y,Z) produces a sphere equal to

% the radius of the earth in kilometers. The continents will be

% displayed.

%

% [X,Y,Z] = EARTH_SPHERE uses N = 50.

%

% EARTH_SPHERE(N) and just EARTH_SPHERE graph the earth as a

% SURFACE and do not return anything.

%

% EARTH_SPHERE(N,'mile') graphs the earth with miles as the unit rather

% than kilometers. Other valid inputs are 'ft' 'm' 'nm' 'miles' and 'AU'

% for feet, meters, nautical miles, miles, and astronomical units

% respectively.

%

% EARTH_SPHERE(AX,...) plots into AX instead of GCA.

%

% Examples:

% earth_sphere('nm') produces an earth-sized sphere in nautical miles

%

% earth_sphere(10,'AU') produces 10 point mesh of the Earth in

% astronomical units

%

% h1 = gca;

% earth_sphere(h1,'mile')

% hold on

% plot3(x,y,z)

% produces the Earth in miles on axis h1 and plots a trajectory from

% variables x, y, and z

% Clay M. Thompson 4-24-1991, CBM 8-21-92.

% Will Campbell, 3-30-2010

% Copyright 1984-2010 The MathWorks, Inc.

%% Input Handling

[cax,args,nargs] = axescheck(varargin{:}); % Parse possible Axes input

error(nargchk(0,2,nargs)); % Ensure there are a valid number of inputs

% Handle remaining inputs.

% Should have 0 or 1 string input, 0 or 1 numeric input

j = 0;

k = 0;

n = 50; % default value

units = 'km'; % default value

for i = 1:nargs

if ischar(args{i})

units = args{i};

j = j+1;

elseif isnumeric(args{i})

n = args{i};

k = k+1;

end

end

if j > 1 || k > 1

error('Invalid input types')

end

%% Calculations

% Scale factors

Scale = {'km' 'm' 'mile' 'miles' 'nm' 'au' 'ft';

1 1000 0.621371192237334 0.621371192237334 0.539956803455724 6.6845871226706e-009 3280.839895};

% Identify which scale to use

try

myscale = 6378.1363*Scale{2,strcmpi(Scale(1,:),units)};

catch %#ok<*CTCH>

error('Invalid units requested. Please use m, km, ft, mile, miles, nm, or AU')

end

% -pi <= theta <= pi is a row vector.

% -pi/2 <= phi <= pi/2 is a column vector.

theta = (-n:2:n)/n*pi;

phi = (-n:2:n)'/n*pi/2;

cosphi = cos(phi); cosphi(1) = 0; cosphi(n+1) = 0;

sintheta = sin(theta); sintheta(1) = 0; sintheta(n+1) = 0;

x = myscale*cosphi*cos(theta);

y = myscale*cosphi*sintheta;

z = myscale*sin(phi)*ones(1,n+1);

%% Plotting

if nargout == 0

cax = newplot(cax);

% Load and define topographic data

load('topo.mat','topo','topomap1');

% Rotate data to be consistent with the Earth-Centered-Earth-Fixed

% coordinate conventions. X axis goes through the prime meridian.

% http://en.wikipedia.org/wiki/Geodetic_system#Earth_Centred_Earth_Fixed_.28ECEF_or_ECF.29_coordinates

%

% Note that if you plot orbit trajectories in the Earth-Centered-

% Inertial, the orientation of the contintents will be misleading.

topo2 = [topo(:,181:360) topo(:,1:180)]; %#ok<NODEF>

% Define surface settings

props.FaceColor= 'texture';

props.EdgeColor = 'none';

props.FaceLighting = 'phong';

props.Cdata = topo2;

% Create the sphere with Earth topography and adjust colormap

surface(x,y,z,props,'parent',cax)

colormap(topomap1)

% Replace the calls to surface and colormap with these lines if you do

% not want the Earth's topography displayed.

% surf(x,y,z,'parent',cax)

% shading flat

% colormap gray

% Refine figure

axis equal

xlabel(['X [' units ']'])

ylabel(['Y [' units ']'])

zlabel(['Z [' units ']'])

view(127.5,30)

else

xx = x; yy = y; zz = z;

end

function [lonE,lat] = ecef2LonLat(rvValuesECEF)

% --------------------------------------------------------------%

% ----------------------- ECEF2LONLAT.M ------------------------%

% --------------------------------------------------------------%

% Compute the Earth-relative longitude and geocentric latitude %

% from the Earth-centered Earth-fixed position. %

% Input: %

% rvValuesECEF: matrix of values of the position expressed %

% in Earth-centered Earth-fixed Cartesian %

% and stored ROW-WISE %

% Outputs: %

% lonE: a column matrix containing the values of the %

% Earth-relative longitude (radians) %

% lat: a column matrix containing the values of the %

% geocentric latitude (radians) %

% --------------------------------------------------------------%

for ii = 1:length(rvValuesECEF)

x = rvValuesECEF(ii,1);

y = rvValuesECEF(ii,2);

z = rvValuesECEF(ii,3);

lonE(ii) = atan2(y,x);

lat(ii) = atan2(z,sqrt(x^2+y^2));

end

end

function rvValuesECEF = eci2ecef(tValues,rvValuesECI,OmegaE)

% --------------------------------------------------------------%

% ------------------------- ECI2ECEF.M -------------------------%

% --------------------------------------------------------------%

% Transform the values of the position along a spacecraft orbit %

% from Earth-centered inertial Cartesian coordinates to %

% Earth-centered Earth-fixed Cartesian coordinates. %

% represents the transformation matrix at the time t %

% Inputs: %

% tValues: column matrix of values of time at which the %

% transformation is to be performed %

% rvValuesECI: matrix of values of the position expressed %

% in Earth-centered inertial Cartesian %

% coordinates and stored ROW-WISE %

% OmegaE: the rotation rate of the Earth %

% Outputs: %

% rvValuesECEF: matrix of values of the position expressed %

% in Earth-centered Earth-fixed Cartesian %

% and stored ROW-WISE %

% --------------------------------------------------------------%

for ii = 1:length(tValues)

TI2E = dcmI2E(tValues(ii), OmegaE);

rvValuesECEF(ii,:) = TI2E * rvValuesECI(ii,:)';

end

end

function [r,v,E,nu] = KeplerSolver_Crouse_James(rPCI0,vPCI0,mu,t0,t)

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

% This function is a general purpose code that determines the position and%

% inertial velocity of a spacecraft in PCI coordinatesat an arbitrary time%

% t, given the position and velocity at a time t0. (0 <= e < 1) %

% %

%Inputs: rPCI0 = column vector containing the position at t = t0 %

% vPCI0 = column vector containing the velocity at t = t0 %

% mu = Gravitational Parameter %

% t0 = initial time %

% t = final time %

% %

%Outputs: r = column vector containing the position at t = t %

% v = column vector containing the velocity at t = t %

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

oe0 = rv2oe_Crouse_James(rPCI0,vPCI0,mu); %calculates orbital elements

a = oe0(1);

e = oe0(2);

nu0 = oe0(6); %initial true anomaly

E0 = 2*atan2(sqrt(1-e)*sin(nu0/2),sqrt(1+e)*cos(nu0/2));%initial eccentric anomaly

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

tp = t0 - (sqrt(a^3/mu)*(E0-(e*sin(E0)))); %time of last periapsis crossing

k = floor(t-tp/tau); %number of periapsis crossings

C = (sqrt(mu/a^3)*(t-t0)) - (2*pi*k) + (E0-(e*sin(E0)));

N = 10*ceil(1/(1-e)); %nuber of iteration depends on e

E = E0 - (e*sin(E0)); %initial guess

for ii = 1:N

E(ii+1) = e*sin(E(ii)) + C;

end

E = E(N+1); %final value of E

nu = 2*atan2(sqrt(1+e)*sin(E/2),sqrt(1-e)*cos(E/2));

oe = [oe0(1) oe0(2) oe0(3) oe0(4) oe0(5) nu]';

[r,v] = oe2rv_Crouse_James(oe,mu);

end

function mercatorDisplay(lonE,lat);

earth = imread('earth.jpg');

% Open a new figure, then run the command "clf"

% Example:

% figure(<next-figure-number-here>);

% clf

% Then run the code below.

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

hold on

plot(lonE*180/pi,lat*180/pi,'w*');

function t_total = NimpulseTime(N,M,t,r0,rf,mu)

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

%inputs:

%N the number of impulse pairs

%M length of position variable during each impulse

%t timevector during impulse

%r0 radius of initial orbit

%rf radius of final orbit

%mu gravitational parameter

%-------------------------------------------------------------------------%

%output

%t_total is the total time required to complete the transfer

%-------------------------------------------------------------------------%

%variables used:

%t_tot_imp is the time during the elliptical transfer between the

% intermediate circular orbits

%t_tot_orb is the time the spacecraft spends in the intermediate

% circular orbits.

%-------------------------------------------------------------------------%

if N == 1

t_total_imp = t(M);

t_total_orb = 0;

t_total = t_total_imp + t_total_orb;

elseif N == 2

t_total_imp = t(M)+t(2*M);

t_total_orb = sqrt(((rf+r0)/2)^3/mu);

t_total = t_total_imp + t_total_orb;

elseif N == 3

t_total_imp = t(M)+t(2*M)+t(3*M);

t_total_orb = sqrt(((rf+r0)/3)^3/mu)+sqrt(((2*(rf+r0))/3)^3/mu);

t_total = t_total_imp + t_total_orb;

elseif N == 4

t_total_imp = t(M)+t(2*M)+t(3*M)+t(4*M);

t_total_orb = sqrt(((rf+r0)/4)^3/mu)+sqrt(((2*(rf+r0))/4)^3/mu)+sqrt(((3*(rf+r0))/4)^3/mu);

t_total = t_total_imp + t_total_orb;

elseif N == 5

t_total_imp = t(M)+t(2*M)+t(3*M)+t(4*M)+t(5*M);

t_total_orb = sqrt(((rf+r0)/5)^3/mu)+sqrt(((2*(rf+r0))/5)^3/mu)+sqrt(((3*(rf+r0))/5)^3/mu)+sqrt(((4*(rf+r0))/5)^3/mu);

t_total = t_total_imp + t_total_orb;

elseif N == 6

t_total_imp = t(M)+t(2*M)+t(3*M)+t(4*M)+t(5*M)+t(6*M);

t_total_orb = sqrt(((rf+r0)/6)^3/mu)+sqrt(((2*(rf+r0))/6)^3/mu)+sqrt(((3*(rf+r0))/6)^3/mu)+sqrt(((4*(rf+r0))/6)^3/mu)+sqrt(((5*(rf+r0))/6)^3/mu);

t_total = t_total_imp + t_total_orb;

end

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

end

function [rPCI,vPCI] = oe2rv_Crouse_James(oe,mu)

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

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

% This function calculates the inirtial position and velocity and %

%expresses it in the Inertial basis. The inputs are the 6 classical %

%orbital elements along with mu, the outputs are the position and velocity%

%expressed in I. %

%Inputs: a = semi-major axis %

% e = eccentricity %

% Omega = Longitude of the Ascending Node %

% inc = Orbital Inclination %

% w = Argument of the Periapsis %

% nu = True Anomaly %

% | a | %

% | e | %

% |Omega| %

% oe = | inc | %

% | w | %

% | nu | %

% | nu | %

%Outputs:rPCI = position expressed in a Planet Centered coordinate system %

% vPCI = velocity expressed in a planet centerd coordinate system %

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

% Calculations %

%-------------------------------------------------------------------------%

%defining variables

a = oe(1);

e = oe(2);

Omega = oe(3);

inc = oe(4);

w = oe(5);

nu = oe(6);

%transformations

T_PQ = [cos(w) -sin(w) 0; sin(w) cos(w) 0; 0 0 1];

T_QN = [1 0 0; 0 cos(inc) -sin(inc); 0 sin(inc) cos(inc)];

T_NI = [cos(Omega) -sin(Omega) 0; sin(Omega) cos(Omega) 0; 0 0 1];

%finding positions

r = (a*(1-e^2))/(1+e*cos(nu));

rP = [r*cos(nu) r*sin(nu) 0]'; %position in P

rPCI = T_NI*T_QN*T_PQ*rP; %inirtial position

%finding velocity

vP = sqrt(mu/(a*(1-e^2)))*[-sin(nu) e+cos(nu) 0]'; %velocity in P

vPCI = T_NI*T_QN*T_PQ*vP; %inirtial velocity

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

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

end

function [tspan,r,v] = propagateOnCircle_Crouse_James(rPCI0,vPCI0,t0,tf,mu,N)

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

% This function uses a spacecrafts position and velocity at an initial time

% and creates a position and velocity array for N = 45 points between the %

% initial time, t0, and a final time, tf for a circular orbit. %

% %

%Inputs: rPCI0 = column vector containing the position at t = t0 %

% vPCI0 = column vector containing the velocity at t = t0 %

% mu = Gravitational Parameter %

% t0 = initial time %

% tf = final time %

% %

%Outputs: r = column vector containing the position at t = t %

% v = column vector containing the velocity at t = t %

% t = columnn vector containing time series from the initial %

% time to final time spaced (tf-t0)/(N-1). %

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

%-------------------------------------------------------------------------

% Solving

r0 = norm(rPCI0);

hv = cross(rPCI0,vPCI0);

h = norm(hv);

%Setting the Units Vectors up: Scalar/Magnitude

w_1 = (rPCI0/r0)'; %basis from propblem statement

w_3 = (hv/h)';

w_2 = cross(w_3,w_1);

a = norm(rPCI0); %The radius is equal to semimajor axis for circle

tspan = (linspace(t0,tf,N))'; %Creates a vector of time from t0 to tf in N evenly spaced intervals

for ii = 1:length(tspan)

%theta and thetadot at each time interval

k = floor((tspan(ii)-t0)/((2*pi)*(sqrt(a^3/mu)))); %number of complete orbits

theta = sqrt(mu/a^3)*(tspan(ii)-t0) - (2*pi)*k;

thetadot = sqrt(mu/a^3);

%position

r(ii,:) = a*cos(theta).*w_1 + a*sin(theta).*w_2;

%velocity

v(ii,:) = -a*thetadot*sin(theta).*w_1 + a*thetadot*cos(theta).*w_2;

end

end

function oe = rv2oe_Crouse_James(rPCI,vPCI,mu)

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

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

% This function calculates the values of the six classical orbital %

% elements from the initial position and velocity vectors. The output is a%

% column vector containing 6 quantities. %

%Inputs: rv = position column vector [x y z]' %

% vv = velocity column vector [vx vy vz]' %

%Outputs: a = semi-major axis %

% e = eccentricity %

% Omega = Longitude of the Ascending Node %

% inc = Orbital Inclination %

% w = Argument of the Periapsis %

% nu = True Anomaly %

% | a | %

% | e | %

% oe = |Omega| %

% | inc | %

% | w | %

% | nu | %

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

% Calculations %

%-------------------------------------------------------------------------%

r = norm(rPCI); %r = p/(1+ecos(nu))

hv = cross(rPCI,vPCI);

h = norm(hv);

ev = (cross(vPCI,hv)/mu)-(rPCI/r);

e = norm(ev);

p = h^2/mu;

a = p/(1-e^2);

Ix = [1 0 0];

Iy = [0 1 0];

Iz = [0 0 1];

nv = cross(Iz,hv);

n = norm(nv);

%Calculating angles

%True Anomaly

hvcrossev = cross(hv,ev);

rvdot_hvcrossev = dot(rPCI,hvcrossev);

rvdotev = dot(rPCI,ev);

nu = atan2(rvdot_hvcrossev,h*rvdotev);

if nu < 0

nu = nu + (2*pi);

end

%Longitude of the Ascending Node

ndotIy = dot(nv,Iy);

ndotIx = dot(nv,Ix);

Omega = atan2(ndotIy, ndotIx);

if Omega < 0

Omega = Omega + (2*pi);

end

%Orbital Inclination

Izcrossnv = cross(Iz,nv);

hvdot_Izcrossnv = dot(hv,Izcrossnv);

hvdotIz = dot(hv,Iz);

if n == 0

inc = 0;

else

inc = atan2(-hvdot_Izcrossnv, n*hvdotIz);

end

%Argument of the Periapsis

hvcrossnv = cross(hv,nv);

evdot_hvcrossnv = dot(ev,hvcrossnv);

evdotnv = dot(ev,nv);

w = atan2(evdot_hvcrossnv, h*evdotnv);

if w < 0

w = w + (2*pi);

end

oe = [a e Omega inc w nu]'; %builds the output column vector

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

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

end

function [R,t,vv1delta,vv2delta] = twoImpulseHohmann_Crouse_James(r0,rf,inc0,incf,Omega0,Omegaf,mu)

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

% two Impulse Hohmann function

% this code solves the case were an orbital transfer needs to be performed

% and the inclination change is all done with the second impulse.

% inputs : r0 = radius of orbit 1

% rf = radius of orbit 2

% inc0 = orbital inclination of orbit 1

% incf = orbital inclination of orbit 2

% Omegao = longitude of the ascending node of orbit 1

% Omegaf = longitude of the ascending node of orbit 2

% mu = gravitational parameter

% outputs : R = Position along transfer orbit (Mx3 vector, M = 500)

% t = time corresponding to R

% vv1delta = change in velocity due to impulse 1

% vv2delta = change in velocity due to impulse 2

%-------------------------------------------------------------------------%

M = 500;

%initial orbit

a0 = r0; %semimajor axis (km)

tau0 = 2*pi*sqrt(a0^3/mu); % orhital period

torbit0_0 = 0;

torbit0_f = tau0;

e0 = 0;

w0 = 0;

nu0 = 0;

oe0 = [a0 e0 Omega0 inc0 w0 nu0]';

[rvorbit0_0,vvorbit0_0] = oe2rv_Crouse_James(oe0,mu);

[tspanorbit0,rvorbit0,vvorbit0] = propagateOnCircle_Crouse_James(rvorbit0_0,vvorbit0_0,torbit0_0,torbit0_f,mu,M);

plot3(rvorbit0(:,1),rvorbit0(:,2),rvorbit0(:,3),'k','Linewidth',2);

hold on

hv0 = cross(rvorbit0_0,vvorbit0_0); %anhular momentum vector

h0 = norm(hv0); %angular momentum

%-------------------------------------------------------------------------%

%final orbit

af = rf; %semimajor axis (km)

tauf = 2*pi*sqrt(af^3/mu); % orhital period

torbitf_0 = 0;

torbitf_f = tauf;

ef = 0;

wf = 0;

nuf = 0;

oef = [af ef Omegaf incf wf nuf]';

[rvorbitf_0,vvorbitf_0] = oe2rv_Crouse_James(oef,mu);

[tspanorbitf,rvorbitf,vvorbitf] = propagateOnCircle_Crouse_James(rvorbitf_0,vvorbitf_0,torbitf_0,torbitf_f,mu,M);

plot3(rvorbitf(:,1),rvorbitf(:,2),rvorbitf(:,3),'c','Linewidth',2);

hold on

hvf = cross(rvorbitf_0,vvorbitf_0); %anhular momentum vector

hf = norm(hv0); %angular momentum

%-------------------------------------------------------------------------%

%impulse 1

% the first impulse occurs at the periapsis of transfer orbit and thus has

% performs no inclination change

l = cross(hv0, hvf)/norm(cross(hv0, hvf)); % unit vector along line of intersection

rv0 = a0*l; % position where first impulse occurs

a1 = (a0 + af)/2; % semimajor axis of transfer orbit

e1 = (af-a0)/(af+a0); % eccentricity of transfer orbit

v1minus = sqrt(mu/a0); % speed just before first impulse

u1 = cross(hv0,l)/h0; % unit vector in the direction of the first impulse

vv1minus = v1minus*u1; % velocity just before first impulse

v1delta = sqrt((2*mu)/a0 - mu/a1) - v1minus; % speed change impulse 1

vv1delta = v1delta*u1; % velocity change impulse 1

vv1plus = vv1minus + vv1delta; % velocity after first impulse

v1plus = norm(vv1plus); % speed after first impulse

%-------------------------------------------------------------------------%

%impulse 2

% the second impulse occurs at apoapsis and performs the inclination change

v2minus = sqrt((2*mu)/af - mu/a1); %speed just before second impulse

vv2minus = -v2minus*u1; % velocity before second impulse

rv2 = -af*l; % position of second impulse

v2plus = sqrt(mu/af); % speed after second impulse

u2 = -cross(hvf,l)/hf; % direction of second impulse

vv2plus = v2plus*u2; % velocity after second impulse

vv2delta = vv2plus - vv2minus; % change in velocity impulse 2

v2delta = norm(vv2delta); %change in speed impulse 2

%-------------------------------------------------------------------------%

% calculating the position and velocity of the transfer orbit

% oe1 = rv2oe_Crouse_James(rv0,vv1plus,mu);

tau1 = 2*pi*sqrt(a1^3/mu); % half of the orbital period of the transfer orb

tspan_transfer = linspace(0,tau1/2,M)'; %time vector of transfer orbit

R(:,1) = rv0; % initial position for kepler solver

V(:,1) = vv1plus; % initial velocity for kepler solver

for ii = 2:M % kepler solver

ti = tspan_transfer(ii-1);

tf = tspan_transfer(ii);

ri = R(:,ii-1);

vi = V(:,ii-1);

[R(:,ii),V(:,ii),E,nu] = KeplerSolver_Crouse_James(ri,vi,mu,ti,tf);

end

R = R'; %Position along the transfer orbit

t = tspan_transfer; % output time vector

plot3(R(:,1),R(:,2),R(:,3),'m','Linewidth',2)

hold on

%optional quiver plots

% pp = quiver3(rv0(1),rv0(2),rv0(3),vv1minus(1)*10000,vv1minus(2)*10000,vv1minus(3)*10000);

% hold on

% pp = quiver3(rv0(1),rv0(2),rv0(3),vv1plus(1)*10000,vv1plus(2)*10000,vv1plus(3)*10000);

% hold on

% pp = quiver3(rv0(1),rv0(2),rv0(3),vv1delta(1)*10000,vv1delta(2)*10000,vv1delta(3)*10000);

% hold on

% pp = quiver3(rv2(1),rv2(2),rv2(3),vv2minus(1)*10000,vv2minus(2)*10000,vv2minus(3)*10000);

% hold on

% pp = quiver3(rv2(1),rv2(2),rv2(3),vv2plus(1)*10000,vv2plus(2)*10000,vv2plus(3)*10000);

% hold on

% pp = quiver3(rv2(1),rv2(2),rv2(3),vv2delta(1)*10000,vv2delta(2)*10000,vv2delta(3)*10000);

% hold on

legend('earth','intermediate orbit','final orbit','transfer orbit')

end

function [R,t,vdeltp,vdelta,vdelt,tr] = twoNImpulseOrbitTransfer_Crouse_James(r0,rf,inc0,incf,Omega0,Omegaf,N,mu)

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

% two N Impulse function

% this function gives the position, time, and magnitude of impulses at

% periapsis apoapsis and total, and the time required to perform an N

% impulse orbit transfer from an initial circular orbit (r0, inc0, Omega0)

% to a final circular orbit (rf, incf, Omegaf).

% inputs : r0 = radius of orbit 1

% rf = radius of orbit 2

% inc0 = orbital inclination of orbit 1

% incf = orbital inclination of orbit 2

% Omegao = longitude of the ascending node of orbit 1

% Omegaf = longitude of the ascending node of orbit 2

% N = number of impulse pairs performed

% mu = gravitational parameter

% outputs : R = Position along transfer orbit ((N*M)x3 vector, M = 500)

% t = time corresponding to R

% v1delta = change in speed due to impulse 1

% v2delta = change in speed due to impulse 2

%-------------------------------------------------------------------------%

for ii = 1:N

rows = ii + (500*(ii-1));

rowe = ii - 1 + (500*ii);

%---------------------------------------------------------------------%

%initial orbit

rn(1) = r0;

Omegan(1) = Omega0;

incn(1) = inc0;

%---------------------------------------------------------------------%

%intermediate orbit quantities

rn(ii+1) = r0 + ii/N*(rf-r0);

Omegan(ii+1) = Omega0 + ii/N*(Omegaf-Omega0);

incn(ii+1) = inc0 + ii/N*(incf-inc0);

%---------------------------------------------------------------------%

%using two impulse hohmann at each intermediate orbit

[Rn,tn,vv1deltan,vv2deltan] = twoImpulseHohmann_Crouse_James(rn(ii),rn(ii+1),incn(ii),incn(ii+1),Omegan(ii),Omegan(ii+1),mu);

%---------------------------------------------------------------------%

%saving the outputs

R(rows:rowe,:) = Rn(1:500,:);

t(rows:rowe,:) = tn(1:500,:);

vdeltp(ii) = norm(vv1deltan');

vdelta(ii) = norm(vv2deltan');

vtotdelt(ii) = vdeltp(ii) + vdelta(ii);

end

hold off

vdelt = sum(vtotdelt);

M = 500;

tr = NimpulseTime(N,M,t,r0,rf,mu);

end