Code covered by the BSD License

# Cowell's Method - MICE Version

### David Eagle (view profile)

PDF document and MATLAB script that demonstrates using Cowellâ€™s method to predict orbital motion.

gravity_mice (t, y)
```function agrav = gravity_mice (t, y)

% N degree and M order
% gravitational acceleration

% input

%  t = simulation time (seconds)
%  y = state vector

% output

%  agrav = ECI gravitational acceleration
%          vector (km/sec/sec)

% Orbital Mechanics with Matlab

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

global emu req j2 j4 etime0

global lgrav mgrav ccoef scoef

% initialize acceleration vector

agrav = zeros(6);

awrk = zeros(6);

% current ephemeris elapsed seconds since J2000

etime = etime0 + t;

% compute true-of-date state vector

xform = cspice_sxform('J2000', 'EARTH_TOD', etime);

ytod = xform * y';

r2 = ytod(1) * ytod(1) + ytod(2) * ytod(2) + ytod(3) * ytod(3);

r1 = sqrt(r2);

r3 = r2 * r1;

if (lgrav == 0 && mgrav == 0)

%%%%%%%%%%%%%%%%%%
% Keplerian motion
%%%%%%%%%%%%%%%%%%

elseif (lgrav == 2 && mgrav == 0)

%%%%%%%%%
% j2 only
%%%%%%%%%

r5 = r2 * r3;

d1 = -1.5 * j2 * req * req * emu / r5;

d2 = 1 - 5 * ytod(3) * ytod(3) / r2;

awrk(1) = awrk(1) + ytod(1) * d1 * d2;

awrk(2) = awrk(2) + ytod(2) * d1 * d2;

awrk(3) = awrk(3) + ytod(3) * d1 * (d2 + 2);

xform = cspice_sxform('EARTH_TOD', 'J2000', etime);

agrav = xform * awrk';

elseif (lgrav == 4 && mgrav == 0)

%%%%%%%%%%%%%%%%
% j2 and j4 only
%%%%%%%%%%%%%%%%

r5 = r2 * r3;

d1 = -1.5d0 * j2 * req * req * emu / r5;

d2 = 1.0d0 - 5.0d0 * ytod(3) * ytod(3) / r2;

aj2(1) = ytod(1) * d1 * d2;

aj2(2) = ytod(2) * d1 * d2;

aj2(3) = ytod(3) * d1 * (d2 + 2.0d0);

r4 = r2 * r2;

r7 = r2 * r5;

r9 = r2 * r7;

d3 = 15.0d0 * j4 * req^4 * emu / (8.0d0 * r7);

aj4(1) = ytod(1) * d3 * (1.0d0 - 14.0d0 * ytod(3) * ytod(3) ...
/ r2 + 21.0d0 * ytod(3)^4 / r4);

aj4(2) = ytod(2) * d3 * (1.0d0 - 14.0d0 * ytod(3) * ytod(3) ...
/ r2 + 21.0d0 * ytod(3)^4 / r4);

aj4(3) = ytod(3) * d3 * (5.0d0 - 70.0d0 * ytod(3) * ytod(3) ...
/ (3.0d0 * r2) + 21.0d0 * ytod(3)^4 / r4);

for i = 1:1:3
awrk(i) = aj2(i) + aj4(i);
end

xform = cspice_sxform('EARTH_TOD', 'J2000', etime);

agrav = xform * awrk';

else

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% user-defined degree and order gravity model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% compute ecf state vector

xform = cspice_sxform('J2000', 'IAU_EARTH', etime);

yecf = xform * y';

ires = 1;

sr2 = ytod(1) * ytod(1) + ytod(2) * ytod(2);

sr1 = sqrt(sr2);

sphi = yecf(3) / r1;

phi = asin(sphi);

lamda = atan3(yecf(2), yecf(1));

im = mgrav;

if (mgrav < lgrav)
im = mgrav + 1;
end

p = legend(lgrav, im, sphi);

[cn, sn, tn] = angles(mgrav, lamda, phi);

e1 = 0;
e2 = 0;
e3 = 0;

d1 = req / r1;
d2 = 1;

for il = 1:1:lgrav
f1 = 0;
f2 = 0;
f3 = 0;

il1 = il + 1;

for im = 0:il
if (im > mgrav)
break;
end

im2 = im + 2;
im1 = im + 1;

d3 = ccoef(il1, im1) * cn(im1) + scoef(il1, im1) * sn(im1);

f1 = f1 + d3 * p(il1, im1);

if (im2 <= il1)
f2 = f2 + d3 * (p(il1, im2) - tn(im1) * p(il1, im1));
end

if (im2 > il1)
f2 = f2 - d3 * tn(im1) * p(il1, im1);
end

if (im ~= 0)
f3 = f3 + im * (scoef(il1, im1) * cn(im1) - ccoef(il1, im1) ...
* sn(im1)) * p(il1, im1);
end
end

d2 = d2 * d1;

e1 = e1 + d2 * il1 * f1;

e2 = e2 + d2 * f2;

e3 = e3 + d2 * f3;
end

d3 = emu / r1;

e1 = -e1 * d3 / r1;

e2 = e2 * d3;

e3 = e3 * d3;

d1 = e1 / r1;

d2 = ytod(3) / (r2 * sr1) * e2;

d3 = e3 / sr2;

awrk(1) = awrk(1) + (d1 - d2) * ytod(1) - d3 * ytod(2);

awrk(2) = awrk(2) + (d1 - d2) * ytod(2) + d3 * ytod(1);

awrk(3) = awrk(3) + d1 * ytod(3) + sr1 * e2 / r2;

% compute gravity vector in eme2000 frame

xform = cspice_sxform('EARTH_TOD', 'J2000', etime);

agrav = xform * awrk';
end

% complete gravity acceleration vector
% (note: keplerian contribution uses j2000 position vector)

for i = 1:1:3
agrav(i) = agrav(i) - emu * y(i) / r3;
end

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

function [cn, sn, tn] = angles (m, a, b)

% multiple angles

cn = zeros(21, 1);
sn = zeros(21, 1);
tn = zeros(21, 1);

cn(1) = 1;
sn(1) = 0;
tn(1) = 0;

if (m == 0)
return;
end

cn(2) = cos(a);
sn(2) = sin(a);
tn(2) = tan(b);

if (m == 1)
return;
end

for i = 2:1:m
cn(i + 1) = 2 * cn(2) * cn(i) - cn(i - 1);
sn(i + 1) = 2 * cn(2) * sn(i) - sn(i - 1);
tn(i + 1) = tn(2) + tn(i);
end

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

function p = legend (n, m, x)

% Legendre polynomials

p(1, 1) = 1;
p(2, 1) = x;

for i = 2:1:n
p(i + 1, 1) = ((2 * i - 1) * x * p(i, 1) -(i - 1) * p(i - 1, 1)) / i;
end

if (m == 0)
return;
end

y = sqrt(1 - x * x);

p(2, 2) = y;

if (m == 1)
% null
else
for i = 2:1:m
p(i + 1, i + 1) = (2 * i - 1) * y * p(i, i);
end
end

for i = 2:1:n

i1 = i - 1;

for j = 1:1:i1

if (j > m)
break;
end

p(i + 1, j + 1) = (2 * i - 1) * y * p(i, j);

if (i - 2 >= j)
p(i + 1, j + 1) = p(i + 1, j + 1) + p(i - 1, j + 1);
end
end
end
```