Code covered by the BSD License

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

### David Eagle (view profile)

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

findrm(jdtdbi, ri, vi)
```function [jdtdb_rm, rm2sc, vm2sc] = findrm(jdtdbi, ri, vi)

% solve for areocentric distance conditions

% input

%  jdtdbi = initial tdb julian date
%  ri     = initial heliocentric position vector (au)
%  vi     = initial heliocentric velocity vector (au/day)

% output

%  jdtdb_rm = tdb julian date at user-defined areocentric distance
%  rm2sc    = mars-to-spacecraft position vector (kilometers)
%  vm2sc    = mars-to-spacecraft velocity vector (kilometers/second)

% NOTE: rm2sc and vm2sc => Mars mean equator and IAU node of epoch

% Orbital Mechanics with MATLAB

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

global aunit jdtdb_wrk

% set up for ode45

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

% define maximum search time (days)

tof = 500.0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve for areocentric distance conditions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

jdtdb_wrk = jdtdbi;

[t, ysol, tevent, yevent, ie] = ode45(@hel_eqm, [0 tof], [ri vi], options);

% tdb julian date at closest approach

jdate = jdtdbi + tevent;

jdtdb_rm = jdate;

% mars heliocentric state vector

svmars = jplephem (jdate, 4, 11);

rmars = svmars(1:3);

vmars = svmars(4:6);

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

rm2sc = yevent(1:3)' - rmars(1:3);

vm2sc = yevent(4:6)' - vmars(1:3);

tmatrix = mme2000(jdate);

rsc = tmatrix * rm2sc;

vsc = tmatrix * vm2sc;

% mars-to-spacecraft state vector (km and km/sec)
% (Mars mean equator and IAU node of epoch)

rm2sc = aunit * rsc;

vm2sc = aunit * vsc / 86400.0;
```