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.

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

% selenocentric equations of motion

% includes perturbations due to:

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

% input

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

% output

%  ydot = integration vector

% Orbital Mechanics with MATLAB

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

global emu smu iearth isun jdtdb0 tmatrix2000

agrav = zeros(3);

aearth = zeros(3);

asun = zeros(3);

% current tdb julian date

jdate = jdtdb0 + t / 86400.0;

% acceleration due to the moon

agrav = gravity_pa (t, y);

if (isun == 1)
    
    % solar point-mass perturbations
    
    svsun_wrk = jplephem(jdate, 11, 10);   
    
    % convert to moon_j2000 coordinates
    
    svsun = tmatrix2000 * svsun_wrk(1:3);

    % compute heliocentric position vector of the spacecraft
    
    rs2sc = y(1:3)' - svsun(1:3);

    % f(q) formulation
        
    vtmp = y(1:3)' - 2.0d0 * svsun(1:3);

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

    dot2 = dot(svsun(1:3), svsun(1:3));

    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 * svsun(1:3)) / d3sun;
        
end

if (iearth == 1)

    % earth point-mass perturbations

    svearth_wrk = jplephem(jdate, 3, 10);
    
    % convert to moon_j2000 coordinates
    
    svearth = tmatrix2000 * svearth_wrk(1:3);
    
    % compute selenocentric position vector of the spacecraft
        
    rm2sc = y(1:3)' - svearth(1:3);

    % f(q) formulation
        
    vtmp = y(1:3)' - 2.0d0 * svearth(1:3);

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

    dot2 = dot(svearth(1:3), svearth(1:3));

    qearth = dot1 / dot2;

    fearth = qearth * ((3.0d0 + 3.0d0 * qearth + qearth * qearth) ...
        / (1.0d0 + (1.0d0 + qearth)^1.5d0));

    d3earth = norm(rm2sc)^3;

    % point-mass gravity of the earth
        
    aearth = -emu * (y(1:3)' + fearth * svearth(1:3)) / d3earth;
        
end

% compute total integration vector

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








Contact us