Code covered by the BSD License  

Highlights from
Optimal Impulsive Orbital Transfer

Optimal Impulsive Orbital Transfer

by

 

27 Nov 2012 (Updated )

MATLAB script for the solution of the one and two impulse orbit transfer between two Earth orbits.

oota_iniz
function oota_iniz

% initialization routine

% required by oota_matlab.m

% Orbital Mechanics with MATLAB

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global rinc

global sma ecc inc argper raan slr

global apeci1 apeci2 teci2op

global aapop1 aapop2 wop1 wop2 ecc1 ecc2

% semiparameters of initial and final orbits

slr(1) = sma(1) * (1.0 - ecc(1) * ecc(1));

slr(2) = sma(2) * (1.0 - ecc(2) * ecc(2));

% -----------------------------------------
% compute eci unit angular momentum vectors
% -----------------------------------------

% initial orbit

weci1(1) = sin(raan(1)) * sin(inc(1));

weci1(2) = -cos(raan(1)) * sin(inc(1));

weci1(3) = cos(inc(1));

% final orbit

weci2(1) = sin(raan(2)) * sin(inc(2));

weci2(2) = -cos(raan(2)) * sin(inc(2));

weci2(3) = cos(inc(2));

% compute relative inclination

crinc = dot(weci1, weci2);

rinc = acos(crinc);

% -------------------------------------------------
% compute orbit plane unit angular momentum vectors
% -------------------------------------------------

% initial orbit

wop1(1) = 0.0;

wop1(2) = -sin(rinc);

wop1(3) = cos(rinc);

% final orbit

wop2(1) = 0.0;

wop2(2) = 0.0;

wop2(3) = 1.0;

% ---------------------------
% compute eci perigee vectors
% ---------------------------

apeci1(1) = cos(argper(1)) * cos(raan(1)) ...
    - sin(argper(1)) * sin(raan(1)) * cos(inc(1));

apeci1(2) = cos(argper(1)) * sin(raan(1)) ...
    + sin(argper(1)) * cos(raan(1)) * cos(inc(1));

apeci1(3) = sin(argper(1)) * sin(inc(1));

apeci2(1) = cos(argper(2)) * cos(raan(2)) ...
    - sin(argper(2)) * sin(raan(2)) * cos(inc(2));

apeci2(2) = cos(argper(2)) * sin(raan(2)) ...
    + sin(argper(2)) * cos(raan(2)) * cos(inc(2));

apeci2(3) = sin(argper(2)) * sin(inc(2));

% compute reference coordinate system-to-eci transformation matrices

tmatrices;

% compute reference coordinate system argument of perigee

apop1 = teci2op * apeci1';

apop2 = teci2op * apeci2';

% first impulse

aapop1 = acos(apop1(1));

if (apop1(3) < 0.0)
    
    aapop1 = -aapop1;
    
end

% second impulse

aapop2 = atan3(apop2(2), apop2(1));

% --------------------------------------------------------
% compute reference coordinate system eccentricity vectors
% --------------------------------------------------------

% initial orbit

ecc1(1) = ecc(1) * cos(aapop1);

ecc1(2) = ecc(1) * sin(aapop1) * cos(rinc);

ecc1(3) = ecc(1) * sin(aapop1) * sin(rinc);

% final orbit

ecc2(1) = ecc(2) * cos(aapop2);

ecc2(2) = ecc(2) * sin(aapop2);

ecc2(3) = 0.0;

%%%%%%%%%%%%%%%%%%

function tmatrices

% reference coordinate system to/from eci transformation matrices

% Orbital Mechanics with MATLAB

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global inc raan

global teci2op top2eci

% eci u vector

ueci(1) = cos(raan(2));

ueci(2) = sin(raan(2));

ueci(3) = 0.0;

% eci w vector - first impulse

weci1(1) = sin(raan(1)) * sin(inc(1));

weci1(2) = -cos(raan(1)) * sin(inc(1));

weci1(3) = cos(inc(1));

% eci w vector - second impulse

weci2(1) = sin(raan(2)) * sin(inc(2));

weci2(2) = -cos(raan(2)) * sin(inc(2));

weci2(3) = cos(inc(2));

% eci n vector

xneci = cross(weci2, weci1);

xnecim = sqrt(dot(xneci, xneci));

% eci n unit vector

if (xnecim == 0.0)
    
    xneci(1) = 1.0;
    
    xneci(2) = 0.0;
    
    xneci(3) = 0.0;
    
else
    
    xneci = xneci / xnecim;
    
end

xndotu = dot(xneci, ueci);

% argument of hinge line

xphi = acos(xndotu);

if (xneci(3) < 0.0)
    
    xphi = -xphi;
    
end

% reference coordinate system-to-eci transformation matrix

top2eci(1, 1) = cos(raan(2)) * cos(xphi) ...
    - sin(raan(2)) * cos(inc(2)) * sin(xphi);

top2eci(1, 2) = -cos(raan(2)) * sin(xphi) ...
    - sin(raan(2)) * cos(inc(2)) * cos(xphi);

top2eci(1, 3) = sin(raan(2)) * sin(inc(2));

top2eci(2, 1) = sin(raan(2)) * cos(xphi) ...
    + cos(raan(2)) * cos(inc(2)) * sin(xphi);

top2eci(2, 2) = -sin(raan(2)) * sin(xphi) ...
    + cos(raan(2)) * cos(inc(2)) * cos(xphi);

top2eci(2, 3) = -cos(raan(2)) * sin(inc(2));

top2eci(3, 1) = sin(inc(2)) * sin(xphi);

top2eci(3, 2) = sin(inc(2)) * cos(xphi);

top2eci(3, 3) = cos(inc(2));

% eci-to-reference coordinate system transformation matrix

teci2op = top2eci';

Contact us