Code covered by the BSD License

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

18 Apr 2013

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');

```