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

 

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