Code covered by the BSD License

# A MATLAB Script for Propagating Interplanetary Trajectories from Earth to Mars

### David Eagle (view profile)

Numerically integrate the orbital equations of motion of an Earth to Mars interplanetary trajectory.

geo_eqm (t, y)
function ydot = geo_eqm (t, y)

% geocentric equations of motion

% version for tcm_e2m.m

% input

%  t = simulation time (seconds)
%  y = state vector (kilometers and kilometers/second)

% output

%  ydot = integration vector

% Orbital Mechanics with MATLAB

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

global mmu smu xmu aunit req

global imoon isun iplanet

global csrp_geo jdtdb_wrk

agrav = zeros(3, 1);

amoon = zeros(3, 1);

asun = zeros(3, 1);

asrp = zeros(3, 1);

accp = zeros(3, 1);

% current tdb julian date

jdate = jdtdb_wrk + t / 86400.0;

% acceleration due to the earth

agrav = gravity (t, y);

if (isun == 1)

%%%%%%%%%%%%%%%%%%%%%
% solar perturbations
%%%%%%%%%%%%%%%%%%%%%

svsun = jplephem(jdate, 11, 3);

rsun = aunit * svsun(1:3);

% compute heliocentric position vector of the spacecraft

rs2sc = y(1:3) - rsun(1:3);

% f(q) formulation

vtmp = y(1:3) - 2.0d0 * rsun(1:3);

dot1 = dot(y(1:3), vtmp);

dot2 = dot(rsun, rsun);

qsun = dot1 / dot2;

fsun = qsun * ((3.0d0 + 3.0d0 * qsun + qsun * qsun) ...
/ (1.0d0 + (1.0d0 + qsun)^1.5d0));

d3sun = norm(rs2sc)^3;

% point-mass gravity of the sun

asun = -smu * (y(1:3) + fsun * rsun(1:3)) / d3sun;

end

if (imoon == 1)

%%%%%%%%%%%%%%%%%%%%%
% lunar perturbations
%%%%%%%%%%%%%%%%%%%%%

svmoon = jplephem(jdate, 10, 3);

rmoon = aunit * svmoon(1:3);

% compute selenocentric position vector of the spacecraft

rm2sc = y(1:3) - rmoon;

% f(q) formulation

vtmp = y(1:3) - 2.0d0 * rmoon;

dot1 = dot(y(1:3), vtmp);

dot2 = dot(rmoon, rmoon);

qmoon = dot1 / dot2;

fmoon = qmoon * ((3.0d0 + 3.0d0 * qmoon + qmoon * qmoon) ...
/ (1.0d0 + (1.0d0 + qmoon)^1.5d0));

d3moon = norm(rm2sc)^3;

% point-mass gravity of the moon

amoon = -mmu * (y(1:3) + fmoon * rmoon(1:3)) / d3moon;

end

if (csrp_geo > 0.0d0)

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

svsun = jplephem(jdate, 11, 3);

rsun = aunit * svsun(1:3);

% geocentric distance of the sun (kilometers)

rmsun = norm(rsun);

% geocentric unit vector of the sun

usun = svsun(1:3) / rmsun;

% geocentric distance of the spacecraft

rscm = norm(y(1:3));

% compute unit position vector of spacecraft

usat = y(1:3) / rscm;

a = usat(2) * usun(3) - usat(3) * usun(2);
b = usat(3) * usun(1) - usat(1) * usun(3);
c = usat(1) * usun(2) - usat(2) * usun(1);

d = sqrt(a * a + b * b + c * c);

e = dot(usat, usun);

u = asin(0.00460983743 / (rmsun / aunit));

p = asin(0.0046951089 / (rmsun / aunit));

if (e > 0.0)

q = -d;

else

q = d;

end

% increase the Earth's radius by 90 km
% to account for the atmosphere

ratm = req + 90.0;

b = asin(ratm / rscm);

v = b - u;
w = b + p;

xx = sin(v);
yy = sin(w);

if (q <= yy && q > xx)

% penumbra

csrp = 0.0;

elseif (q <= xx && q >= 0)

% umbra

csrp = 0.0;

else

% sunlight

csrp = csrp_geo;

end

% compute heliocentric position vector and magnitude of the spacecraft

rs2sc = y(1:3) - rsun;

rs2scm = norm(rs2sc);

rs2scm3 = rs2scm * rs2scm * rs2scm;

% compute acceleration vector (km / sec^2)

asrp = 0.001d0 * csrp * aunit^2 * rs2sc / rs2scm3;

end

if (iplanet == 1)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% planetary point-mass perturbations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i = 1:1:9

svplanet = jplephem(jdate, i, 11);

rplanet = svplanet(1:3);

for j = 1:1:3

rp(i, j) = aunit * rplanet(j);

end

end

% compute planet-centered position vectors of spacecraft (au)

for i = 1:1:9

rp2sc(i, 1) = y(1) - rp(i, 1);

rp2sc(i, 2) = y(2) - rp(i, 2);

rp2sc(i, 3) = y(3) - rp(i, 3);

end

% compute f(q) functions for each planet

for k = 1:1:9

qp(k) = dot(y(1:3), y(1:3) - 2.0 * rp(k, :)') / dot(rp(k, :), rp(k, :));

fp(k) = qp(k) * ((3.0 + 3.0 * qp(k) + qp(k) * qp(k)) / (1.0 + (1.0 + qp(k))^1.5));

d3p(k) = norm(rp2sc(k, :)) * norm(rp2sc(k, :)) * norm(rp2sc(k, :));

end

% compute planetary perturbations

for j = 1:1:3

accp(j) = 0.0;

for k = 1:1:9

xmuwrk = aunit^3 * xmu(k + 1) / 86400.0d0 ^2;

% note: exclude Earth perturbation

if (k ~= 3)

accp(j) = accp(j) - xmuwrk * (y(j) + fp(k) * rp(k, j)) / d3p(k);

end

end

end

end

% compute total integration vector

ydot = [ y(4)
y(5)
y(6)
agrav(1) + amoon(1) + asun(1) + accp(1) + asrp(1)
agrav(2) + amoon(2) + asun(2) + accp(2) + asrp(2)
agrav(3) + amoon(3) + asun(3) + accp(3) + asrp(3)];