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.

pprop_e2m.m
% pprop_e2m.m        April 15, 2013

% propagation of Earth-to-Mars interplanetary trajectories

% user-defined final condition choices

%  (1) Mars closest approach
%  (2) user-defined areocentric distance
%  (3) user-defined final epoch

% JPL ephemeris and NOVAS sidereal time

% Orbital Mechanics with MATLAB

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

clear all;

global aunit iephem km ephname dtr rtd

global emu smu xmu xmu_mars

global jdtdb_wrk

% ecliptic-to-equatorial transformation matrix

eq2000(1, 1) = 1.0d0;
eq2000(1, 2) = 0.0d0;
eq2000(1, 3) = 0.0d0;

eq2000(2, 1) = 0.0d0;
eq2000(2, 2) = 0.917482062069182d0;
eq2000(2, 3) = -0.397777155931914d0;

eq2000(3, 1) = 0.0d0;
eq2000(3, 2) = 0.397777155931914d0;
eq2000(3, 3) = 0.917482062069182d0;

% compute equatorial-to-ecliptic transformation matrix

eqwrk = eq2000';

% angular conversion factors

rtd = 180.0d0 / pi;

dtr = pi / 180.0d0;

% initialize JPL ephemeris

ephname = 'de421.bin';

iephem = 1;

km = 0;

% **************************************
% gravitational constants (au**3/day**2)
% DE 421 values
% **************************************

% sun

xmu(1) = 0.295912208285591100d-3;

% mercury

xmu(2) = 0.491254957186794000D-10;

% venus

xmu(3) = 0.724345233269844100D-09;

% earth - computed from user-defined value

% mars

xmu(5) = 0.954954869562239000D-10;

% jupiter

xmu(6) = 0.282534584085505000D-06;

% saturn

xmu(7) = 0.845970607330847800D-07;

% uranus

xmu(8) = 0.129202482579265000D-07;

% neptune

xmu(9) = 0.152435910924974000D-07;

% pluto

xmu(10) = 0.217844105199052000D-11;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% read simulation definition data file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[filename, pathname] = uigetfile('*.in', 'Please select the input file to read');

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

% sun gravitational constant (km**3/sec**2)

smu = xmu(1) * aunit^3 / 86400.0d0^2;

% earth gravitational constant (au**3/day**2)

xmu(4) = emu * 86400.0d0^2 / aunit^3;

% begin simulation

clc; home;

fprintf('\npprop_e2m - Earth-to-Mars trajectory propagation');
fprintf('\n------------------------------------------------\n');

fprintf('\nsimulation definition data file ==> ');

disp(filename);

fprintf('\nconstants data file ==> ');

disp(cmfile);

if (iptype == 1)
  
    fprintf('\nMars closest approach\n');
    
end

if (iptype == 2)
   
    fprintf('\nuser-defined areocentric distance\n');
    
end

if (iptype == 3)
   
    fprintf('\nuser-defined final epoch\n');
    
end

fprintf('\nDE421 ephemeris\n');

fprintf('\n\ndeparture hyperbola characteristics');
fprintf('\n(Earth mean equator and equinox of J2000)');
fprintf('\n-----------------------------------------\n');

display1(jdtdb_tip, emu, r_tip, v_tip);

[c3, rla, dla] = rv2hyper (emu, r_tip, v_tip);

fprintf('\nspecific orbital energy     %12.6f  (km/sec)**2\n', c3);

fprintf('\nasymptote right ascension   %12.6f  degrees\n', rtd * rla);

fprintf('\nasymptote declination       %12.6f  degrees\n\n', rtd * dla);

fprintf('\nEarth heliocentric coordinates at departure');
fprintf('\n(Earth mean equator and equinox of J2000)');
fprintf('\n-----------------------------------------\n');

svearth = jplephem(jdtdb_tip, 3, 11);

rearth = svearth(1:3);

vearth = svearth(4:6);

display1(jdtdb_tip, smu, aunit * rearth, aunit * vearth / 86400.0d0);

fprintf('\nEarth heliocentric coordinates at departure');
fprintf('\n(Earth mean ecliptic and equinox of J2000)');
fprintf('\n-----------------------------------------\n');

recl = eqwrk * rearth;

vecl = eqwrk * vearth;

display1(jdtdb_tip, smu, aunit * recl, aunit * vecl / 86400.0d0);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute time and conditions at the SOI
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

jdtdb_wrk = jdtdb_tip;

% set up for ode45

options = odeset('RelTol', 1.0e-8, 'AbsTol', 1.0e-10, 'Events', @soi_event);

% define maximum search time (seconds)

tof = 25.0 * 86400.0;

[t, ysol, tevent, yevent, ie] = ode45(@geo_eqm, [0 tof], [r_tip v_tip], options);

jdtdb_soi = jdtdb_tip + tevent / 86400.0;

% spacecraft geocentric state vector at soi (km & km/sec)

rgeo_soi = yevent(1:3);

vgeo_soi = yevent(4:6);

% heliocentric state vector of the spacecraft at soi (au & au/day)

svsun = jplephem(jdtdb_soi, 3, 11);

rsun = svsun(1:3);

vsun = svsun(4:6);

for i = 1:1:3
    
    rhelio_soi(i) = aunit * rsun(i) + yevent(i);
    
    vhelio_soi(i) = aunit * vsun(i) / 86400.0d0 + yevent(i + 3);
    
end

fprintf('\ntime and geocentric conditions at Earth SOI');
fprintf('\n(Earth mean equator and equinox of J2000)');
fprintf('\n-----------------------------------------\n');

display1(jdtdb_soi, emu, rgeo_soi, vgeo_soi);

fprintf('\ntime and heliocentric conditions at Earth SOI');
fprintf('\n(Earth mean equator and equinox of J2000)');
fprintf('\n-----------------------------------------\n');

display1(jdtdb_soi, smu, rhelio_soi, vhelio_soi);

if (iptype == 1)
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % find closest approach conditions
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    jdtdb_wrk = jdtdb_soi;
    
    ri = rhelio_soi / aunit;
    
    vi = 86400.0d0 * vhelio_soi / aunit;
    
    % compute closest approach conditions
    
    [jdtdb_ca, rm2sc, vm2sc] = findca(jdtdb_soi, ri, vi);
    
    fprintf('\ntime and conditions at Mars closest approach');
    fprintf('\n(areocentric mean equator and IAU node of epoch)');
    fprintf('\n------------------------------------------------\n');
    
    display1(jdtdb_ca, xmu_mars, rm2sc, vm2sc);
    
    % compute b-plane coordinates
    
    [bplane, ibperr] = rv2bp2(xmu_mars, rm2sc, vm2sc);
    
    % display results
    
    fprintf('\nb-plane coordinates of incoming hyperbola');
    fprintf('\n(areocentric mean equator and IAU node of epoch)');
    fprintf('\n------------------------------------------------\n');
    
    fprintf ('\nb-magnitude                %14.6f  kilometers', sqrt(bplane(7)^2 + bplane(8)^2));
    
    fprintf ('\nb dot t                    %14.6f  kilometers', bplane(7));
    
    fprintf ('\nb dot r                    %14.6f  kilometers', bplane(8));
    
    fprintf ('\nb-plane angle              %14.6f  degrees', rtd * bplane(6));
    
    fprintf ('\nv-infinity                 %14.6f  meters/second', 1000.0d0 * bplane(1));
    
    fprintf ('\nr-periapsis                %14.6f  kilometers', bplane(5));
    
    fprintf ('\ndecl-asymptote             %14.6f  degrees', rtd * bplane(2));
    
    fprintf ('\nrasc-asymptote             %14.6f  degrees', rtd * bplane(3));
    
    % flight path angle at entry interface
    
    sfpa = rm2sc' * vm2sc /(norm(rm2sc) * norm(vm2sc));
    
    fprintf('\n\nareocentric flight path angle     %12.10e  degrees\n\n', rtd * asin(sfpa));
    
    fprintf('\nMars heliocentric coordinates at closest approach');
    fprintf('\n(Earth mean equator and equinox of J2000)');
    fprintf('\n-----------------------------------------\n');
    
    svmars = jplephem(jdtdb_ca, 4, 11);
    
    rmars = svmars(1:3);
    
    vmars = svmars(4:6);
    
    display1(jdtdb_ca, smu, aunit * rmars, aunit * vmars / 86400.0d0);
    
    fprintf('\nMars heliocentric coordinates at closest approach');
    fprintf('\n(Earth mean ecliptic and equinox of J2000)');
    fprintf('\n------------------------------------------\n');
    
    recl = eqwrk * rmars;
    
    vecl = eqwrk * vmars;
    
    display1(jdtdb_ca, smu, aunit * recl, aunit * vecl / 86400.0d0);
    
    fprintf ('\nsequence of events summary');
    fprintf ('\n--------------------------');
    
    fprintf ('\n\ndeparture to SOI time         %12.6f  days', jdtdb_soi - jdtdb_tip);
    
    fprintf ('\n\nSOI to closest approach time  %12.6f  days', jdtdb_ca - jdtdb_soi);
    
    fprintf ('\n\nEarth to Mars transfer time   %12.6f  days\n', jdtdb_ca - jdtdb_tip);
    
end

if (iptype == 2)
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % compute conditions at user-defined areocentric distance
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    jdtdb_wrk = jdtdb_soi;
    
    ri = rhelio_soi / aunit;
    
    vi = 86400.0d0 * vhelio_soi / aunit;
    
    % compute areocentric distance conditions
    
    [jdtdb_rm, rm2sc, vm2sc] = findrm(jdtdb_soi, ri, vi);
    
    fprintf('\ntime and conditions at user-defined areocentric distance');
    fprintf('\n(areocentric mean equator and IAU node of epoch)');
    fprintf('\n------------------------------------------------\n');
    
    display1(jdtdb_rm, xmu_mars, rm2sc, vm2sc);
    
    fprintf('\nMars heliocentric coordinates at areocentric distance');
    fprintf('\n(Earth mean equator and equinox of J2000)');
    fprintf('\n-----------------------------------------\n');
    
    svmars = jplephem(jdtdb_rm, 4, 11);
    
    rmars = svmars(1:3);
    
    vmars = svmars(4:6);
    
    display1(jdtdb_rm, smu, aunit * rmars, aunit * vmars / 86400.0d0);
    
    fprintf('\nMars heliocentric coordinates at areocentric distance');
    fprintf('\n(Earth mean ecliptic and equinox of J2000)');
    fprintf('\n------------------------------------------\n');
    
    recl = eqwrk * rmars;
    
    vecl = eqwrk * vmars;
    
    display1(jdtdb_rm, smu, aunit * recl, aunit * vecl / 86400.0d0);
    
    fprintf ('\nsequence of events summary');
    fprintf ('\n--------------------------');
    
    fprintf ('\n\ndeparture to SOI time         %12.6f  days', jdtdb_soi - jdtdb_tip);
    
    fprintf ('\n\nSOI to areocentric distance   %12.6f  days', jdtdb_rm - jdtdb_soi);
    
    fprintf ('\n\nEarth to Mars transfer time   %12.6f  days\n', jdtdb_rm - jdtdb_tip);
    
end

if (iptype == 3)
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % compute conditions at user-defined final epoch
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    % set up for ode45
    
    options = odeset('RelTol', 1.0e-8, 'AbsTol', 1.0e-10);
    
    % define propagation duration (days)
    
    tof = jdtdb_user - jdtdb_soi;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%
    % set up initial conditions
    %%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    jdtdb_wrk = jdtdb_soi;
    
    ri = rhelio_soi / aunit;
    
    vi = 86400.0d0 * vhelio_soi / aunit;
 
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % solve for final conditions
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    [t, ysol] = ode45(@hel_eqm, [0 tof], [ri vi], options);
     
    % final heliocentric state vector (au and au/day)
    
    rf = ysol(end,(1:3));
     
    vf = ysol(end,(4:6));

    % compute mars-centered state vector of the spacecraft
    % (Mars mean equator and IAU node of epoch)

    svmars = jplephem(jdtdb_user, 4, 11);
    
    rmars = svmars(1:3);
    
    vmars = svmars(4:6);
    
    rm2sc = rf(1:3)' - rmars(1:3);

    vm2sc = vf(1:3)' - vmars(1:3);

    tmatrix = mme2000(jdtdb_user);

    rsc = tmatrix * rm2sc;

    vsc = tmatrix * vm2sc;
    
    rm2sc = aunit * rsc;
 
    vm2sc = aunit * vsc / 86400.0;

    fprintf('\ntime and conditions at user-defined final epoch');
    fprintf('\n(areocentric mean equator and IAU node of epoch)');
    fprintf('\n------------------------------------------------\n');
    
    display1(jdtdb_user, xmu_mars, rm2sc, vm2sc);
    
    fprintf('\nMars heliocentric coordinates at final epoch');
    fprintf('\n(Earth mean equator and equinox of J2000)');
    fprintf('\n-----------------------------------------\n');
        
    display1(jdtdb_user, smu, aunit * rmars, aunit * vmars / 86400.0d0);
    
    fprintf('\nMars heliocentric coordinates at final epoch');
    fprintf('\n(Earth mean ecliptic and equinox of J2000)');
    fprintf('\n------------------------------------------\n');
    
    recl = eqwrk * rmars;
    
    vecl = eqwrk * vmars;
    
    display1(jdtdb_user, smu, aunit * recl, aunit * vecl / 86400.0d0);
    
    fprintf ('\nsequence of events summary');
    fprintf ('\n--------------------------');
    
    fprintf ('\n\ndeparture to SOI time         %12.6f  days', jdtdb_soi - jdtdb_tip);
    
    fprintf ('\n\nSOI to final epoch            %12.6f  days', jdtdb_user - jdtdb_soi);
    
    fprintf ('\n\nEarth to Mars transfer time   %12.6f  days\n', jdtdb_user - jdtdb_tip);
    
end

fprintf('\n\n');

close('all');

Contact us