Code covered by the BSD License  

Highlights from
A MATLAB Script for Propagating Interplanetary Trajectories from Earth to Mars

A MATLAB Script for Propagating Interplanetary Trajectories from Earth to Mars

by

David Eagle (view profile)

 

Numerically integrate the orbital equations of motion of an Earth to Mars interplanetary trajectory.

readdata(filename)
function [fid, r_tip, v_tip, jdtdb_tip, jdtdb_user, iptype, cmfile] = readdata(filename)

% read simulation definition, gravity model and constants data files

% version for pprop_e2m.m

% Orbital Mechanics with MATLAB

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

global aunit emu req mmu xmu_mars req_mars

global isun imoon iplanet j2 j3 j4

global csrp_hel csrp_geo rmag_user

global dutet rsoi_target lgrav mgrav ccoef scoef

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!!');
    
    pause;
    
    return;
    
end

% read 96 lines of data file

for i = 1:1:96
    
    cline = fgetl(fid);
    
    switch i
        
        case 8
            
            % name of constants data file
            
            cmfile = cline;
            
        case 12
            
            % ET - UTC (seconds)
            
            dutet = str2double(cline);
            
        case 19
            
            % type of propagation
            
            iptype = str2num(cline);
            
        case 24
            
            % user-defined final epoch 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 29
            
            % user-defined final epoch UTC
            
            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 user-defined final epoch tdb julian date
            
            tdb_sec = utc_sec + dutet;
            
            day = day + utc_hr / 24.0 + utc_min / 1440.0 + tdb_sec / 86400.0;
            
            jdtdb_user = julian(month, day, year);
            
        case 32
            
            % user-defined areocentric distance (lilometers)
            
            rmag_user = str2double(cline);
            
        case 37
            
            % Earth departure 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 42
            
            % Earth departure 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 departure tdb julian date
            
            tdb_sec = utc_sec + dutet;
            
            day = day + utc_hr / 24.0 + utc_min / 1440.0 + tdb_sec / 86400.0;
            
            jdtdb_tip = julian(month, day, year);
            
        case 47
            
            % x-component of departure position vector (kilometers)
            
            r_tip(1) = str2double(cline);
            
        case 48
            
            % y-component of departure position vector (kilometers)
            
            r_tip(2) = str2double(cline);
            
        case 49
            
            % z-component of departure position vector (kilometers)
            
            r_tip(3) = str2double(cline);
            
        case 54
            
            % x-component of departure velocity vector (kilometers/second)
            
            v_tip(1) = str2double(cline);
            
        case 55
            
            % y-component of departure velocity vector (kilometers/second)
            
            v_tip(2) = str2double(cline);
            
        case 56
            
            % z-component of departure velocity vector (kilometers/second)
            
            v_tip(3) = str2double(cline);
            
        case 64
            
            % solar point-mass gravity flag
            
            isun = str2double(cline);
            
        case 68
            
            % lunar point-mass gravity flag
            
            imoon = str2double(cline);
            
        case 72
            
            % planets point-mass gravity flag
            
            iplanet = str2double(cline);
            
        case 76
            
            % name of earth gravity model data file
            
            grav_fname = cline;
            
        case 80
            
            % order of gravity model
            
            lgrav = str2double(cline);
            
        case 84
            
            % degree of gravity model
            
            mgrav = str2double(cline);
            
        case 88
            
            % spacecraft mass (kilograms)
            
            scmass = str2double(cline);
            
        case 92
            
            % srp reference area (square meters)
            
            srp_area = str2double(cline);
            
        case 96
            
            % srp reflectivity coefficient
            
            reflect = str2double(cline);
    end
end

status = fclose(fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% open and read constants and models data file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fid = fopen(cmfile, 'r');

% check for file open error

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

% read 39 lines of data file

for i = 1:1:39
    
    cline = fgetl(fid);
    
    switch i
        
        case 7
            
            % astronomical unit (kilometers)
            
            aunit = str2double(cline);
            
        case 11
            
            % speed of light (meters/second)
            
            clight = str2double(cline);
            
        case 15
            
            % solar flux at 1 AU (watts/meter**2)
            
            psr = str2double(cline);
            
        case 19
            
            % Earth gravitational constant (km**3/sec**2)
            
            emu = str2double(cline);
            
        case 23
            
            % Earth equatorial radius (kilometers)
            
            req = str2double(cline);
            
        case 27
            
            % Earth sphere-of-influence value (kilometers)
            
            rsoi_target = str2double(cline);
            
        case 31
            
            % Moon gravitational constant (km**3/sec**2)
            
            mmu = str2double(cline);
            
        case 35
            
            % Mars gravitational constant (km**3/sec**2)
            
            xmu_mars = str2double(cline);
            
        case 39
            
            % Mars equatorial radius (kilometers)
            
            req_mars = str2double(cline);
    end
    
end

status = fclose(fid);

% read earth gravity model data file

[ccoef, scoef] = readgm(grav_fname);

j2 = -ccoef(3, 1);

j3 = -ccoef(4, 1);

j4 = -ccoef(5, 1);

% solar flux at 1 au (N/m**2)

ps = psr / clight;

% conversion factor (meters/second**2 to au/day**2)

xmult = 86400.0d0^2 / (1000.0d0 * aunit);

% heliocentric solar radiation constant (au/day^2)

if (scmass ~= 0.0d0)
    
    csrp_hel = xmult * reflect * ps * srp_area / scmass;
    
else
    
    csrp_hel = 0.0;
    
end

% geocentric solar radiation constant (meters/second^2)

if (scmass ~= 0.0d0)
    
    csrp_geo = reflect * ps * srp_area / scmass;
    
else
    
    csrp_geo = 0.0;
    
end

Contact us