Code covered by the BSD License  

Highlights from
A MATLAB Script for Predicting the Evolution of Lunar Orbits

A MATLAB Script for Predicting the Evolution of Lunar Orbits

by

 

Script for propagating lunar orbits subject to non-spherical lunar gravity and third-body gravity.

gravity_pa (t, sv)
function agrav = gravity_pa (t, sv)

% first order equations of selenocentric orbital motion

% input

%  t  = mission elapsed time (seconds)
%  sv = state vector

% output

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

% Orbital Mechanics with MATLAB

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

global mmu req jdtdb0

global lgrav mgrav ccoef scoef

% initialize acceleration vector

agrav = zeros(3, 1);

awrk = zeros(3, 1);

rwrk1 = zeros(3, 1);

% load "working" position vector

rwrk1 = sv(1:3);

% current tdb julian date

jdate = jdtdb0 + t / 86400.0d0;

% compute moon_j2000 to moon_pa transformation matrix

tmatrix1 = moon_pa(jdate);

% convert position vector to moon_pa coordinates

rwrk = tmatrix1 * rwrk1';

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

r1 = sqrt(r2);

r3 = r2 * r1;

% user-defined degree and order gravity model

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

sr1 = sqrt(sr2);

sphi = rwrk(3) / r1;

phi = asin(sphi);

% east longitude of the spacecraft

lamda = mod(atan3(rwrk(2), rwrk(1)), 2.0 * pi);

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:1: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 = mmu / r1;

e1 = -e1 * d3 / r1;

e2 = e2 * d3;

e3 = e3 * d3;

d1 = e1 / r1;

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

d3 = e3 / sr2;

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

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

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

% convert gravity acceleration vector back to moon_j2000 frame

agrav = tmatrix1' * awrk;

% compute gravity acceleration vector

% (note: keplerian contribution uses moon_j2000 position vector)

for i = 1:1:3
    
    agrav(i) = agrav(i) - mmu * sv(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

Contact us