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.

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

% read data file for lprop_matlab.m

% NOTE: all angular elements are returned in radians

% Orbital Mechanics with MATLAB

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

dtr = pi / 180.0;

% open data file

fid = fopen(filename, 'r');

% check for file open error

if (fid == -1)
    
    clc; home;
    
    fprintf('\n\n error: cannot find this file!!');
    
    keycheck;
    
    return;
    
end

% read 85 lines of data file

for i = 1:1:85
    
    cline = fgetl(fid);
    
    switch i
        
        case 13
            
            % initial UTC calendar date
            
            tl = size(cline);
            
            ci = findstr(cline, ',');
            
            % extract month, day and year
            
            month = str2double(cline(1:ci(1)-1));
            
            day = str2double(cline(ci(1)+1:ci(2)-1));
            
            year = str2double(cline(ci(2)+1:tl(2)));
            
        case 17
            
            % initial UTC epoch
            
            utstr = cline;
            
            tl = size(utstr);
            
            ci = findstr(utstr, ',');
            
            % extract hours, minutes and seconds
            
            utc_hr = str2num(utstr(1:ci(1)-1));
            
            utc_min = str2num(utstr(ci(1)+1:ci(2)-1));
            
            utc_sec = str2num(utstr(ci(2)+1:tl(2)));
            
            % compute utc julian date
            
            day = day + utc_hr / 24.0 + utc_min / 1440.0 + utc_sec / 86400.0;
            
            jdutc0 = julian(month, day, year);
            
        case 21
            
            % semimajor axis (kilometers)
            
            oev1(1) = str2double(cline);
            
        case 25
            
            % orbital eccentricity (non-dimensional)
            
            oev1(2) = str2double(cline);
            
        case 29
            
            % orbital inclination (radians)
            
            oev1(3) = dtr * str2double(cline);
            
        case 33
            
            % argument of periapsis (radians)
            
            oev1(4) = dtr * str2double(cline);
            
        case 37
            
            % raan (radians)
            
            oev1(5) = dtr * str2double(cline);
            
        case 41
            
            % true anomaly (radians)
            
            oev1(6) = dtr * str2double(cline);
            
        case 45
            
            % propagation duration (days)
            
            tdur = str2double(cline);
            
        case 53
            
            % order of gravity model
            
            lgrav = str2double(cline);
            
        case 57
            
            % degree of gravity model
            
            mgrav = str2double(cline);
            
        case 61
            
            % earth perturbation
            
            iearth = str2double(cline);
            
        case 65
            
            % solar perturbation
            
            isun = str2double(cline);
            
        case 73
            
            % truncation error tolerance
            
            tetol = str2double(cline);
            
        case 77
            
            % name of lunar gravity model data file
            
            gravfname = cline;
            
        case 81
            
            % name of output data file
            
            outputfname = cline;
            
        case 85
            
            % data file step size (minutes)
            
            dtcsv = str2double(cline);
            
    end
    
end

status = fclose(fid);

Contact us