Code covered by the BSD License

# The Long-term Evolution of Geosynchronous Transfer Orbits

### David Eagle (view profile)

Interactive MATLAB script that predicts the long-term evolution of geosynchronous transfer orbits.

ceqm1 (t, y)
```function ydot = ceqm1 (t, y)

% first order form of Cowell's equations of orbital motion

% includes perturbations due to:

%   non-spherical earth gravity
%   point mass sun and moon gravity

% input

%  t = current simulation time (seconds)
%  y = current eci state vector (km and km/sec)

% output

%  ydot = eci acceleration vector (km/sec^2)

% Orbital Mechanics with MATLAB

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

global isun imoon jdate0 smu mmu

agrav = zeros(3, 1);

asun = zeros(3, 1);

amoon = zeros(3, 1);

ydot = zeros(1, 6);

% compute geopotential perturbation

agrav = gravity(t, y);

if (isun == 1)

% solar perturbations

jdate = jdate0 + t / 86400;

[rasc, decl, rsun] = sun1(jdate);

% compute heliocentric position vector of the spacecraft

for i = 1:3

rs2sc(i) = y(i) - rsun(i);

end

% f(q) formulation

for i = 1:1:3

vtmp(i) = y(i) - 2.0d0 * rsun(i);

end

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

for i = 1:1:3

asun(i) = -smu * (y(i) + fsun * rsun(i)) / d3sun;

end

end

if (imoon == 1)

% lunar perturbations

jdate = jdate0 + t / 86400;

[rasc, decl, rmoon] = moon(jdate);

% compute selenocentric position vector of the spacecraft

for i = 1:1:3

rm2sc(i) = y(i) - rmoon(i);

end

% f(q) formulation

for i = 1:1:3

vtmp(i) = y(i) - 2.0d0 * rmoon(i);

end

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

for i = 1:1:3

amoon(i) = -mmu * (y(i) + fmoon * rmoon(i)) / d3moon;

end

end

% total acceleration vector

ydot(1) = y(4);
ydot(2) = y(5);
ydot(3) = y(6);

for i = 1:1:3

ydot(i + 3) = agrav(i) + asun(i) + amoon(i);

end

```