A MATLAB Script for Interplanetary Patched-Conic, Gravity-Assist

Trajectory Design and Optimization

This document describes a MATLAB script named flyby_matlab that can be used to design and
optimize interplanetary patched-conic trajectories that include a single gravity assist maneuver. The
user specifies the launch, flyby and destination planets along with initial guesses for the launch, flyby
and arrival calendar dates, and lower and upper bounds for the flyby altitude. This script searches for a
patched-conic gravity-assist trajectory that satisfies the flyby mission constraints (V matching and
flyby altitude) and minimizes the launch, arrival or total impulsive heliocentric delta-v for the mission.
The type of trajectory optimization is specified by the user.

The planet positions and velocities are modeled using the JPL DE 421 ephemeris. Information about
JPL ephemerides can be found at http://naif.jpl.nasa.gov/naif/. The trajectory optimization is performed
with the SNOPT nonlinear programming (NLP) algorithm. SNOPT mex files for several computer
platforms can be found at http://ccom.ucsd.edu/~optimizers/.

The flyby_matlab software provides a simple two-dimensional graphic display of the heliocentric
planet orbits and the interplanetary transfer trajectory. It also provides a three-dimensional graphic
display of the flyby trajectory.

Typical user interaction

The following is typical user interaction with this MATLAB application. This example is an Earth with
Venus gravity-assist to Mars mission. The user inputs for this example are in bold font.

The software allows the user to specify an initial guess for the launch, flyby and arrival calendar dates,
and lower and upper bounds on the actual launch, flyby and arrival calendar dates found during the
optimization process. For any guess for launch time t L and user-defined launch time lower and upper
bounds tl and tu , the launch time t is constrained as follows:

tL  tl  t  tL  tu

Likewise, for any guess for arrival time t A and user-defined arrival time bounds tl and tu , the arrival
time t is constrained as follows:
t A  tl  t  t A  tu

For fixed launch, flyby or arrival times, the lower and upper bounds should be set to 0.001.

The following is a typical user interaction with the flyby_matlab script. It defines one of the classical
JPL Earth-Venus-Mars missions analyzed during the 1960s.

program flyby_matlab

gravity-assist trajectory analysis

departure - calendar date guess

please input the calendar date

page 1
(1 <= month <= 12, 1 <= day <= 31, year = all digits!)
? 8,12,1970

please input the departure date search boundary in days

? 30

flyby - calendar date guess

please input the calendar date

(1 <= month <= 12, 1 <= day <= 31, year = all digits!)
? 12, 20, 1970

please input the flyby date search boundary in days

? 30

arrival - calendar date guess

please input the calendar date

(1 <= month <= 12, 1 <= day <= 31, year = all digits!)
? 6, 13, 1971

please input the arrival date search boundary in days

? 30

planet menu

<1> Mercury
<2> Venus
<3> Earth
<4> Mars
<5> Jupiter
<6> Saturn
<7> Uranus
<8> Neptune
<9> Pluto

please select the departure planet

? 3

planet menu

<1> Mercury
<2> Venus
<3> Earth
<4> Mars
<5> Jupiter
<6> Saturn
<7> Uranus
<8> Neptune
<9> Pluto

please select the flyby planet

? 2

planet menu

<1> Mercury
<2> Venus
<3> Earth
page 2
<4> Mars
<5> Jupiter
<6> Saturn
<7> Uranus
<8> Neptune
<9> Pluto

please select the arrival planet

? 4

please input the lower bound for the flyby altitude (kilometers)
? 500

please input the upper bound for the flyby altitude (kilometers)
? 5000

optimization menu

<1> minimize departure delta-v

<2> minimize arrival delta-v

<3> minimize total delta-v

selection (1, 2 or 3)
? 1

Optimal solution and trajectory graphics

This section summarizes the program output for this example. The software provides the heliocentric
orbital elements and state vectors of each leg of the transfer trajectory in the Earth mean ecliptic and
equinox of J2000 coordinate system. These results also include the characteristics of the hyperbolic
flyby trajectory with respect to the flyby planet as well as the trajectory characteristics at the time of
entrance to the sphere-of-influence (SOI) and closest approach relative to the flyby planet. The time
scale is Barycentric Dynamical Time (TDB). See Appendix A for a description of this information.

program flyby_matlab
--------------------

patched-conic gravity assist trajectory analysis

minimize departure delta-v

departure planet 'Earth'

flyby planet 'Venus'

arrival planet 'Mars'

LAUNCH CONDITIONS
=================

departure calendar date 12-Aug-1970

departure TDB time 10:26:30.851

departure TDB julian date 2440810.935079

heliocentric launch delta-v vector and magnitude

(Earth mean ecliptic and equinox of J2000)
------------------------------------------

page 3
launch delta-vx -1520.218236 meters/second
launch delta-vy -2163.716837 meters/second
launch delta-vz 1903.009119 meters/second

launch delta-v 3257.940722 meters/second

launch hyperbola characteristics

(Earth mean equator and equinox of J2000)
-----------------------------------------

asymptote right ascension 240.996444 degrees

asymptote declination 15.767592 degrees

launch energy 10.614178 (km/sec)^2

post-impulse heliocentric orbital elements and state vector of the spacecraft

(Earth mean ecliptic and equinox of J2000)
------------------------------------------

sma (km) eccentricity inclination (deg) argper (deg)

+1.28621576607637e+08 +1.78516079861672e-01 +4.07175310704317e+00 +1.79283757715307e+02

raan (deg) true anomaly (deg) arglat (deg) period (days)

+3.19739740270847e+02 +1.80687567096999e+02 +3.59971324812306e+02 +2.91193081269169e+02

rx (km) ry (km) rz (km) rmag (km)

+1.15624492017452e+08 -9.80180660555501e+07 -5.38665739078075e+03 +1.51580224490466e+08

vx (kps) vy (kps) vz (kps) vmag (kps)

+1.72450080678064e+01 +2.04506760424211e+01 +1.90428162500954e+00 +2.68187759295258e+01

heliocentric orbital elements and state vector of the departure planet

(Earth mean ecliptic and equinox of J2000)
------------------------------------------

sma (km) eccentricity inclination (deg) argper (deg)

-4.67266870810208e+02 +3.24380560886372e+05 +3.19314885605842e-03 +3.20986692182320e+02

raan (deg) true anomaly (deg) arglat (deg)

+3.59327717669158e+02 +3.59396727566710e+02 +3.20383419749030e+02

rx (km) ry (km) rz (km) rmag (km)

+1.15624492017452e+08 -9.80180660555501e+07 -5.38665739078075e+03 +1.51580224490466e+08

vx (kps) vy (kps) vz (kps) vmag (kps)

+1.87652263038200e+01 +2.26143928798734e+01 +1.27250579574323e-03 +2.93861274273783e+01

PATCHED-CONIC FLYBY CONDITIONS

==============================

flyby calendar date 19-Dec-1970

flyby TDB time 17:27:19.126

flyby TDB julian date 2440940.227305

launch-to-flyby time 129.292225 days

v-infinity in 5471.917891 meters/second

v-infinity out 5473.022379 meters/second

flyby altitude 3523.012352 kilometers

maximum turn angle 79.872068 degrees

actual turn angle 64.173912 degrees

heliocentric delta-v 5814.014808 meters/second

max heliocentric delta-v 7326.580266 meters/second

page 4
planet-centered orbital elements and state vector of the spacecraft at periapsis
(Earth mean ecliptic and equinox of J2000)
------------------------------------------

sma (km) eccentricity inclination (deg) argper (deg)

-1.08496373748032e+04 +1.88250989604897e+00 +1.12470605615357e+02 +2.46605055719858e+02

raan (deg) true anomaly (deg) arglat (deg)

+2.16827139290049e+02 +0.00000000000000e+00 +2.46605055719858e+02

rx (km) ry (km) rz (km) rmag (km)

+5.05645806848405e+03 -4.09655391642531e+02 -8.12055175505567e+03 +9.57491235180651e+03

vx (kps) vy (kps) vz (kps) vmag (kps)

-6.36549461062597e+00 -6.64168576216376e+00 -3.62857784674308e+00 +9.88929161219911e+00

b-plane coordinates of the spacecraft at periapsis

(Earth mean ecliptic and equinox of J2000)
------------------------------------------

b-magnitude 17304.554325 kilometers

b dot r 13977.297028 kilometers
b dot t -10202.096263 kilometers
b-plane angle 126.125886 degrees
v-infinity 5471.917891 meters/second
r-periapsis 9574.912352 kilometers
decl-asymptote -49.586859 degrees
rasc-asymptote 245.889904 degrees

flight path angle -0.000000 degrees

heliocentric orbital elements and state vector of the spacecraft

(Earth mean ecliptic and equinox of J2000)
------------------------------------------

sma (km) eccentricity inclination (deg) argper (deg)

+1.28621576607637e+08 +1.78516079861672e-01 +4.07175310704317e+00 +1.79283757715307e+02

raan (deg) true anomaly (deg) arglat (deg) period (days)

+3.19739740270847e+02 +3.32351926052002e+02 +1.51635683767308e+02 +2.91193081269169e+02

rx (km) ry (km) rz (km) rmag (km)

-3.92720261161839e+07 +1.00025789550060e+08 +3.62699388823714e+06 +1.07520257138124e+08

vx (kps) vy (kps) vz (kps) vmag (kps)

-3.41657184360208e+01 -1.62321242286620e+01 -2.45354380902881e+00 +3.79051190533283e+01

heliocentric orbital elements and state vector of the flyby planet

(Earth mean ecliptic and equinox of J2000)
------------------------------------------

sma (km) eccentricity inclination (deg) argper (deg)

+1.08209176115514e+08 +6.76716812755798e-03 +3.39503633361855e+00 +5.46696519770352e+01

raan (deg) true anomaly (deg) arglat (deg) period (days)

+7.67589526806317e+01 +3.40054415742044e+02 +3.47240677190790e+01 +2.24701744670303e+02

rx (km) ry (km) rz (km) rmag (km)

-3.92720261161870e+07 +1.00025789550058e+08 +3.62699388823691e+06 +1.07520257138124e+08

vx (kps) vy (kps) vz (kps) vmag (kps)

-3.27166304255047e+01 -1.29941782834181e+01 +1.71271775284402e+00 +3.52442900022291e+01

ARRIVAL CONDITIONS
==================

arrival calendar date 18-Jun-1971

arrival TDB time 15:02:15.491

page 5
arrival TDB julian date 2441121.126568

flyby-to-arrival time 180.899263 days

heliocentric arrival delta-v vector and magnitude

(Earth mean ecliptic and equinox of J2000)
------------------------------------------

arrival delta-vx 4025.678393 meters/second

arrival delta-vy 5266.112483 meters/second
arrival delta-vz 974.578869 meters/second

arrival delta-v 6699.838146 meters/second

pre-impulse heliocentric orbital elements and state vector of the spacecraft

(Earth mean ecliptic and equinox of J2000)
------------------------------------------

sma (km) eccentricity inclination (deg) argper (deg)

+1.61763582089572e+08 +3.36451361464569e-01 +4.04871782743126e+00 +3.51927214079591e+01

raan (deg) true anomaly (deg) arglat (deg) period (days)

+8.29558001291680e+01 +1.66944008491318e+02 +2.02136729899277e+02 +4.10706372959098e+02

rx (km) ry (km) rz (km) rmag (km)

+5.53636002907808e+07 -2.06006926232007e+08 -5.67732366320330e+06 +2.13392159876350e+08

vx (kps) vy (kps) vz (kps) vmag (kps)

+2.02946344314881e+01 +3.09799699568943e+00 -1.39874653515199e+00 +2.05773240233672e+01

post-impulse heliocentric orbital elements and state vector of the spacecraft

(Earth mean ecliptic and equinox of J2000)
------------------------------------------

sma (km) eccentricity inclination (deg) argper (deg)

+2.27941975487972e+08 +9.32907670879080e-02 +1.85199968562863e+00 +2.86269612640230e+02

raan (deg) true anomaly (deg) arglat (deg) period (days)

+4.96468871616997e+01 +3.09140095836467e+02 +2.35409708476697e+02 +6.86984234438126e+02

rx (km) ry (km) rz (km) rmag (km)

+5.53636002907808e+07 -2.06006926232007e+08 -5.67732366320330e+06 +2.13392159876350e+08

vx (kps) vy (kps) vz (kps) vmag (kps)

+2.43203128247463e+01 +8.36410947856076e+00 -4.24167665977691e-01 +2.57218945933563e+01

heliocentric orbital elements and state vector of the arrival planet

(Earth mean ecliptic and equinox of J2000)
------------------------------------------

sma (km) eccentricity inclination (deg) argper (deg)

+2.27941975487972e+08 +9.32907670879081e-02 +1.85199968562863e+00 +2.86269612640230e+02

raan (deg) true anomaly (deg) arglat (deg) period (days)

+4.96468871616997e+01 +3.09140095836467e+02 +2.35409708476697e+02 +6.86984234438126e+02

rx (km) ry (km) rz (km) rmag (km)

+5.53636002907808e+07 -2.06006926232007e+08 -5.67732366320330e+06 +2.13392159876350e+08

vx (kps) vy (kps) vz (kps) vmag (kps)

+2.43203128247463e+01 +8.36410947856076e+00 -4.24167665977691e-01 +2.57218945933563e+01

MISSION SUMMARY
===============

total delta-v 9957.778867 meters/second

total energy 55.502009 (km/sec)^2

page 6
total mission duration 310.191489 days

PATCHED-CONIC SOI ENTRANCE CONDITIONS

=====================================

calendar date 18-Dec-1970

TDB time 10:09:33.081

TDB julian date 2440938.923300

planet-centered orbital elements and state vector of the spacecraft

(Earth mean ecliptic and equinox of J2000)
------------------------------------------

sma (km) eccentricity inclination (deg) argper (deg)

-1.12497057136791e+04 +1.00092122453233e+00 +1.09325179245321e+02 +5.12951225154509e+01

raan (deg) true anomaly (deg) arglat (deg)

+9.02461540050350e+01 +1.82502908292374e+02 +2.33798030807825e+02

rx (km) ry (km) rz (km) rmag (km)

-1.63005311375096e+05 -3.64697523781747e+05 -4.69278685541624e+05 +6.16277129295502e+05

vx (kps) vy (kps) vz (kps) vmag (kps)

+1.44622920763433e+00 +3.24096878808079e+00 +4.16363772539961e+00 +5.47095391268281e+00

b-plane coordinates at sphere-of-influence

(Earth mean ecliptic and equinox of J2000)
------------------------------------------

b-magnitude 482.990851 kilometers

b dot r 415.416715 kilometers
b dot t -246.392198 kilometers
b-plane angle 120.673051 degrees
v-infinity 5373.739415 meters/second
r-periapsis 10.363505 kilometers
decl-asymptote 49.555999 degrees
rasc-asymptote 65.952276 degrees

flight path angle -89.955894 degrees

NUMERICALLY INTEGRATED CLOSEST APPROACH CONDITIONS

==================================================

calendar date 19-Dec-1970

TDB time 15:50:49.735

TDB julian date 2440940.160298

planet-centered orbital elements and state vector of the spacecraft

(Earth mean ecliptic and equinox of J2000)
------------------------------------------

sma (km) eccentricity inclination (deg) argper (deg)

-1.12452080579129e+04 +1.00015973695730e+00 +8.24925749879297e+01 +4.91468615457492e+01

raan (deg) true anomaly (deg) arglat (deg)

+5.69886738092402e+01 +1.77501624375001e-05 +4.91468792959116e+01

rx (km) ry (km) rz (km) rmag (km)

+4.91276439771532e-01 +1.08201034283853e+00 +1.34703604287792e+00 +1.79627531940300e+00

vx (kps) vy (kps) vz (kps) vmag (kps)

-2.90947905631954e+02 -3.53477931154337e+02 +3.90043606367791e+02 +6.01441058187053e+02

page 7
b-plane coordinates at closest approach
(Earth mean ecliptic and equinox of J2000)
------------------------------------------

b-magnitude 201.003000 kilometers

b dot r 196.878860 kilometers
b dot t 40.508276 kilometers
b-plane angle 78.373515 degrees
v-infinity 5374.813953 meters/second
r-periapsis 1.796275 kilometers
decl-asymptote 49.585583 degrees
rasc-asymptote 65.891966 degrees

flight path angle 0.000009 degrees

After the solution is displayed, the software will ask the user if he or she would like to create a graphics
display of the transfer and flyby trajectories with the following prompt:
would you like to create graphics for this mission (y = yes, n = no)
?

If the user’s response is y for yes, the script will request a plot step size with
please input the plot step size (days)
?

A plot step size between 5 and 10 days is recommended.

The following is the heliocentric graphics display for this example. This plot is a north ecliptic view
where we are looking down on the ecliptic plane from the north celestial pole. The vernal equinox
direction is the labeled line pointing to the right, the launch planet is labeled with an L, the flyby planet
is labeled with an F and the arrival planet is labeled with an A. The locations of the launch, flyby and
arrival planets at the launch time are marked with an asterisk. Please note that the x and y axes are
labeled in Astronomical Units (AU).

page 8
The following is the flyby graphics display for this example. The direction of the sun is the yellow line
and the black line points in the direction of the heliocentric velocity of the flyby planet. This display is
labeled with a planet-centered, inertial mean ecliptic and equinox of J2000 coordinate system. The x-
axis of this system is red, the y-axis is green and the z-axis is blue. The flyby hyperbola is colored
magenta with the incoming leg marked with an asterisk and periapsis labeled with a small circle. Please
note that the axes are labeled in radii of the flyby planet (PR).

This script will also create disk files of these graphic displays. These files are color Postscript images
with a TIFF preview. The name and other characteristics of these files can be edited in the following
lines of MATLAB source code;
print -depsc -tiff -r300 flyby_matlab1.eps

print -depsc -tiff -r300 flyby_matlab2.eps

Technical discussion

The computational steps used to create an initial guess for the spacecraft state, and the launch and arrival
delta-v characteristics are as follows:

(1) compute the state vector of the departure planet at the departure date initial guess

(2) compute the state vector of the flyby planet at the flyby date initial guess

(3) compute the state vector of the arrival planet at the arrival date initial guess

(4) solve Lambert's problem for the departure-to-flyby leg and determine the initial velocity vector
of this leg
page 9
(5) compute the launch delta-v vector from the launch planet’s velocity vector and the initial
velocity vector of the first leg of the transfer trajectory

(6) solve Lambert's problem for the second heliocentric leg and determine the initial and final
velocity vectors of the second leg

(7) compute the arrival delta-v vector from the arrival planet’s velocity vector and the final velocity
vector of the second leg of the transfer trajectory

(8) estimate the flyby incoming v-infinity vector from the flyby planet’s velocity vector and the final
velocity vector of the first leg of the transfer trajectory

Lambert’s Problem

Lambert’s problem is concerned with the determination of an orbit that passes between two positions
within a specified time-of-flight. This classic astrodynamics problem is also known as the orbital two-
point boundary value problem (TPBVP).

The time to traverse a trajectory depends only upon the length of the semimajor axis a of the transfer
trajectory, the sum ri  rf of the distances of the initial and final positions relative to a central body, and
the length c of the chord joining these two positions. This relationship can be stated as follows:

tof  tof  ri  rf , c, a 

From the following form of Kepler’s equation

a3
t  t0   E  e sin E 


we can write
a3
t  E  E0  e  sin E  sin E0 
 

where E is the eccentric anomaly associated with radius r, E0 is the eccentric anomaly at r0 , and t  0
when r  r0 .

At this point we need to introduce the following trigonometric sun and difference identities:

 a
sin   sin   2sin cos
2 2
 a
cos   cos   2sin sin
2 2
 a
cos   cos   2 cos cos
2 2

page 10
If we let E   and E0   and substitute the first trig identity into the second equation above, we have
the following equation:

a3  E  E0  E  E0  
t  E  E0  2sin  e cos 
  2  2 

With the two substitutions given by

E  E0 
e cos  cos
2 2
E  E0 
sin  sin
2 2
the time equation becomes

a3   
t      2sin cos 
 2 2 

From the elliptic relationships given by

r  a 1  e cos E 

x  a  cos E  e 

y  a sin E 1  e 2

and some more manipulation, we have the following equations:

 r  r0  c r  r0  c s
cos    1    1  1
 2a  2a 2a a

 r  r0  c r  r0  c sc
sin    1    1  1
 2a  2a 2a a

This part of the derivation makes use of the following three relationships:

  r  r0
cos cos  1
2 2 2

  E  E0 E  E0 
2

sin sin  sin 1   e cos 
2 2 2  2 

       x  x0   y  y0   c 
2 2 2 2

 sin sin       
 2 2   2a   2a   2a 

With the use of the half angle formulas given by

page 11
 s  sc
sin  sin 
2 2a 2 2a

and several additional substitutions, we have the time-of-flight form of Lambert’s theorem

a3
t       sin   sin  
 

A discussion about the angles  and  can be found in “Geometrical Interpretation of the Angles  and
 in Lambert’s Problem” by J. E. Prussing, AIAA Journal of Guidance and Control, Volume 2,
Number 5, Sept.-Oct. 1979, pages 442-443.

The Lambert algorithm used in this MATLAB script is based on the method described in “A Procedure
for the Solution of Lambert’s Orbital Boundary-Value Problem” by R. H. Gooding, Celestial Mechanics
and Dynamical Astronomy 48: 145-165, 1990. This iterative solution is valid for elliptic, parabolic and
hyperbolic transfer orbits which may be either posigrade or retrograde, and involve one or more
revolutions about the central body.

Gravity-assist flight mechanics

The vector relationship between the incoming v-infinity vector v  , the outgoing v-infinity vector v  and
the two legs of the heliocentric transfer orbit are as follows:

v   v f b  v to1

v   v to2  v f b
where

v f b  heliocentric velocity vector of the flyby planet at the flyby date

vto1  heliocentric velocity vector of the first transfer orbit at the flyby date
vto2  heliocentric velocity vector of the second transfer orbit at the flyby date

The turn angle of the planet-centered trajectory during the flyby is determined from

1 1  v   v  
 
  1 
   sin  2sin 1 
2 2  v v 
   1  rp v /  
2
   

where rp is the periapsis radius of the flyby hyperbola, v is the magnitude of the incoming (or
outgoing) v-infinity vector and  is the gravitational constant of the flyby planet.

The maximum turn angle possible during a gravity assist flyby occurs when the spacecraft just gazes the
planet’s surface. It is given by

page 12
 1 
 max  2sin 1  
 1  rp v2 /  
 

The semimajor axis and orbital eccentricity of the flyby hyperbola are given by

a    / v     / v 
2 2

rp rp v2
e  1/ cos   1   1
a 

where   is the true anomaly at infinity which is determined from the following expression:

1 1  v   v  
 

   sin
2 2  v  v  
 

The periapsis radius of the flyby hyperbola is determined from the expression r  a 1  e  and the flyby
altitude is h  r  rs where rs is the radius of the flyby planet..

The heliocentric speed gained during the flyby and the heliocentric delta-v vector caused by the close
encounter can be determined from the following two equations:

v fb  2v / e

v fb  v h  v h

where e is the orbital eccentricity of the hyperbolic flyby trajectory.

In the second equation v h is the heliocentric velocity vector of the spacecraft prior to the flyby and v h is
the heliocentric velocity vector after the flyby. For any planet it can be shown that the maximum
heliocentric delta-v possible from a closes flyby is given by the expression


vmax 
rs

This corresponds to a “grazing” flyby at the planet’s surface and is equal to the “local circular velocity”
for the flyby planet at the surface.

During the optimization analysis, the software enforces the following two nonlinear equality point
constraints:
v  v  0
hl  h f b  hu

page 13
The first equation is the v-infinity matching equality constraint and the second equation is the flyby
altitude inequality constraint. In the second expression h f b is the actual flyby altitude computed by the
software, hl is the user-defined lower bound, and hu is the user-defined upper bound for the flyby
altitude. To enforce a “fixed” flyby altitude, set hl  hu .

Departure trajectory

The orientation of the departure hyperbola is specified in terms of the right ascension and declination of
the asymptote of the launch hyperbola. These coordinates can be calculated using the components of the
planet-centered departure unit v-infinity vector.

The inertial right ascension of the asymptote is determined from


  tan 1 vˆ , vˆ
y x

and the geocentric declination of the asymptote is given by

  900  cos1  vˆ z


where vˆx , vˆ y , vˆz are the Cartesian components of the unit v-infinity vector.

In a typical “targeting spec” for an interplanetary mission, the right ascension is called RLA and the
declination is called DLA. The corresponding launch energy is called C3 which is equal to twice the
specific (per unit mass) orbital energy of the trajectory. The actual geocentric hyperbolic injection
delta-v depends on the size, orientation and other characteristics of the Earth park orbit.

In this MATLAB script the heliocentric planetary coordinates and therefore the v vectors are
computed in the J2000 mean ecliptic and equinox coordinate system. In order to determine the
orientation of the geocentric departure hyperbola, the heliocentric delta-v vector must be transformed to
the Earth mean equator and equinox of J2000 (EME2000) frame.

The required transformation is given by

1 0 0  1 0 0 

v eq  0 cos   sin  v ec  0 0.917482062069182 -0.397777155931914 v ec
 
   
0 sin  cos   0 0.397777155931914 0.917482062069182 

where v ec is the delta-velocity vector in the ecliptic frame, v eq is the delta-velocity vector in the
Earth equatorial frame and  is the mean obliquity of the ecliptic at J2000.

The J2000 value of the mean obliquity of the ecliptic is equal to   23 2621.448 .

The conversion of coordinates from the equatorial system to the ecliptic system involves the transpose
of this fundamental matrix.

page 14
The B-plane

The derivation of B-plane coordinates is described in the classic JPL reports, “A Method of Describing
Miss Distances for Lunar and Interplanetary Trajectories” and “Some Orbital Elements Useful in Space
Trajectory Calculations”, both by William Kizner. The following diagram illustrates the fundamental
geometry of the B-plane coordinate system.

The arrival asymptote unit vector Ŝ is given by

cos   cos   
 
Sˆ   cos   sin   
 sin  
  

where   and   are the declination and right ascension of the asymptote of the incoming or outgoing
hyperbola.

The following computational steps summarize the calculation of the predicted B-plane vector from a
planet-centered position vector r and velocity vector v at the sphere-of-influence or closest approach.

angular momentum vector

h rv
radius rate
r r v r

page 15
semiparameter
p  h2 

semimajor axis
r
a
 r v2 
 
 
2


orbital eccentricity
e  1 p a

true anomaly
pr rh
cos   sin  
er e

B-plane magnitude
B pa

fundamental vectors
r v  rr
zˆ 
h

pˆ  cos rˆ  sin  zˆ qˆ  sin  rˆ  cos zˆ

S vector
a b
S pˆ  qˆ
a 2  b2 a 2  b2

B vector
b2 ab
B pˆ  qˆ
a 2  b2 a 2  b2

T vector
 S ,  S ,0
2 2 T

T
y x

S x2  S y2

R vector
R  S  T    SzTy , SzTx , S xTy  S yTx 
T

page 16
Calculating Sphere-of-Influence Entrance and Closest Approach

This section describes the numerical methods used to predict the time at which the spacecraft enters the
sphere-of-influence (SOI) of the flyby planet and the time of closest approach to the flyby planet. These
calculations can be used to assess the difference between the patched-conic and numerically integrated
trajectory for the flyby hyperbola.

The sphere-of-influence objective function is given by

  rp sc  rSOI

where rp  sc is the radius vector of the spacecraft relative to the flyby planet and rSOI is the SOI radius of
the flyby planet. The SOI radius is a function of the gravity of the flyby planet.

During the search for SOI conditions, the heliocentric equations of motion of the spacecraft subject to
the point-mass gravity of the sun are defined by

rs  sc
a   s 3
rs  sc

where  s is the gravitational constant of the sun and rs  sc is the heliocentric position vector from the sun
to the spacecraft.

The following is the MATLAB source code from the main script that predicts SOI entrance. It uses the
root-finder or events option of the built-in ODE45 function to evaluate the two-body heliocentric
equations of motion defined in a MATLAB function named twobody_eqm while searching for the time
at which the flyby planet-to-spacecraft distance equals the SOI radius of the flyby planet.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% two-body integration of the trajectory from
% first impulse to SOI entrance at flyby planet
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% set up for ode45

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

% define maximum search time (seconds)

tof = 86400.0 * (jdtdb2 - jdtdb1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve for flyby planet SOI entrance conditions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[~, ~, tevent, ~, ie] = ode45(@twobody_heqm, [0 tof], [rito1 vito1], options);

The following is the MATLAB source code that computes the scalar value of the SOI objective
function.

page 17
function [value, isterminal, direction] = soi_event(t, y)

% entrance to SOI of flyby planet event function

% input

% t = time since departure (seconds)

% y = spacecraft heliocentric state vector at time t (km & km/sec)

% output

% value = delta-SOI (kilometers)

% Orbital Mechanics with MATLAB

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

global ip2 rsoi jdtdb1 rp2sc vp2sc

% current TDB julian date

jdtdb = jdtdb1 + t / 86400.0d0;

% compute the position and velocity of the flyby planet (km & km/sec)

[rfbp, vfbp] = p2000(ip2, jdtdb);

rsc = y(1:3);

vsc = y(4:6);

% state vector from flyby planet to spacecraft (km & km/sec)

rp2sc = rfbp - rsc;

vp2sc = vfbp - vsc;

% objective function (delta-SOI distance; kilometers)

value = norm(rp2sc) - rsoi(ip2);

isterminal = 1;

direction = [];

The software uses a similar technique to predict the spacecraft flight conditions at closest approach to
the flyby planet. The flyby planet-centered initial conditions for the search are determined from the
state vector of the spacecraft at SOI entrance.

For close approach prediction, the script uses the following objective function

rp  sc v p  sc
  sin 
rp  sc v p  sc

which is simply the sine of the flight path angle  . Since the spacecraft’s orbit within the SOI is
hyperbolic, periapsis occurs whenever the flight path angle is zero.
page 18
During the search for closest approach conditions, the flyby planet-centered equations of motion subject
to the point-mass gravity of the sun and the flyby planet are computed using

 r r 
rp  sc  
a  a p  as    p 3
 s  s sc 3  p 3 
rp  sc  r
 s sc rp  

In this equation, a is the total acceleration on the spacecraft and rs  sc is the heliocentric position vector
of the spacecraft given by
rssc  rpsc  rp

where rp is the heliocentric position vector of the flyby planet.

These calculations are performed by a MATLAB function named peqm. Here is the source code.

function ydot = peqm (t, y)

% first-order flyby-planet-centered equations of motion

% NOTE: includes point-mass gravity of the sun

% input

% t = time since SOI (seconds)

% y = planet-centered state vector at time t (km & km/sec)

% output

% ydot = integration vector (km/sec^2)

% Orbital Mechanics with MATLAB

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

global ip2 smu pmu jdtdb_soi

% current TDB julian date

jdtdb = jdtdb_soi + t / 86400.0d0;

% get current flyby planet heliocentric state vector (km & km/sec)

[rs2p, vs2p] = p2000(ip2, jdtdb);

% calculate position vector of sun relative to flyby planet (kilometers)

rp2s(1) = -rs2p(1);
rp2s(2) = -rs2p(2);
rp2s(3) = -rs2p(3);

rp2sm = norm(rp2s);

% compute heliocentric position vector of spacecraft (kilometers)

for i = 1:1:3
page 19
rs2sc(i) = y(i) + rs2p(i);

end

rs2scm = norm(rs2sc);

rmag = norm(y(1:3));

r3 = rmag * rmag * rmag;

% acceleration due to the flyby planet (km/sec^2)

accp(1) = -pmu(ip2) * y(1) / r3;

accp(2) = -pmu(ip2) * y(2) / r3;

accp(3) = -pmu(ip2) * y(3) / r3;

% acceleration due to the sun (km/sec^2)

accs(1) = -smu * (rs2sc(1) / rs2scm^3 + rp2s(1) / rp2sm^3);

accs(2) = -smu * (rs2sc(2) / rs2scm^3 + rp2s(2) / rp2sm^3);

accs(3) = -smu * (rs2sc(3) / rs2scm^3 + rp2s(3) / rp2sm^3);

% compute total flyby planet-centered acceleration vector (km/sec^2)

ydot(1) = y(4);

ydot(2) = y(5);

ydot(3) = y(6);

ydot(4) = accs(1) + accp(1);

ydot(5) = accs(2) + accp(2);

ydot(6) = accs(3) + accp(3);

ydot = ydot';

Here’s the MATLAB source code that performs the closest approach calculations.

% set up for ode45

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

% define maximum search time (seconds)

tof = 10.0 * 86400.0;

[~, ~, tevent, yevent, ie] = ode45(@peqm, [0 tof], [rp2sc vp2sc], options);

In this code, rp2sc and vp2sc represent the flyby planet-centered position and velocity vectors at
entrance to the SOI.

page 20
The MATLAB source code that evaluates the flight path event function is as follows.

function [value, isterminal, direction] = fpa_event(~, y)

% flight path angle wrt flyby planet event function

% input

% y = spacecraft planet-centered state vector (km & km/sec)

% output

% value = sine of flight path angle wrt flyby planet

% NOTE: units for state vector must be consistent

% Orbital Mechanics with MATLAB

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

% extract the spacecraft position and velocity

% vectors relative to the flyby planet (km & km/sec)

rsc = y(1:3);

vsc = y(4:6);

% sine of the flight path angle

value = dot(rsc, vsc) / (norm(rsc) * norm(vsc));

isterminal = 1;

direction = [];

SNOPT algorithm implementation

This section provides details about the part of the flyby_matlab MATLAB script that solves this
nonlinear programming (NLP) problem using the SNOPT algorithm. In this classic trajectory
optimization problem, the launch, flyby and arrival calendar dates are the control variables and the user-
specified scalar V is the objective function or performance index.

The SNOPT algorithm requires an initial guess for the control variables. For this problem they are given
by the following code

% control variables initial guess

xg(1) = jdtdbi1 - jdtdb0;

xg(2) = jdtdbi2 - jdtdb0;

xg(3) = jdtdbi3 - jdtdb0;

xg = xg';

page 21
where jdtdbi1, jdtdbi2 and jdtdbi3 are the initial user-provided launch, flyby and arrival calendar
date guesses, and jdtdb0 is a reference Julian date equal to 2451544.5 (January 1, 2000 TDB). This
offset value is used to scale the control variables.

The algorithm also requires lower and upper bounds for the control variables. These are determined
from the user-defined initial guesses and search boundaries as follows:

% bounds on control variables

xlwr(1) = xg(1) - ddays1;

xupr(1) = xg(1) + ddays1;

xlwr(2) = xg(2) - ddays2;

xupr(2) = xg(2) + ddays2;

xlwr(3) = xg(3) – ddays3;

xupr(3) = xg(3) + ddays3;

xlwr = xlwr';
xupr = xupr';

xlwr = xlwr';
xupr = xupr';

where ddays1, ddays2 and ddays3 are the user-defined launch, flyby and arrival search boundaries,
respectively.

The algorithm requires lower and upper bounds on the objective function. For this problem these
bounds are given by
% bounds on objective function

flow(1) = 0.0d0;
fupp(1) = +Inf;

The following MATLAB source code enforces the flyby constraints:

% bounds on flyby altitude

flow(2) = fbalt_lwr;
fupp(2) = fbalt_upr;

% bounds on v-infinity matching (equality) constraint

flow(3) = 0.0;
fupp(3) = 0.0;

In this code fbalt, vinfm_in, and vinfm_out are the current computed flyby altitude, incoming and
outgoing v-infinity speeds, respectively. fbalt_lwr and fbalt_upr are the user-defined lower and
upper bounds on the flyby altitude.

The actual call to the SNOPT MATLAB interface function is as follows

[x, f, inform, xmul, fmul] = snopt(xg, xlwr, xupr, xmul, xstate, ...
flow, fupp, fmul, fstate, 'fbyfunc');

page 22
where fbyfunc is the name of the MATLAB function that solves Lambert’s problem for both legs of
the interplanetary transfer trajectory and computes the current value of the objective function.

Algorithm resources

“Gravity-Assist Trajectory Design with MATLAB”, C. David Eagle, AAS 99-138, AAS/AIAA Space
Flight Mechanics Meeting, Breckenridge, CO, February 7-10, 1999.

“Optimal Interplanetary Orbit Transfers by Direct Transcription”, John T. Betts, The Journal of the
Astronautical Sciences, Vol. 42, No. 3, July-September 1994, pp. 247-268.

“Multiple Gravity Assist Interplanetary Trajectories”, A. V. Labunsky, O. V. Papkov, and K. G.

Sukhanov, Gordon and Breach Science Publishers, 1998.

“Modern Astrodynamics”, Victor R. Bond and Mark C. Allman, Princeton University Press, 1996.

“Automated Design of Gravity-Assist Trajectories to Mars and the Outer Planets”, J.M. Longuski and
S.N. Williams, Celestial Mechanics and Dynamical Astronomy, 52: 207-220, 1991.

“Gravitational Assist in Celestial Mechanics – A Tutorial”, J. A. Van Allen, American Journal of

Physics, 71 (5), May 2003.

“Notes for the Gravitational Assisted Trajectories Lectures”, E. Barrabes, G. Gomez, and J. Rodriquez-
Canabal, Advanced Topics in Astrodynamics Summer Course, Barcelona, July 2004.

“Error Analysis of Multiple Planet Trajectories”, F. M. Sturms Jr., JPL Space Programs Summary, No.
37-27, Vol. IV.

“Trajectory Analysis of a 1970 Mission to Mercury via a Close Encounter with Venus”, E. Cutting, F.
M. Sturms Jr., AIAA Journal of Spacecraft and Rockets, Vol.3, No. 5, pp. 624-631, 1966.

“A Graphical Method for Gravity-assist Trajectory Design”, Nathan Strange, James Longuski,
Astrodynamics Specialist Conference, 2000.

“Mars Free Returns via Gravity Assist from Venus”, Masataka Okutsu, James Longuski, Astrodynamics
Specialist Conference, 2000.

“Application of Tisserand’s Criterion to the Design of Gravity Assist Trajectories”, James Miller,
Connie Weeks, AIAA/AAS Astrodynamics Specialist Conference and Exhibit, 2002.

“Fast Mars Free-Returns via Venus Gravity Assist”, Kyle M. Hughes, Peter J. Edelman, James
Longuski, Mike Loucks, John P. Carrico, Dennis A. Tito, AIAA/AAS Astrodynamics Specialist
Conference, 2014.

page 23
APPENDIX A
Contents of the Simulation Summary

This appendix is a brief summary of the information contained in the simulation summary screen
displays produced by the flyby_matlab software.

The simulation summary screen display contains the following information:

calendar date = calendar date of trajectory event

TDB time = TDB time of trajectory event

TDB julian date = julian date of trajectory event on TDB time scale

sma (au) = semimajor axis in astronomical unit

eccentricity = orbital eccentricity (non-dimensional)

inclination (deg) = orbital inclination in degrees

argper (deg) = argument of periapsis in degrees

raan (deg) = right ascension of the ascending node in degrees

true anomaly (deg) = true anomaly in degrees

arglat (deg) = argument of latitude in degrees. The argument of latitude is the sum of
true anomaly and argument of perigee.

period (days) = orbital period in days

rx (km) = x-component of the spacecraft’s position vector in kilometers

ry (km) = y-component of the spacecraft’s position vector in kilometers

rz (km) = z-component of the spacecraft’s position vector in kilometers

rmag (km) = scalar magnitude of the spacecraft’s position vector in kilometers

vx (kps) = x-component of the spacecraft’s velocity vector in kilometers per second

vy (kps) = y-component of the spacecraft’s velocity vector in kilometers per second

vz (ksp) = z-component of the spacecraft’s velocity vector in kilometers per second

vmag (kps) = scalar magnitude of the spacecraft’s velocity vector in kilometers per
second

deltav-x = x-component of the impulsive TCM velocity vector in meters/second

deltav-y = y-component of the impulsive TCM velocity vector in meters/second

deltav-z = z-component of the impulsive TCM velocity vector in meters/second

delta-v = scalar magnitude of the impulsive TCM delta-v in meters/seconds

b-magnitude = magnitude of the b-plane vector

b dot r = dot product of the b-vector and r-vector

b dot t = dot product of the b-vector and t-vector

page 24
theta = orientation of the b-plane vector in degrees

v-infinity = magnitude of incoming v-infinity vector in kilometers/second

r-periapsis = periapsis radius of incoming hyperbola in kilometers

decl-asy = declination of incoming v-infinity vector in degrees

rasc-asy = right ascension of incoming v-infinity vector in degrees

fpa = flight path angle in degrees

page 25
APPENDIX B
Classical Orbital Elements of a Hyperbolic Flyby Orbit

This Appendix describes the MATLAB function and mathematical equations used to determine the
classical orbital elements of a hyperbolic flyby orbit at periapsis of the trajectory. The inputs to this
MATLAB function are the incoming and outgoing v-infinity velocity vectors, the periapsis radius of the
flyby hyperbola, and the gravitational constant of the flyby planet.

The syntax of this MATLAB function is

function oev = fbhyper(mu, vinfi, vinfo, rp)

% classical orbital elements of a planet-centered

% flyby hyperbola at periapsis passage

% input

% mu = flyby planet gravitational constant (km**3/sec**2)

% vinfi = incoming v-infinity vector (kilometers/second)
% vinfo = outgoing v-infinity vector (kilometers/second)
% rp = planet-centered periapsis distance (kilometers)

% output

% oev(1) = semimajor axis (kilometers)

% oev(2) = orbital eccentricity (non-dimensional)
% oev(3) = orbital inclination (radians)
% oev(4) = argument of perigee (radians)
% oev(5) = right ascension of ascending node (radians)
% oev(6) = true anomaly (radians)

Algorithm equations

The semimajor axis can be determined from the gravitational constant of the flyby planet and the scalar
magnitude of the flyby v-infinity vector according to the expression


a
v2

The orbital eccentricity of the flyby hyperbola is given by

rp rp v2
e  1  1
a 

A unit vector parallel to the incoming asymptote is given by this next equation

v 
sˆ 
v 

page 26
A unit vector perpendicular to the orbit plane of the flyby hyperbola is calculated from

v   v 
ˆ 
w
v   v 

The orbital inclination relative to the ecliptic plane is determined from the z-component of this unit
vector as follows
i  cos1  wz 

A unit vector along the z-axis of the ecliptic coordinate system is given by

uˆ z   0,1,1
T

From a unit vector in the direction of the ascending node given by

nˆ  uˆ z  w
ˆ

the right ascension of the ascending node can be computed using a four quadrant inverse tangent
function as follows
  tan 1  n y , nx 

A unit vector in the direction of the periapsis of the flyby hyperbola is given by

pˆ   cossˆ  sin bˆ

where
cos   1/ e

sin    1  1/ e2

and bˆ  sˆ  w
ˆ . In these expressions   is the orbital true anomaly at infinity.

The sine and cosine of the argument of periapsis are given by

cos  nˆ  pˆ

sin   1  cos2 

If pz  0 then sin    sin  . Finally, the argument of periapsis is determined from a four quadrant
inverse tangent function as   tan 1 sin  ,cos   .

page 27