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.

lprop_matlab.m
% lprop_matlab.m      February 21, 2013

% long-term evolution of lunar orbits

% includes perturbations due to:

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

% Orbital Mechanics with MATLAB

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

clear all;

global dtr rtd tmatrix2 tmatrix2000

global emu smu mmu aunit omega

global lgrav mgrav req

global jdtdb0 isun iearth

global rkcoef ccoef scoef

global iephem km ephname

% moon_j2000-to-j2000 transformation matrix

tmatrix2(1, 1) = 0.998496505205088;
tmatrix2(1, 2) = 4.993572939853833E-002;
tmatrix2(1, 3) = -2.260867140418499E-002;

tmatrix2(2, 1) = -5.481540926807404E-002;
tmatrix2(2, 2) = 0.909610125238044;
tmatrix2(2, 3) = -0.411830900942612;

tmatrix2(3, 1) = 0.000000000000000E+000;
tmatrix2(3, 2) = 0.412451018902688;
tmatrix2(3, 3) = 0.910979778593430;

% lunar mean equator and IAU node of j2000 (moon_j2000)
% to Earth mean equator and equinox of j2000 (EME2000)
% transformation matrix

tmatrix2000(1, 1) = 0.998496505205088d0;
tmatrix2000(1, 2) = -5.481540926807404E-002;
tmatrix2000(1, 3) = 0.000000000000000E+000;

tmatrix2000(2, 1) = 4.993572939853833E-002;
tmatrix2000(2, 2) = 0.909610125238044;
tmatrix2000(2, 3) = 0.412451018902688;

tmatrix2000(3, 1) = -2.260867140418499E-002;
tmatrix2000(3, 2) = -0.411830900942612;
tmatrix2000(3, 3) = 0.910979778593430;

% read astrodynamic and utility constants

om_constants;

% angular conversion constants

dtr = pi / 180.0;

rtd = 180.0 / pi;

% read leaps seconds data file

readleap;

% initialize JPL ephemeris

iephem = 1;

km = 1;

ephname = 'de421.bin';

% initialize rkf78 integrator

rkcoef = 1;

% define number of differential equations

neq = 6;

% begin simulation

clc; home;

fprintf('\nprogram lprop_matlab\n');

fprintf('\n< long-term evolution of lunar orbits >\n\n');

% request simulation data file

[filename, pathname] = uigetfile('*.in', 'please select an input data file');

[fid, jdutc0, oevi, tdur, lgrav mgrav, iearth, isun, gravfname, ...
    tetol, outputfname, dtcsv] = readdata(filename);

% compute initial tdb julian date

jdtdb0 = utc2tdb (jdutc0);

% read lunar gravity model data file

[ccoef, scoef] = readgm(gravfname);

% determine initial eci state vector

[ri, vi] = orb2eci(mmu, oevi);

% create initial data

ti = 0;

% create ascii csv data file

fid = fopen(outputfname, 'w');

% write initial mission elapsed time and orbital elements to data file

rp = oevi(1) * (1.0d0 - oevi(2));

ra = oevi(1) * (1.0d0 + oevi(2));

hp = rp - req;

ha = ra - req;

fprintf(fid, '%+16.14e, %+16.14e, %+16.14e, %+16.14e, %+16.14e, %+16.14e, %+16.14e, %+16.14e, %+16.14e\n', ...
    ti / 86400.0, oevi(1), oevi(2), rtd * oevi(3), rtd * oevi(4), rtd * oevi(5), rtd * oevi(6), hp, ha);

% load initial position and velocity vectors

for i = 1:1:3
    
    yi(i) = ri(i);
    
    yi(i + 3) = vi(i);
    
end

% print initial conditions

[cdstr0, utstr0] = jd2str(jdutc0);

fprintf('\ninitial UTC calendar date   ');

disp(cdstr0);

fprintf('\ninitial UTC time            ');

disp(utstr0);

[cdstr0, utstr0] = jd2str(jdtdb0);

fprintf('\ninitial TDB calendar date   ');

disp(cdstr0);

fprintf('\ninitial TDB time            ');

disp(utstr0);

fprintf('\ninitial orbital elements and state vector');
fprintf('\n(lunar mean equator and IAU node of J2000)');
fprintf('\n------------------------------------------\n');

oeprint1(mmu, oevi, 1);

svprint(ri, vi);

% initialize initial and final times and total duration

ti = 0.0d0;

tf = 0.0d0;

tsim = 86400.0d0 * tdur;

fprintf('\nplease wait, computing data ...\n\n');

while(1)
    
    % step size guess (seconds)
    
    h = 10.0;
    
    ti = tf;
    
    tf = ti + 60.0d0 * dtcsv;
    
    % integrate from ti to tf
    
    yfinal = rkf78('sel_eqm', neq, ti, tf, h, tetol, yi);
    
    % compute current state vector
    
    for i = 1:1:3
        
        rf(i) = yfinal(i);
        
        vf(i) = yfinal(i + 3);
        
    end
    
    % compute current orbital elements
    
    oevf = eci2orb1(mmu, rf, vf);
    
    % write mission elapsed time and orbital elements to data file
    
    rp = oevf(1) * (1.0d0 - oevf(2));
    
    ra = oevf(1) * (1.0d0 + oevf(2));
    
    hp = rp - req;
    
    ha = ra - req;
    
    fprintf(fid, '%+16.14e, %+16.14e, %+16.14e, %+16.14e, %+16.14e, %+16.14e, %+16.14e, %+16.14e, %+16.14e\n', ...
        tf / 86400.0, oevf(1), oevf(2), rtd * oevf(3), rtd * oevf(4), rtd * oevf(5), rtd * oevf(6), hp, ha);
    
    yi = yfinal;
    
    % check for end of simulation
    
    if (tf >= tsim)
        
        break;
        
    end
    
end

% compute final state vector and orbital elements

for i = 1:1:3
    
    rf(i) = yfinal(i);
    
    vf(i) = yfinal(i + 3);
    
end

oevf = eci2orb1(mmu, rf, vf);

% print results

jdutcf = jdutc0 + tf / 86400;

[cdstrf, utstrf] = jd2str(jdutcf);

fprintf('\nfinal UTC calendar date     ');

disp(cdstrf);

fprintf('\nfinal UTC time              ');

disp(utstrf);

jdtdbf = jdtdb0 + tf / 86400;

[cdstrf, utstrf] = jd2str(jdtdbf);

fprintf('\nfinal TDB calendar date     ');

disp(cdstrf);

fprintf('\nfinal TDB time              ');

disp(utstrf);

fprintf('\nfinal orbital elements and state vector');
fprintf('\n(lunar mean equator and IAU node of J2000)');
fprintf('\n------------------------------------------\n');

oeprint1(mmu, oevf, 1);

svprint(rf, vf);

fprintf('\ndegree of gravity model    %2i', lgrav);

fprintf('\norder of gravity model     %2i \n', mgrav);

fprintf('\n\n');

fclose('all');

Contact us