P. 1
lambert

# lambert

|Views: 63|Likes:

See more
See less

02/06/2012

pdf

text

original

# Orbital Mechanics with MATLAB

page 1
Lambert’s Problem

This document describes four MATLAB scripts that demonstrate how to solve the Earth orbit,
interplanetary, and
2
J -perturbed form of 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 astrodynamic problem is also known as the orbital two-point boundary value
problem (TPBVP) or the flyby and rendezvous problems.

Lambert’s theorem

The time to traverse a trajectory depends only upon the length of the semimajor axis a of the
transfer trajectory, the sum r r
i f
+ 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:
( )
, ,
i f
tof tof r r c a = +

From the following form of Kepler’s equation

( )
3
0
sin
a
t t E e E
µ
÷ = ÷

we can write
( )
3
0 0
sin sin
a
t E E e E E
µ
= ÷ ÷ ÷ (
¸ ¸

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

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

sin sin 2sin cos
2 2
cos cos 2sin sin
2 2
cos cos 2cos cos
2 2
a
a
a
o | |
o |
o | |
o |
o | |
o |
÷ +
÷ =
÷ +
÷ = ÷
÷ +
+ =

If we let E = o and E
0
= | and substitute the first trig identity into the second equation above,
we have the following equation:

3
0 0
0
2sin cos
2 2
a E E E E
t E E e
µ
÷ + ¦ ¹
| |
= ÷ ÷
´ `
|
\ . ¹ )

Orbital Mechanics with MATLAB
page 2
With the two substitutions given by

0
0
cos cos
2 2
sin sin
2 2
E E
e
E E
o |
o |
+ +
=
÷ ÷
=

the time equation becomes

( )
3
2sin cos
2 2
a
t
o | o |
o |
µ
÷ + ¦ ¹
= ÷ ÷
´ `
¹ )

From the elliptic relationships given by

( )
( )
2
1 cos
cos
sin 1
r a e E
x a E e
y a E e
= ÷
= ÷
= ÷

and some more manipulation, we have the following equations:

0 0
0 0
cos 1 1 1
2 2 2
sin 1 1 1
2 2 2
r r c r r c s
a a a a
r r c r r c s c
a a a a
o
|
+ + +
| |
= ÷ ÷ = ÷ = ÷
|
\ .
+ + ÷ ÷
| |
= ÷ + = ÷ = ÷
|
\ .

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

0
2
0 0
2 2 2 2
0 0
cos cos 1
2 2 2
sin sin sin 1 cos
2 2 2 2
sin sin
2 2 2 2 2
r r
E E E E
e
x x y y c
a a a
o | o |
o | o |
o | o |
÷ + +
= ÷
÷ + ÷ +
| |
= ÷
|
\ .
÷ + ÷ ÷
| | | | | | | |
= + =
| | | |
\ . \ . \ . \ .

With the use of the half angle formulas given by

sin sin
2 2 2 2
s s c
a a
o | ÷
= =
Orbital Mechanics with MATLAB
page 3
and several additional substitutions, we have the time-of-flight form of Lambert’s theorem

( ) ( )
3
sin sin
a
t o | o |
µ
= ÷ ÷ ÷ (
¸ ¸

A discussion about the angles o and | can be found in “Geometrical Interpretation of the
Angles o 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 algorithm used in these MATLAB scripts 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.

Primer Vector Analysis

This section summarizes the primer vector analysis included with the lambert1.m MATLAB
script. The term primer vector was invented by Derek F. Lawden and represents the adjoint
vector for velocity. A technical discussion about primer theory can be found in Lawden’s classic
text, Optimal Trajectories for Space Navigation, Butterworths, London, 1963. Another excellent
resource is “Primer Vector Theory and Applications”, Donald J. Jezewski, NASA TR R-454,
November 1975, along with “Optimal, Multi-burn, Space Trajectories”, also by Jezewski.

As shown by Lawden, the following four necessary conditions must be satisfied in order for an
impulsive orbital transfer to be locally optimal:

(1) the primer vector and its first derivative are everywhere continuous
(2) whenever a velocity impulse occurs, the primer is a unit vector aligned with the impulse
and has unit magnitude
( )
ˆ ˆ and 1
T
= = = p p u p
(3) the magnitude of the primer vector may not exceed unity on a coasting arc
( )
1 p = s p
(4) at all interior impulses (not at the initial or final times) 0 = p p  ; therefore, 0 d dt = p
at the intermediate impulses

Furthermore, the scalar magnitudes of the primer vector derivative at the initial and final
impulses provide information about how to improve the nominal transfer trajectory by changing
the endpoint times and/or moving the impulse times. These four cases for non-zero slopes are
summarized as follows;

- If
0
0 p >  and 0
f
p < ÷  perform an initial coast before the first impulse and add a final
coast after the second impulse
Orbital Mechanics with MATLAB
page 4
- If
0
0 p >  and 0
f
p > ÷  perform an initial coast before the first impulse and move the
second impulse to a later time
- If
0
0 p <  and 0
f
p < ÷  perform the first impulse at an earlier time and add a final coast
after the second impulse
- If
0
0 p <  and 0
f
p > ÷  perform the first impulse at an earlier time and move the second
impulse to a later time

The primer vector analysis of a two impulse orbital transfer involves the following steps.

First partition the two-body state transition matrix as follows:

( )
0 0 11 12
0
21 22
0 0
,
rr rv
vr vv
t t
c c
(
(
c c u u u u
( (
(
u = = =
( (
u u u u c c (
¸ ¸ ¸ ¸
(
c c
¸ ¸
r r
r v
v v
r v

where
0 0 0
11 0 0 0
0
0 0 0
/ / /
/ / /
/ / /
x x x y x z
y x y y y z
z x z y z z
c c c c c c
(
( c
(
u = = c c c c c c
(
(
c
¸ ¸
c c c c c c (
¸ ¸
r
r

and so forth.

The value of the primer vector at any time t along a two body trajectory is given by

( ) ( ) ( )
11 0 0 12 0 0
, , t t t t t = u + u p p p

and the value of the primer vector derivative is

( ) ( ) ( )
21 0 0 22 0 0
, , t t t t t = u + u p p p  

which can also be expressed as

( )
0
0
0
, t t
¦ ¹ ¦ ¹
= u
´ ` ´ `
¹ ) ¹ )
p p
p p  

The primer vector boundary conditions at the initial and final impulses are as follows:

( )
0
0 0
0
t
A
= =
A
V
p p
V

Orbital Mechanics with MATLAB
page 5

( )
f
f f
f
t
A
= =
A
V
p p
V

These two conditions illustrate that at the locations of velocity impulses, the primer vector is a
unit vector in the direction of the impulses.

The value of the primer vector derivative at the initial time is

( ) ( ) ( ) { }
1
0 0 12 0 11 0 0
, ,
f f f
t t t t t
÷
= = u ÷ u p p p p  

provided the
12
u sub-matrix is non-singular.

The scalar magnitude of the derivative of the primer vector can be determined from

( )
2
d d
dt dt
= =
p p p
p p
p
 

lambert1.m – Earth orbit solution

This MATLAB application demonstrates how to solve the two-body form of Lambert’s problem
for a satellite in Earth orbit. The following is a typical user interaction with this script.

program lambert1

< Earth orbit lambert problem >

initial orbit

please input the semimajor axis (kilometers)
(semimajor axis > 0)
? 8000

please input the orbital eccentricity (non-dimensional)
(0 <= eccentricity < 1)
? 0

please input the orbital inclination (degrees)
(0 <= inclination <= 180)
? 28.5

please input the right ascension of the ascending node (degrees)
(0 <= raan <= 360)
? 100

please input the true anomaly (degrees)
(0 <= true anomaly <= 360)
? 0

Orbital Mechanics with MATLAB
page 6
final orbit

please input the semimajor axis (kilometers)
(semimajor axis > 0)
? 8000

please input the orbital eccentricity (non-dimensional)
(0 <= eccentricity < 1)
? 0

please input the orbital inclination (degrees)
(0 <= inclination <= 180)
? 28.5

please input the right ascension of the ascending node (degrees)
(0 <= raan <= 360)
? 100

please input the true anomaly (degrees)
(0 <= true anomaly <= 360)
? 170

please input the transfer time in minutes
? 56

orbital direction

selection (1 or 2)
? 1

please input the maximum number of transfer orbits around the Earth
? 0

The following is the output created by this application.

program lambert1

< Earth orbit lambert problem >

orbital elements of the initial orbit

sma (km) eccentricity inclination (deg) argper (deg)
+8.00000000000000e+003 +0.00000000000000e+000 +2.85000000000000e+001 +0.00000000000000e+000

raan (deg) true anomaly (deg) arglat (deg) period (min)
+1.00000000000000e+002 +0.00000000000000e+000 +0.00000000000000e+000 +1.18684693004297e+002

orbital elements of the transfer orbit after the initial delta-v

sma (km) eccentricity inclination (deg) argper (deg)
+8.00047140990639e+003 +6.70937482986916e-004 +2.85000000000000e+001 +8.49999999999880e+001

raan (deg) true anomaly (deg) arglat (deg) period (min)
Orbital Mechanics with MATLAB
page 7
+1.00000000000000e+002 +2.75000000000012e+002 +0.00000000000000e+000 +1.18695183622590e+002

orbital elements of the transfer orbit prior to the final delta-v

sma (km) eccentricity inclination (deg) argper (deg)
+8.00047140990639e+003 +6.70937482986917e-004 +2.85000000000000e+001 +8.49999999999975e+001

raan (deg) true anomaly (deg) arglat (deg) period (min)
+1.00000000000000e+002 +8.50000000000025e+001 +1.70000000000000e+002 +1.18695183622590e+002

orbital elements of the final orbit

sma (km) eccentricity inclination (deg) argper (deg)
+8.00000000000000e+003 +0.00000000000000e+000 +2.85000000000000e+001 +0.00000000000000e+000

raan (deg) true anomaly (deg) arglat (deg) period (min)
+1.00000000000000e+002 +1.70000000000000e+002 +1.70000000000000e+002 +1.18684693004297e+002

initial delta-v vector and magnitude

x-component of delta-v 0.640619 meters/second
y-component of delta-v -4.677599 meters/second
z-component of delta-v 0.098476 meters/second

delta-v magnitude 4.722290 meters/second

final delta-v vector and magnitude

x-component of delta-v -0.279892 meters/second
y-component of delta-v 4.704815 meters/second
z-component of delta-v -0.293925 meters/second

delta-v magnitude 4.722290 meters/second

total delta-v 9.444579 meters/second

transfer time 56.000000 minutes

The graphical primer vector analysis for this example is shown below. These plots illustrate the
behavior of the scalar magnitudes of the primer vector and its derivative as a function of the orbit
transfer time.
0 500 1000 1500 2000 2500 3000 3500
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Primer Vector Analysis
simulation time (seconds)
p
r
i
m
e
r

v
e
c
t
o
r

m
a
g
n
i
t
u
d
e

Orbital Mechanics with MATLAB
page 8
0 500 1000 1500 2000 2500 3000 3500
−1.5
−1
−0.5
0
0.5
1
1.5
x 10
−3 Primer Vector Analysis
simulation time (seconds)
p
r
i
m
e
r

d
e
r
i
v
a
t
i
v
e

m
a
g
n
i
t
u
d
e

lambert2.m – interplanetary solution

This MATLAB script demonstrates how to solve the two-body interplanetary Lambert problem.
The following is a typical user interaction with this script.

program lambert2

< interplanetary lambert problem >

departure conditions

(1 <= month <= 12, 1 <= day <= 31, year = all digits!)
? 9,1,1998

(0 <= hours <= 24, 0 <= minutes <= 60, 0 <= seconds <= 60)
? 0,0,0

arrival conditions

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

(0 <= hours <= 24, 0 <= minutes <= 60, 0 <= seconds <= 60)
? 0,0,0

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

? 3

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

? 4

The following is the program output for this example.

program lambert2

< interplanetary lambert problem >

departure planet 'Earth'

departure calendar date 01-Sep-1998
departure universal time 00:00:00.000

departure julian date 2451057.500000

arrival planet 'Mars'

arrival calendar date 15-Aug-1999
arrival universal time 00:00:00.000

arrival julian date 2451405.500000

transfer time 348.000000 days

heliocentric ecliptic orbital elements of the departure planet

sma (km) eccentricity inclination (deg) argper (deg)
1.4959802229e+008 1.6709181047e-002 0.0000000000e+000 1.0291439752e+002

raan (deg) true anomaly (deg) arglat (deg) period (days)
0.0000000000e+000 2.3546050380e+002 3.3837490132e+002 3.6525745091e+002

Orbital Mechanics with MATLAB
page 10
heliocentric ecliptic orbital elements of the transfer orbit after the initial
delta-v

sma (km) eccentricity inclination (deg) argper (deg)
1.7145607013e+008 3.3064571581e-001 1.4025406657e+000 8.8020190015e+001

raan (deg) true anomaly (deg) arglat (deg) period (days)
3.3837490132e+002 2.7197980998e+002 0.0000000000e+000 4.4816672004e+002

heliocentric ecliptic orbital elements of the transfer orbit prior to the final
delta-v

sma (km) eccentricity inclination (deg) argper (deg)
1.7145607013e+008 3.3064571581e-001 1.4025406657e+000 8.8020190015e+001

raan (deg) true anomaly (deg) arglat (deg) period (days)
3.3837490132e+002 2.0668220696e+002 2.9470239697e+002 4.4816672004e+002

heliocentric ecliptic orbital elements of the arrival planet

sma (km) eccentricity inclination (deg) argper (deg)
2.2793918413e+008 9.3400274417e-002 1.8497282956e+000 2.8649805839e+002

raan (deg) true anomaly (deg) arglat (deg) period (days)
4.9555144147e+001 2.9704552666e+002 2.2354358506e+002 6.8697161038e+002

initial delta-v vector and magnitude

x-component of delta-v -8563.836300 meters/second
y-component of delta-v 3718.454469 meters/second
z-component of delta-v 729.797537 meters/second

delta-v magnitude 9364.763759 meters/second

energy 87.698800 km^2/sec^2

final delta-v vector and magnitude

x-component of delta-v 4608.052564 meters/second
y-component of delta-v -2097.558698 meters/second
z-component of delta-v -856.066401 meters/second

delta-v magnitude 5134.856434 meters/second

energy 26.366751 km^2/sec^2

After the script computes the numerical data, it will ask the user if he or she would like to create
a graphics display of the trajectory. This prompt appears as

would you like to plot this trajectory (y = yes, n = no)
?

If the user responds with y for yes, the script will ask for the plot step size with

please input the plot step size (days)
?
Orbital Mechanics with MATLAB
page 11
The following is a typical graphics display created with this MATLAB script. The 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 and the arrival planet is labeled with an A. The location of the launch and arrival planets at
the launch time is marked with an asterisk. The plot step size for this example is 5 days.

−2 −1.5 −1 −0.5 0 0.5 1 1.5
−1
−0.5
0
0.5
1
1.5
ϒ
L
A
L
A
Launch
Earth
01−Sep−1998
Arrival
Mars
15−Aug−1999
X coordinate (AU)
Y

c
o
o
r
d
i
n
a
t
e

(
A
U
)
Interplanetary Lambert Problem

lambert3.m – perturbed motion solution – shooting method with state transition matrix updates

This MATLAB script demonstrates how to solve the J
2
-perturbed Earth orbit Lambert problem.
However, more sophisticated equations of motion can easily be implemented. The algorithm
solves this problem using a simple shooting technique.

An initial guess for this algorithm is created by first solving the two-body form of Lambert’s
problem. At each shooting iteration, the initial delta-velocity vector is updated according to

| |
1
12
÷
A = u A v r

where the error in the final position vector Ar is determined from the difference between the
two body final position vector r
tb
and the final position vector predicted by numerical integration
r
int
of the orbital equations of motion as follows:

int tb
A = ÷ r r r

Orbital Mechanics with MATLAB
page 12
The new initial velocity vector can now be calculated from

1 n n +
= + A v v v

The sub-matrix u
12
of the full state transition matrix is as follows:

0 0 0
12 0 0 0
0
0 0 0
/ / /
/ / /
/ / /
x x x y x z
y x y y y z
z x z y z z
c c c c c c
(
( c
(
u = = c c c c c c
(
(
c
¸ ¸
c c c c c c (
¸ ¸
r
v
  
  
  

This sub-matrix consists of the partial derivatives of the rectangular components of the final
position vector with respect to the initial velocity vector.

The following is a typical user interaction with this script.

program lambert3

< j2 perturbed Earth orbit lambert problem >

initial orbit

please input the semimajor axis (kilometers)
(semimajor axis > 0)
? 8000

please input the orbital eccentricity (non-dimensional)
(0 <= eccentricity < 1)
? 0

please input the orbital inclination (degrees)
(0 <= inclination <= 180)
? 28.5

please input the right ascension of the ascending node (degrees)
(0 <= raan <= 360)
? 100

please input the true anomaly (degrees)
(0 <= true anomaly <= 360)
? 0

final orbit

please input the semimajor axis (kilometers)
(semimajor axis > 0)
? 8000

please input the orbital eccentricity (non-dimensional)
(0 <= eccentricity < 1)
? 0
Orbital Mechanics with MATLAB
page 13

please input the orbital inclination (degrees)
(0 <= inclination <= 180)
? 28.5

please input the right ascension of the ascending node (degrees)
(0 <= raan <= 360)
? 100

please input the true anomaly (degrees)
(0 <= true anomaly <= 360)
? 170

please input the transfer time in minutes
? 56

The following is the program output for this example. Please note that the program displays both
the Keplerian (two body) and perturbed solutions for the transfer orbit.

program lambert3

j2 perturbed Earth orbit lambert problem

shooting method with state transition matrix updates

orbital elements of the initial orbit

sma (km) eccentricity inclination (deg) argper (deg)
+8.00000000000000e+003 +0.00000000000000e+000 +2.85000000000000e+001 +0.00000000000000e+000

raan (deg) true anomaly (deg) arglat (deg) period (min)
+1.00000000000000e+002 +0.00000000000000e+000 +0.00000000000000e+000 +1.18684693788431e+002

orbital elements of the final orbit

sma (km) eccentricity inclination (deg) argper (deg)
+8.00000000000000e+003 +0.00000000000000e+000 +2.85000000000000e+001 +0.00000000000000e+000

raan (deg) true anomaly (deg) arglat (deg) period (min)
+1.00000000000000e+002 +1.70000000000000e+002 +1.70000000000000e+002 +1.18684693788431e+002

keplerian transfer orbit

sma (km) eccentricity inclination (deg) argper (deg)
+8.00047141376779e+003 +6.70942937076629e-004 +2.85000000000000e+001 +8.50000000000164e+001

raan (deg) true anomaly (deg) arglat (deg) period (min)
+1.00000000000000e+002 +2.74999999999984e+002 +0.00000000000000e+000 +1.18695184492725e+002

j2 perturbed transfer orbit

sma (km) eccentricity inclination (deg) argper (deg)
+8.00723489628623e+003 +9.68258357992949e-004 +2.89460861696644e+001 +2.10862530931085e+001

raan (deg) true anomaly (deg) arglat (deg) period (min)
+1.00000000000000e+002 +3.38913746906891e+002 +0.00000000000000e+000 +1.18845731080656e+002

delta-v vector and magnitude

Orbital Mechanics with MATLAB
page 14
x-component of delta-v 23.6892 meters/second
y-component of delta-v 1.6813 meters/second
z-component of delta-v 49.7371 meters/second

total delta-v 55.1161 meters/second

transfer time 56.0000 minutes

final position vector error components and magnitude

x-component of delta-r 0.00000999 meters
y-component of delta-r 0.00000171 meters
z-component of delta-r 0.00000460 meters

delta-r magnitude 0.00001113 meters

lambert4.m – perturbed motion solution – NLP solution

This MATLAB application demonstrates how to solve the J
2
-perturbed Earth orbit Lambert
problem. However, more sophisticated equations of motion can easily be implemented. The
algorithm solves this problem using a nonlinear programming technique. This script can solve
both the flyby and rendezvous problems. For the flyby problem, the program attempts to match
all three components of the position vector. For the rendezvous problem, the script attempts to
match all three components of both the target position and velocity vectors.

SNOPT algorithm implementation

This section provides details about the part of the lambert4 script that solves this nonlinear
programming (NLP) problem using the SNOPT 6.0 algorithm. In this classic trajectory
optimization problem, the components of the initial and final delta-v vectors are the control
variables and the scalar magnitude of the flyby or rendezvous V A is the objective function or
performance index.

MATLAB versions of SNOPT 6.0 for several computer platforms can be found at Professor
Philip Gill’s web site which is located at http://cam.ucsd.edu/~peg/Software.html.

The SNOPT algorithm requires an initial guess for the control variables. For this example they
are determined from the two-body solution of Lambert’s problem. The algorithm also requires
lower and upper bounds for the control variables. These are determined from the initial guesses
as follows:

% define lower and upper bounds for components of delta-v vectors
(kilometers/second)

for i = 1:1:3
xlwr(i) = min(-1.1 * norm(xg(1:3)), -75.0);

xupr(i) = max(+1.1 * norm(xg(1:3)), +75.0);
end

if (otype == 2)
for i = 4:1:6
xlwr(i) = min(-1.1 * norm(xg(4:6)), -75.0);
Orbital Mechanics with MATLAB
page 15

xupr(i) = max(+1.1 * norm(xg(4:6)), +75.0);
end
end

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;

Finally, the NLP algorithm also requires the following state vector equality constraints.

% enforce final position vector equality constraints

flow(2) = 0.0d0;
fupp(2) = 0.0d0;

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

flow(4) = 0.0d0;
fupp(4) = 0.0d0;

if (otype == 2)
% enforce final velocity vector equality constraints

flow(5) = 0.0d0;
fupp(5) = 0.0d0;

flow(6) = 0.0d0;
fupp(6) = 0.0d0;

flow(7) = 0.0d0;
fupp(7) = 0.0d0;
end

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

[x, f, inform, xmul, fmul] = snopt(xg, xlwr, xupr, flow, fupp, 'tpbvp');

where tpbvp is the name of the MATLAB function that solves Lambert’s problem and
computes the current value of the objective function and equality constraints. The following is
the MATLAB source code for this example.

function [f, g] = tpbvp(x)

% two point boundary value objective function
% and state vector constraints

Orbital Mechanics with MATLAB
page 16
% input

% x = current delta-v vector

% output

% f(1) = objective function (delta-v magnitude)
% f(2) = rx constraint delta
% f(3) = ry constraint delta
% f(4) = rz constraint delta
% f(5) = vx constraint delta
% f(6) = vy constraint delta
% f(7) = vz constraint delta

% Orbital Mechanics with Matlab

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

global otype neq tetol

global ri vi tof rtarget vtarget drf dvf

% load current state vector of transfer orbit

xi(1) = ri(1);
xi(2) = ri(2);
xi(3) = ri(3);

xi(4) = vi(1) + x(1);
xi(5) = vi(2) + x(2);
xi(6) = vi(3) + x(3);

% initial guess for step size (seconds)

h = 10.0;

% initial time (seconds)

ti = 0.0;

% final time (seconds)

tf = tof;

% integrate equations of motion

xf = rkf78 ('j2eqm', neq, ti, tf, h, tetol, xi);

% objective function (delta-v magnitude)

if (otype == 1)
% initial delta-v only (flyby)

Orbital Mechanics with MATLAB
page 17
f(1) = norm(x);
else
% total delta-v (rendezvous)

f(1) = norm(x(1:3)) + norm(x(4:6));
end

% final position vector equality constraints

f(2) = rtarget(1) - xf(1);

f(3) = rtarget(2) - xf(2);

f(4) = rtarget(3) - xf(3);

if (otype == 2)
% final velocity vector

vf(1) = xf(4) + x(4);
vf(2) = xf(5) + x(5);
vf(3) = xf(6) + x(6);
end

if (otype == 2)
% enforce final velocity vector constraints

f(5) = vtarget(1) - vf(1);

f(6) = vtarget(2) - vf(2);

f(7) = vtarget(3) - vf(3);
end

% save state vector deltas for print summary

for i = 1:1:3
drf(i) = f(i + 1);

if (otype == 2)
% rendezvous

dvf(i) = f(i + 4);
end
end

% transpose objective function/constraints vector

f = f';

% no derivatives

g = [];

Orbital Mechanics with MATLAB
page 18
The following is a typical user interaction with this MATLAB script.

program lambert4

< perturbed Earth orbit Lambert problem >

trajectory type (1 = flyby, 2 = rendezvous)
? 1

classical orbital elements of the initial orbit

please input the semimajor axis (kilometers)
(semimajor axis > 0)
? 8000

please input the orbital eccentricity (non-dimensional)
(0 <= eccentricity < 1)
? 0

please input the orbital inclination (degrees)
(0 <= inclination <= 180)
? 28.5

please input the right ascension of the ascending node (degrees)
(0 <= raan <= 360)
? 100

please input the true anomaly (degrees)
(0 <= true anomaly <= 360)
? 0

classical orbital elements of the final orbit

please input the semimajor axis (kilometers)
(semimajor axis > 0)
? 8000

please input the orbital eccentricity (non-dimensional)
(0 <= eccentricity < 1)
? 0

please input the orbital inclination (degrees)
(0 <= inclination <= 180)
? 28.5

please input the right ascension of the ascending node (degrees)
(0 <= raan <= 360)
? 100

please input the true anomaly (degrees)
(0 <= true anomaly <= 360)
? 170

please input the transfer time in minutes
? 56
Orbital Mechanics with MATLAB
page 19
The following is the script output for this example. The first part of the display includes the two-
body solution for the initial guess and the optimization summary from SNOPT.

two-body guess for initial delta-v vector and magnitude

x-component of delta-v 0.640625 meters/second
y-component of delta-v -4.677637 meters/second
z-component of delta-v 0.098476 meters/second

delta-v magnitude 4.722328 meters/second

Major Minors Step nCon Feasible Optimal MeritFunction nS Penalty
0 2 1 2.8E+01 1.5E-02 4.7223280E-03 r
1 0 3.6E-03 2 2.8E+01 1.5E-02 4.7179225E-03 n r
2 0 1.4E-03 3 2.8E+01 5.4E+00 6.4112054E-03 3.0E-06 s
3 0 1.0E+00 4 1.6E+00 4.0E-02 3.7031413E-02 8.9E-04
4 0 1.0E+00 5 1.2E-03 7.7E-05 5.5116083E-02 1.2E-02 m
5 0 1.0E+00 6 (2.6E-09)(1.2E-10) 5.5116073E-02 1.2E-02 n c
5 0 1.0E+00 6 (2.6E-09)(1.2E-10) 5.5116073E-02 1.2E-02 n c

EXIT -- optimal solution found

program lambert4

< perturbed Earth orbit Lambert problem >

orbital elements and state vector of the initial orbit
------------------------------------------------------

sma (km) eccentricity inclination (deg) argper (deg)
+8.00000000000000e+003 +0.00000000000000e+000 +2.85000000000000e+001 +0.00000000000000e+000

raan (deg) true anomaly (deg) arglat (deg) period (min)
+1.00000000000000e+002 +0.00000000000000e+000 +0.00000000000000e+000 +1.18684693788431e+002

rx (km) ry (km) rz (km) rmag (km)
-1.38918542133544e+003 +7.87846202409766e+003 +0.00000000000000e+000 +8.00000000000000e+003

vx (kps) vy (kps) vz (kps) vmag (kps)
-6.10905247177800e+000 -1.07719077733820e+000 +3.36811407992746e+000 +7.05868645918807e+000

orbital elements and state vector of the transfer orbit after the initial delta-v
---------------------------------------------------------------------------------

sma (km) eccentricity inclination (deg) argper (deg)
+8.00723489628470e+003 +9.68258357783872e-004 +2.89460861700788e+001 +2.10862530925057e+001

raan (deg) true anomaly (deg) arglat (deg) period (min)
+1.00000000000000e+002 +3.38913746907494e+002 +0.00000000000000e+000 +1.18845731080622e+002

rx (km) ry (km) rz (km) rmag (km)
-1.38918542133544e+003 +7.87846202409766e+003 +0.00000000000000e+000 +8.00000000000000e+003

vx (kps) vy (kps) vz (kps) vmag (kps)
-6.08536330087357e+000 -1.07550945881513e+000 +3.41785116756283e+000 +7.06187465926933e+000

orbital elements and state vector of the transfer orbit prior to the final delta-v
----------------------------------------------------------------------------------

sma (km) eccentricity inclination (deg) argper (deg)
+8.00712131217672e+003 +9.65333285374587e-004 +2.89453516511177e+001 +1.47239268484716e+002

raan (deg) true anomaly (deg) arglat (deg) period (min)
+9.98377817600496e+001 +2.29029908765670e+001 +1.70142259361283e+002 +1.18843202316575e+002

Orbital Mechanics with MATLAB
page 20
rx (km) ry (km) rz (km) rmag (km)
+1.65787953981183e+002 -7.97076711063437e+003 +6.62861993417606e+002 +7.99999999999834e+003

vx (kps) vy (kps) vz (kps) vmag (kps)
+6.20553203642353e+000 -1.53599271890635e-001 -3.36706800686018e+000 +7.06182466181549e+000

orbital elements and state vector of the final orbit
----------------------------------------------------

sma (km) eccentricity inclination (deg) argper (deg)
+7.99999999999668e+003 +4.93755404440110e-013 +2.85000000000013e+001 +0.00000000000000e+000

raan (deg) true anomaly (deg) arglat (deg) period (min)
+1.00000000000015e+002 +1.69999999999961e+002 +1.69999999999961e+002 +1.18684693788357e+002

rx (km) ry (km) rz (km) rmag (km)
+1.65787953981183e+002 -7.97076711063437e+003 +6.62861993417606e+002 +7.99999999999834e+003

vx (kps) vy (kps) vz (kps) vmag (kps)
+6.22908767828918e+000 -1.46280648235324e-001 -3.31694485894215e+000 +7.05868645918807e+000

initial delta-v vector and magnitude
------------------------------------

x-component of delta-v 23.689171 meters/second
y-component of delta-v 1.681319 meters/second
z-component of delta-v 49.737088 meters/second

delta-v magnitude 55.116073 meters/second

final position vector error components and magnitude
----------------------------------------------------

x-component of delta-r 0.00000271 meters
y-component of delta-r -0.00000181 meters
z-component of delta-r -0.00000246 meters

delta-r magnitude 0.00000409 meters

transfer time 56.000000 minutes

Here’s the rendezvous option solution for the same initial conditions.

two-body guess for initial delta-v vector and magnitude

x-component of delta-v 0.640625 meters/second
y-component of delta-v -4.677637 meters/second
z-component of delta-v 0.098476 meters/second

delta-v magnitude 6.678380 meters/second

two-body guess for final delta-v vector and magnitude

x-component of delta-v -0.279895 meters/second
y-component of delta-v 4.704854 meters/second
z-component of delta-v -0.293927 meters/second

delta-v magnitude 4.722328 meters/second

total delta-v 9.444656 meters/second

Major Minors Step nCon Feasible Optimal MeritFunction nS Penalty
0 2 1 2.8E+01 1.5E-02 9.4446559E-03 r
1 0 3.2E-03 2 2.8E+01 1.5E-02 9.4374421E-03 n r
2 0 1.4E-03 3 2.8E+01 9.2E-01 1.2614431E-02 5.5E-06 s
Orbital Mechanics with MATLAB
page 21
3 0 1.0E+00 4 1.5E+00 9.6E-02 7.2963247E-02 1.8E-03
4 0 1.0E+00 5 1.2E-03 3.1E-04 1.1097988E-01 2.5E-02 m
5 0 1.0E+00 6 (4.9E-09)(1.5E-09) 1.1097984E-01 2.5E-02 n c
5 0 1.0E+00 6 (4.9E-09)(1.5E-09) 1.1097984E-01 2.5E-02 n c

EXIT -- optimal solution found

program lambert4

< perturbed Earth orbit Lambert problem >

orbital elements and state vector of the initial orbit
------------------------------------------------------

sma (km) eccentricity inclination (deg) argper (deg)
+8.00000000000000e+003 +0.00000000000000e+000 +2.85000000000000e+001 +0.00000000000000e+000

raan (deg) true anomaly (deg) arglat (deg) period (min)
+1.00000000000000e+002 +0.00000000000000e+000 +0.00000000000000e+000 +1.18684693788431e+002

rx (km) ry (km) rz (km) rmag (km)
-1.38918542133544e+003 +7.87846202409766e+003 +0.00000000000000e+000 +8.00000000000000e+003

vx (kps) vy (kps) vz (kps) vmag (kps)
-6.10905247177800e+000 -1.07719077733820e+000 +3.36811407992746e+000 +7.05868645918807e+000

orbital elements and state vector of the transfer orbit after the initial delta-v
---------------------------------------------------------------------------------

sma (km) eccentricity inclination (deg) argper (deg)
+8.00723489628503e+003 +9.68258357800548e-004 +2.89460861700790e+001 +2.10862530883712e+001

raan (deg) true anomaly (deg) arglat (deg) period (min)
+1.00000000000000e+002 +3.38913746911629e+002 +0.00000000000000e+000 +1.18845731080629e+002

rx (km) ry (km) rz (km) rmag (km)
-1.38918542133544e+003 +7.87846202409766e+003 +0.00000000000000e+000 +8.00000000000000e+003

vx (kps) vy (kps) vz (kps) vmag (kps)
-6.08536330087376e+000 -1.07550945881474e+000 +3.41785116756291e+000 +7.06187465926947e+000

orbital elements and state vector of the transfer orbit prior to the final delta-v
----------------------------------------------------------------------------------

sma (km) eccentricity inclination (deg) argper (deg)
+8.00712131217702e+003 +9.65333285316214e-004 +2.89453516511178e+001 +1.47239268486159e+002

raan (deg) true anomaly (deg) arglat (deg) period (min)
+9.98377817600496e+001 +2.29029908751038e+001 +1.70142259361262e+002 +1.18843202316582e+002

rx (km) ry (km) rz (km) rmag (km)
+1.65787953978664e+002 -7.97076711063496e+003 +6.62861993419040e+002 +7.99999999999899e+003

vx (kps) vy (kps) vz (kps) vmag (kps)
+6.20553203642318e+000 -1.53599271892845e-001 -3.36706800685980e+000 +7.06182466181505e+000

orbital elements and state vector of the final orbit
----------------------------------------------------

sma (km) eccentricity inclination (deg) argper (deg)
+7.99999999998343e+003 +2.05098256580623e-012 +2.84999999999968e+001 +0.00000000000000e+000

raan (deg) true anomaly (deg) arglat (deg) period (min)
+1.00000000000018e+002 +1.69999999999938e+002 +1.69999999999938e+002 +1.18684693788062e+002

rx (km) ry (km) rz (km) rmag (km)
+1.65787953978664e+002 -7.97076711063496e+003 +6.62861993419040e+002 +7.99999999999899e+003
Orbital Mechanics with MATLAB
page 22

vx (kps) vy (kps) vz (kps) vmag (kps)
+6.22908767828381e+000 -1.46280648236246e-001 -3.31694485893854e+000 +7.05868645918165e+000

initial delta-v vector and magnitude
------------------------------------

x-component of delta-v 23.689171 meters/second
y-component of delta-v 1.681319 meters/second
z-component of delta-v 49.737088 meters/second

delta-v magnitude 78.476379 meters/second

final delta-v vector and magnitude
----------------------------------

x-component of delta-v 23.555642 meters/second
y-component of delta-v 7.318624 meters/second
z-component of delta-v 50.123148 meters/second

delta-v magnitude 55.863767 meters/second

total delta-v 110.979840 meters/second

final position vector error components and magnitude
----------------------------------------------------

x-component of delta-r 0.00000523 meters
y-component of delta-r -0.00000123 meters
z-component of delta-r -0.00000390 meters

delta-r magnitude 0.00000664 meters

final velocity vector error components and magnitude
----------------------------------------------------

x-component of delta-v 0.00000001 meters/second
y-component of delta-v 0.00000000 meters/second
z-component of delta-v -0.00000000 meters/second

delta-v magnitude 0.00000001 meters/second

transfer time 56.000000 minutes

Lambert Functions

This section describes two MATLAB functions that solve the two-body form of Lambert’s
boundary value problem.

glambert.m – Gooding’s solution of Lambert’s problem

This two-body Lambert function has the following syntax.

function [vi, vf] = glambert(cbmu, sv1, sv2, tof, nrev)

% Gooding's solution of Lambert's problem

% input

Orbital Mechanics with MATLAB
page 23
% cbmu = central body gravitational constant
% sv1 = initial 6-element state vector (position + velocity)
% sv2 = final 6-element state vector (position + velocity)
% nrev = number of full revolutions
% (positive for long period orbit,
% negative for short period orbit)

% output

% vi = initial velocity vector of the transfer orbit
% vf = final velocity vector of the transfer orbit

lambfunc.m – Geodeon’s solution of Lambert’s problem

The algorithm used in this MATLAB function is based on the method described in “A Practical
Note on the Use of Lambert’s Equation” by Geza Gedeon, AIAA Journal, Volume 3, Number 1,
1965, pages 149-150. 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. Additional information can also be found in G. S. Gedeon, “Lambertian
Mechanics”, Proceedings of the 12
th
International Astronautical Congress, Vol. I, 172-190.

The elliptic form of the general Lambert Theorem is

( ) ( ) ( )
3
1 sin sin
a
t k m k t o o | |
µ
= ÷ + ÷ ÷ (
¸ ¸

where k may be either +1 (posigrade) or –1 (retrograde), and m is the number of revolutions

The Gedeon algorithm introduces the following variable

2
s
z
a
=

and solves the problem with a Newton-Raphson procedure. In this equation, a is the semimajor
axis of the transfer orbit and

1 2
2
r r c
s
+ +
=

This algorithm also makes use of the following constant:

1
c
w
s
= ± ÷

The function to be solved iteratively is given by:
Orbital Mechanics with MATLAB
page 24
( ) ( ) ( )
1/ 2
1/ 2 1/ 2 1/ 2 1/ 2 1/ 2 1/ 2
2
1/ 2
1/ 2
1 1
1 1
2
2
k
N z m k z z z w z w z w z w z
z z
t
÷
¦ ¹
(
(
= + ÷ ÷ ÷ ÷ ÷ ÷
´ `
( ¸ ¸
¸ ¸
¹ )

The Newton-Raphson algorithm also requires the derivative of this equation given by

( )
( )
( )
( )
3
1/ 2 1/ 2 1/ 2 1/ 2
2
3 1
2 2
1
1
N z dN k w
N z
dz z
z
w z
¦ ¹
¦ ¦
' = = ÷ ÷
´ `
÷
÷
¦ ¦
¹ )

The iteration for z is as follows:

( )
( )
1
n
n n
n
N z
z z
N z
+
= ÷
'

This Lambert function has the following syntax.

function [statev, nsol] = lambfunc(ri, rf, tof, direct, revmax)

% solve Lambert's orbital two point boundary value problem

% input

% ri = initial ECI position vector (kilometers)
% rf = final ECI position vector (kilometers)
% tof = time of flight (seconds)
% revmax = maximum number of complete orbits

% output

% nsol = number of solutions
% statev = matrix of state vector solutions of the
% transfer trajectory after the initial delta-v

% statev(1, sn) = position vector x component
% statev(2, sn) = position vector y component
% statev(3, sn) = position vector z component
% statev(4, sn) = velocity vector x component
% statev(5, sn) = velocity vector y component
% statev(6, sn) = velocity vector z component
% statev(7, sn) = semimajor axis
% statev(8, sn) = orbital eccentricity
% statev(9, sn) = orbital inclination
% statev(10, sn) = argument of perigee
% statev(11, sn) = right ascension of the ascending node
% statev(12, sn) = true anomaly

% where sn is the solution number

Please note the value of the central body gravitational constant (mu) should be passed to this
function with a global statement located in the main MATLAB script.

Orbital Mechanics with MATLAB

With the two substitutions given by

 E  E0  cos 2 2  E  E0  sin sin 2 2
e cos

the time equation becomes
t a3    cos      2sin   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 s  r  r0  c cos    1   1  1  a 2a  2a 2a  r  r0  c sc  r  r0  c sin    1   1  1  a 2 a  2a 2a  This part of the derivation makes use of the following three relationships:
cos


2

cos


2

 1

r  r0 2
2

sin

 
2

sin


2

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

     x  x0   y  y0   c     sin  sin        2 2   2a   2a   2 a  
With the use of the half angle formulas given by
sin

2

2

s 2a

sin

2

sc 2a

page 2

1979. the primer is a unit vector aligned with the impulse ˆ ˆ and has unit magnitude  p  p  uT and p  1 (3) the magnitude of the primer vector may not exceed unity on a coasting arc p  p  1  (4) at all interior impulses (not at the initial or final times) pp  0 . the scalar magnitudes of the primer vector derivative at the initial and final impulses provide information about how to improve the nominal transfer trajectory by changing the endpoint times and/or moving the impulse times. Donald J. The term primer vector was invented by Derek F. November 1975. Jezewski. Primer Vector Analysis This section summarizes the primer vector analysis included with the lambert1. we have the time-of-flight form of Lambert’s theorem t a3       sin   sin      A discussion about the angles  and  can be found in “Geometrical Interpretation of the Angles  and  in Lambert’s Problem” by J. 1963. Volume 2. Multi-burn. Butterworths. London. E. NASA TR R-454. Optimal Trajectories for Space Navigation. Celestial Mechanics and Dynamical Astronomy 48: 145-165. Lawden and represents the adjoint vector for velocity. Another excellent resource is “Primer Vector Theory and Applications”. d p dt  0 at the intermediate impulses Furthermore. Number 5. Prussing. These four cases for non-zero slopes are summarized as follows. Sept. 1990.m MATLAB script. This iterative solution is valid for elliptic. H. Gooding. parabolic and hyperbolic transfer orbits which may be either posigrade or retrograde. and involve one or more revolutions about the central body. The algorithm used in these MATLAB scripts is based on the method described in “A Procedure for the Solution of Lambert’s Orbital Boundary-Value Problem” by R. As shown by Lawden. along with “Optimal. therefore. pages 442-443. the following four necessary conditions must be satisfied in order for an impulsive orbital transfer to be locally optimal: (1) the primer vector and its first derivative are everywhere continuous (2) whenever a velocity impulse occurs. Space Trajectories”. also by Jezewski.Orbital Mechanics with MATLAB and several additional substitutions. A technical discussion about primer theory can be found in Lawden’s classic text.    If p0  0 and p f  0  perform an initial coast before the first impulse and add a final coast after the second impulse page 3 .-Oct. AIAA Journal of Guidance and Control.

t0  p 0   22  t . First partition the two-body state transition matrix as follows:  r  r 0   t . t0  p 0 and the value of the primer vector derivative is   p  t    21  t .Orbital Mechanics with MATLAB      If p0  0 and p f  0  perform an initial coast before the first impulse and move the second impulse to a later time   If p0  0 and p f  0  perform the first impulse at an earlier time and add a final coast after the second impulse   If p0  0 and p f  0  perform the first impulse at an earlier time and move the second impulse to a later time The primer vector analysis of a two impulse orbital transfer involves the following steps. The value of the primer vector at any time t along a two body trajectory is given by  p  t   11  t . t0    v   r  0 where  x / x0  r   11     y / x0   r0   z / x 0  x / y0 y / y0 z / y0 x / z0  y / z0   z / z0   r  v0   11 12    rr   v   21  22   vr   v0   rv   vv   and so forth. t0  p 0 which can also be expressed as p0  p       t . t0      p  p0  The primer vector boundary conditions at the initial and final impulses are as follows: p  t0   p0  V0 V0 page 4 . t0  p0  12  t .

the primer vector is a unit vector in the direction of the impulses. t0  p0   provided the 12 sub-matrix is non-singular.5 please input the right ascension of the ascending node (degrees) (0 <= raan <= 360) ? 100 please input the true anomaly (degrees) (0 <= true anomaly <= 360) ? 0 page 5 . program lambert1 < Earth orbit lambert problem > initial orbit please input the semimajor axis (kilometers) (semimajor axis > 0) ? 8000 please input the orbital eccentricity (non-dimensional) (0 <= eccentricity < 1) ? 0 please input the orbital inclination (degrees) (0 <= inclination <= 180) ? 28. The following is a typical user interaction with this script. The value of the primer vector derivative at the initial time is 1   p  t0   p0  12  t f .Orbital Mechanics with MATLAB p t f   p f  V f V f These two conditions illustrate that at the locations of velocity impulses. The scalar magnitude of the derivative of the primer vector can be determined from d p  d p p 2   p p   dt dt p lambert1.m – Earth orbit solution This MATLAB application demonstrates how to solve the two-body form of Lambert’s problem for a satellite in Earth orbit. t0  p f  11  t f .

00000000000000e+000 period (min) +1.85000000000000e+001 arglat (deg) +0.49999999999880e+001 period (min) page 6 .70937482986916e-004 true anomaly (deg) inclination (deg) +2.Orbital Mechanics with MATLAB final orbit please input the semimajor axis (kilometers) (semimajor axis > 0) ? 8000 please input the orbital eccentricity (non-dimensional) (0 <= eccentricity < 1) ? 0 please input the orbital inclination (degrees) (0 <= inclination <= 180) ? 28.00000000000000e+000 argper (deg) +0.00047140990639e+003 raan (deg) eccentricity +6.00000000000000e+002 eccentricity +0.00000000000000e+000 inclination (deg) +2.00000000000000e+000 true anomaly (deg) +0. program lambert1 < Earth orbit lambert problem > orbital elements of the initial orbit sma (km) +8.18684693004297e+002 orbital elements of the transfer orbit after the initial delta-v sma (km) +8.85000000000000e+001 arglat (deg) argper (deg) +8.00000000000000e+003 raan (deg) +1.5 please input the right ascension of the ascending node (degrees) (0 <= raan <= 360) ? 100 please input the true anomaly (degrees) (0 <= true anomaly <= 360) ? 170 please input the transfer time in minutes ? 56 orbital direction <1> posigrade <2> retrograde selection (1 or 2) ? 1 please input the maximum number of transfer orbits around the Earth ? 0 The following is the output created by this application.

640619 -4.70937482986917e-004 true anomaly (deg) +8.279892 4.7 0.293925 4.85000000000000e+001 arglat (deg) +1.18695183622590e+002 orbital elements of the final orbit sma (km) +8.9 0.704815 -0.00000000000000e+002 +2.8 primer vector magnitude 0.50000000000025e+001 inclination (deg) +2.00047140990639e+003 raan (deg) +1.00000000000000e+002 eccentricity +6.00000000000000e+000 +1.098476 4.2 0.722290 9.444579 56.6 0.Orbital Mechanics with MATLAB +1.18684693004297e+002 initial delta-v vector and magnitude x-component of delta-v y-component of delta-v z-component of delta-v delta-v magnitude 0.85000000000000e+001 arglat (deg) +1.677599 0.4 0.00000000000000e+002 eccentricity +0. Primer Vector Analysis 1 0.18695183622590e+002 orbital elements of the transfer orbit prior to the final delta-v sma (km) +8.75000000000012e+002 +0.00000000000000e+000 period (min) +1. These plots illustrate the behavior of the scalar magnitudes of the primer vector and its derivative as a function of the orbit transfer time.49999999999975e+001 period (min) +1.722290 meters/second meters/second meters/second meters/second final delta-v vector and magnitude x-component of delta-v y-component of delta-v z-component of delta-v delta-v magnitude total delta-v transfer time -0.70000000000000e+002 argper (deg) +8.5 0.70000000000000e+002 inclination (deg) +2.00000000000000e+003 raan (deg) +1.000000 meters/second meters/second meters/second meters/second meters/second minutes The graphical primer vector analysis for this example is shown below.3 0.00000000000000e+000 true anomaly (deg) +1.70000000000000e+002 argper (deg) +0.1 0 500 1000 1500 2000 2500 3000 3500 simulation time (seconds) page 7 .

1998 please input the universal time (0 <= hours <= 24. year = all digits!) ? 9.0 arrival conditions please input the calendar date (1 <= month <= 12. 1 <= day <= 31. 0 <= minutes <= 60. program lambert2 < interplanetary lambert problem > departure conditions please input the calendar date (1 <= month <= 12. 0 <= minutes <= 60.1999 please input the universal time (0 <= hours <= 24.1.5 x 10 −3 1 primer derivative magnitude 0. 1 <= day <= 31.5 −1 −1.5 0 −0.5 0 500 1000 1500 2000 2500 3000 3500 simulation time (seconds) lambert2.0.m – interplanetary solution This MATLAB script demonstrates how to solve the two-body interplanetary Lambert problem.Orbital Mechanics with MATLAB Primer Vector Analysis 1.0 planetary menu <1> Mercury page 8 . 0 <= seconds <= 60) ? 0. year = all digits!) ? 8. The following is a typical user interaction with this script.15.0. 0 <= seconds <= 60) ? 0.

6709181047e-002 inclination (deg) 0.6525745091e+002 true anomaly (deg) arglat (deg) 2.0000000000e+000 argper (deg) 1.Orbital Mechanics with MATLAB <2> <3> <4> <5> <6> <7> <8> <9> Venus Earth Mars Jupiter Saturn Uranus Neptune Pluto please select the departure planet ? 3 planetary menu <1> <2> <3> <4> <5> <6> <7> <8> <9> Mercury Venus Earth Mars Jupiter Saturn Uranus Neptune Pluto please select the arrival planet ? 4 The following is the program output for this example. program lambert2 < interplanetary lambert problem > departure planet departure calendar date departure universal time departure julian date arrival planet arrival calendar date arrival universal time arrival julian date transfer time 'Earth' 01-Sep-1998 00:00:00.3546050380e+002 3.000000 days heliocentric ecliptic orbital elements of the departure planet sma (km) 1.0000000000e+000 eccentricity 1.0291439752e+002 period (days) 3.500000 'Mars' 15-Aug-1999 00:00:00.000 2451405.000 2451057.4959802229e+008 raan (deg) 0.500000 348.3837490132e+002 page 9 .

8697161038e+002 true anomaly (deg) arglat (deg) 2.698800 meters/second meters/second meters/second meters/second km^2/sec^2 final delta-v vector and magnitude x-component of delta-v y-component of delta-v z-component of delta-v delta-v magnitude energy 4608.052564 -2097.8020190015e+001 period (days) 4.836300 3718. it will ask the user if he or she would like to create a graphics display of the trajectory.8020190015e+001 period (days) 4.7197980998e+002 0.4816672004e+002 true anomaly (deg) arglat (deg) 2.2354358506e+002 initial delta-v vector and magnitude x-component of delta-v y-component of delta-v z-component of delta-v delta-v magnitude energy -8563.3064571581e-001 inclination (deg) 1.066401 5134.856434 26.3400274417e-002 inclination (deg) 1.3837490132e+002 eccentricity 3. the script will ask for the plot step size with please input the plot step size (days) ? page 10 .4025406657e+000 argper (deg) 8.Orbital Mechanics with MATLAB heliocentric ecliptic orbital elements of the transfer orbit after the initial delta-v sma (km) 1.0668220696e+002 2. n = no) ? If the user responds with y for yes.2793918413e+008 raan (deg) 4.4816672004e+002 true anomaly (deg) arglat (deg) 2.4025406657e+000 argper (deg) 8.3837490132e+002 eccentricity 3.454469 729.366751 meters/second meters/second meters/second meters/second km^2/sec^2 After the script computes the numerical data.797537 9364.7145607013e+008 raan (deg) 3.763759 87.9470239697e+002 heliocentric ecliptic orbital elements of the arrival planet sma (km) 2.8649805839e+002 period (days) 6. This prompt appears as would you like to plot this trajectory (y = yes.0000000000e+000 heliocentric ecliptic orbital elements of the transfer orbit prior to the final delta-v sma (km) 1.9704552666e+002 2.3064571581e-001 inclination (deg) 1.9555144147e+001 eccentricity 9.558698 -856.7145607013e+008 raan (deg) 3.8497282956e+000 argper (deg) 2.

The plot is a north ecliptic view where we are looking down on the ecliptic plane from the north celestial pole.5 0 0. The algorithm solves this problem using a simple shooting technique. the initial delta-velocity vector is updated according to v  12  r 1 where the error in the final position vector r is determined from the difference between the two body final position vector rtb and the final position vector predicted by numerical integration rint of the orbital equations of motion as follows: r  rtb  rint page 11 . An initial guess for this algorithm is created by first solving the two-body form of Lambert’s problem. The plot step size for this example is 5 days.5 0 ϒ L −0. the launch planet is labeled with an L and the arrival planet is labeled with an A. The vernal equinox direction is the labeled line pointing to the right.5 X coordinate (AU) lambert3. At each shooting iteration.Orbital Mechanics with MATLAB The following is a typical graphics display created with this MATLAB script.5 −1 −0. However.5 1 1.m – perturbed motion solution – shooting method with state transition matrix updates This MATLAB script demonstrates how to solve the J2 -perturbed Earth orbit Lambert problem. Interplanetary Lambert Problem 1. more sophisticated equations of motion can easily be implemented.5 A Launch Earth 01−Sep−1998 −1 Arrival Mars 15−Aug−1999 A −2 −1.5 L 1 Y coordinate (AU) 0. The location of the launch and arrival planets at the launch time is marked with an asterisk.

5 please input the right ascension of the ascending node (degrees) (0 <= raan <= 360) ? 100 please input the true anomaly (degrees) (0 <= true anomaly <= 360) ? 0 final orbit please input the semimajor axis (kilometers) (semimajor axis > 0) ? 8000 please input the orbital eccentricity (non-dimensional) (0 <= eccentricity < 1) ? 0 page 12 . program lambert3 < j2 perturbed Earth orbit lambert problem > initial orbit please input the semimajor axis (kilometers) (semimajor axis > 0) ? 8000 please input the orbital eccentricity (non-dimensional) (0 <= eccentricity < 1) ? 0 please input the orbital inclination (degrees) (0 <= inclination <= 180) ? 28. The following is a typical user interaction with this script.Orbital Mechanics with MATLAB The new initial velocity vector can now be calculated from v n 1  v n  v The sub-matrix 12 of the full state transition matrix is as follows:   x / x0  r    12     y / x0  v 0   z / x 0    x / y0  y / y0  z / y0  x / z0   y / z0     z / z0  This sub-matrix consists of the partial derivatives of the rectangular components of the final position vector with respect to the initial velocity vector.

00000000000000e+000 inclination (deg) +2. Please note that the program displays both the Keplerian (two body) and perturbed solutions for the transfer orbit.50000000000164e+001 period (min) +1.18684693788431e+002 orbital elements of the final orbit sma (km) +8.68258357992949e-004 true anomaly (deg) +3.70000000000000e+002 argper (deg) +0.89460861696644e+001 arglat (deg) +0. program lambert3 j2 perturbed Earth orbit lambert problem shooting method with state transition matrix updates orbital elements of the initial orbit sma (km) +8.18845731080656e+002 delta-v vector and magnitude page 13 .00000000000000e+002 eccentricity +6.00000000000000e+000 argper (deg) +2.18684693788431e+002 keplerian transfer orbit sma (km) +8.00000000000000e+003 raan (deg) +1.5 please input the right ascension of the ascending node (degrees) (0 <= raan <= 360) ? 100 please input the true anomaly (degrees) (0 <= true anomaly <= 360) ? 170 please input the transfer time in minutes ? 56 The following is the program output for this example.85000000000000e+001 arglat (deg) +0.00000000000000e+000 argper (deg) +0.74999999999984e+002 inclination (deg) +2.00000000000000e+000 argper (deg) +8.38913746906891e+002 inclination (deg) +2.00723489628623e+003 raan (deg) +1.10862530931085e+001 period (min) +1.00000000000000e+002 eccentricity +0.00000000000000e+000 true anomaly (deg) +0.00000000000000e+000 period (min) +1.85000000000000e+001 arglat (deg) +0.00000000000000e+002 eccentricity +9.00000000000000e+002 eccentricity +0.00000000000000e+003 raan (deg) +1.00000000000000e+000 true anomaly (deg) +1.85000000000000e+001 arglat (deg) +1.00047141376779e+003 raan (deg) +1.00000000000000e+000 period (min) +1.Orbital Mechanics with MATLAB please input the orbital inclination (degrees) (0 <= inclination <= 180) ? 28.18695184492725e+002 j2 perturbed transfer orbit sma (km) +8.70942937076629e-004 true anomaly (deg) +2.70000000000000e+002 inclination (deg) +2.

0).0).0 for several computer platforms can be found at Professor Philip Gill’s web site which is located at http://cam.m – perturbed motion solution – NLP solution This MATLAB application demonstrates how to solve the J2 -perturbed Earth orbit Lambert problem.html. the script attempts to match all three components of both the target position and velocity vectors. The algorithm also requires lower and upper bounds for the control variables. the program attempts to match all three components of the position vector.1 * norm(xg(1:3)).ucsd.00000171 0. page 14 .1 * norm(xg(1:3)). -75. xupr(i) = max(+1. MATLAB versions of SNOPT 6.Orbital Mechanics with MATLAB x-component of delta-v y-component of delta-v z-component of delta-v total delta-v transfer time 23.6892 1. However.0 algorithm. -75.7371 55. This script can solve both the flyby and rendezvous problems. SNOPT algorithm implementation This section provides details about the part of the lambert4 script that solves this nonlinear programming (NLP) problem using the SNOPT 6.0000 meters/second meters/second meters/second meters/second minutes final position vector error components and magnitude x-component of delta-r y-component of delta-r z-component of delta-r delta-r magnitude 0. +75. These are determined from the initial guesses as follows: % define lower and upper bounds for components of delta-v vectors (kilometers/second) for i = 1:1:3 xlwr(i) = min(-1.00001113 meters meters meters meters lambert4. For the flyby problem. In this classic trajectory optimization problem. more sophisticated equations of motion can easily be implemented.0).edu/~peg/Software. The algorithm solves this problem using a nonlinear programming technique. The SNOPT algorithm requires an initial guess for the control variables.6813 49. For this example they are determined from the two-body solution of Lambert’s problem.1 * norm(xg(4:6)).00000460 0. end if (otype == 2) for i = 4:1:6 xlwr(i) = min(-1.1161 56. the components of the initial and final delta-v vectors are the control variables and the scalar magnitude of the flyby or rendezvous V is the objective function or performance index.00000999 0. For the rendezvous problem.

0d0. g] = tpbvp(x) % two point boundary value objective function % and state vector constraints page 15 .0d0.0d0. flow(7) = 0. if (otype == 2) % enforce final velocity vector equality constraints flow(5) = 0. fupp(7) = 0. flow(6) = 0. Finally. fupp(4) = 0.Orbital Mechanics with MATLAB xupr(i) = max(+1.0d0.1 * norm(xg(4:6)). fupp(5) = 0. For this problem these bounds are given by % bounds on objective function flow(1) = 0. end The actual call to the SNOPT MATLAB interface function is as follows [x. flow(4) = 0. % enforce final position vector equality constraints flow(2) = 0. fmul] = snopt(xg. inform. flow(3) = 0.0d0.0d0.0d0. The following is the MATLAB source code for this example. xupr. +75. fupp(6) = 0. the NLP algorithm also requires the following state vector equality constraints. function [f. flow. where tpbvp is the name of the MATLAB function that solves Lambert’s problem and computes the current value of the objective function and equality constraints.0d0. xlwr. end end The algorithm requires lower and upper bounds on the objective function.0d0. fupp. xmul.0). fupp(3) = 0.0d0. fupp(2) = 0. 'tpbvp'). fupp(1) = +Inf.0d0.0d0.0d0. f.

xi(6) = vi(3) + x(3). xi(5) = vi(2) + x(2). % initial time (seconds) ti = 0. neq. xi). % integrate equations of motion xf = rkf78 ('j2eqm'.0. ti. % objective function (delta-v magnitude) if (otype == 1) % initial delta-v only (flyby) page 16 . tetol. % initial guess for step size (seconds) h = 10. xi(2) = ri(2). h.Orbital Mechanics with MATLAB % input % x = current delta-v vector % output % % % % % % % f(1) f(2) f(3) f(4) f(5) f(6) f(7) = = = = = = = objective function (delta-v magnitude) rx constraint delta ry constraint delta rz constraint delta vx constraint delta vy constraint delta vz constraint delta % Orbital Mechanics with Matlab %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% global otype neq tetol global ri vi tof rtarget vtarget drf dvf % load current state vector of transfer orbit xi(1) = ri(1). xi(4) = vi(1) + x(1). tf.0. % final time (seconds) tf = tof. xi(3) = ri(3).

xf(1). end end % transpose objective function/constraints vector f = f'. f(4) = rtarget(3) .vf(1). end if (otype == 2) % enforce final velocity vector constraints f(5) = vtarget(1) . end % final position vector equality constraints f(2) = rtarget(1) .vf(2).xf(2). end % save state vector deltas for print summary for i = 1:1:3 drf(i) = f(i + 1). page 17 . % no derivatives g = []. if (otype == 2) % rendezvous dvf(i) = f(i + 4). vf(3) = xf(6) + x(6). vf(2) = xf(5) + x(5).xf(3).vf(3). f(7) = vtarget(3) . f(6) = vtarget(2) . else % total delta-v (rendezvous) f(1) = norm(x(1:3)) + norm(x(4:6)). f(3) = rtarget(2) . if (otype == 2) % final velocity vector vf(1) = xf(4) + x(4).Orbital Mechanics with MATLAB f(1) = norm(x).

0E+00 1.098476 4.00000000000000e+000 true anomaly (deg) +0.0E+00 1.00000000000000e+003 vmag (kps) +7.06187465926933e+000 orbital elements and state vector of the transfer orbit prior to the final delta-v ---------------------------------------------------------------------------------sma (km) +8.8E+01 5.6E-03 1.65333285374587e-004 true anomaly (deg) +2.07550945881513e+000 inclination (deg) +2.29029908765670e+001 inclination (deg) +2.2E-02 n nCon Feasible Optimal 1 2.00000000000000e+000 rz (km) +0.00000000000000e+000 vz (kps) +3.89453516511177e+001 arglat (deg) +1.00723489628470e+003 raan (deg) +1.87846202409766e+003 vy (kps) -1.87846202409766e+003 vy (kps) -1.5E-02 2 2.41785116756283e+000 argper (deg) +2.6E-09)(1.00000000000000e+000 vz (kps) +3.0E-06 s 8.8E+01 1.85000000000000e+001 arglat (deg) +0.00000000000000e+000 ry (km) +7.36811407992746e+000 argper (deg) +0.8E+01 1.2E-10) c c EXIT -.98377817600496e+001 eccentricity +9.4E-03 1.677637 0.38913746907494e+002 ry (km) +7.00000000000000e+003 vmag (kps) +7.5E-02 3 2.2E-10) 6 (2.18843202316575e+002 page 19 .7223280E-03 4.00712131217672e+003 raan (deg) +9.6E-09)(1.5116073E-02 nS Penalty r n r 3.38918542133544e+003 vx (kps) -6.00000000000000e+003 raan (deg) +1.0E+00 1.00000000000000e+002 rx (km) -1.4112054E-03 3.00000000000000e+002 rx (km) -1.2E-03 7.0E-02 5 1.0E+00 0.00000000000000e+000 period (min) +1.2E-02 n 1.89460861700788e+001 arglat (deg) +0.18684693788431e+002 rmag (km) +8.47239268484716e+002 period (min) +1.08536330087357e+000 eccentricity +9.722328 meters/second meters/second meters/second meters/second MeritFunction 4.10862530925057e+001 period (min) +1. The first part of the display includes the twobody solution for the initial guess and the optimization summary from SNOPT. two-body guess for initial delta-v vector and magnitude x-component of delta-v y-component of delta-v z-component of delta-v delta-v magnitude Major Minors 0 2 1 0 2 0 3 0 4 0 5 0 5 0 Step 3.4E+00 4 1.05868645918807e+000 orbital elements and state vector of the transfer orbit after the initial delta-v --------------------------------------------------------------------------------sma (km) +8.68258357783872e-004 true anomaly (deg) +3.2E-02 m 1.Orbital Mechanics with MATLAB The following is the script output for this example.7179225E-03 6.7E-05 6 (2.70142259361283e+002 argper (deg) +1.optimal solution found program lambert4 < perturbed Earth orbit Lambert problem > orbital elements and state vector of the initial orbit -----------------------------------------------------sma (km) +8.5116073E-02 5.18845731080622e+002 rmag (km) +8.7031413E-02 5.6E+00 4.00000000000000e+000 rz (km) +0.9E-04 1.5116083E-02 5.640625 -4.38918542133544e+003 vx (kps) -6.07719077733820e+000 inclination (deg) +2.10905247177800e+000 eccentricity +0.

00000000000000e+000 period (min) +1.69999999999961e+002 ry (km) -7.62861993417606e+002 vz (kps) -3.5E-06 s page 20 .737088 55.5E-02 9.8E+01 meters/second meters/second meters/second meters/second meters/second Optimal 1.36706800686018e+000 rmag (km) +7.000000 meters meters meters meters minutes Here’s the rendezvous option solution for the same initial conditions.00000246 0.62861993417606e+002 vz (kps) -3.293927 4.69999999999961e+002 rz (km) +6.8E+01 3 2.2E-01 MeritFunction 9.678380 meters/second meters/second meters/second meters/second two-body guess for final delta-v vector and magnitude x-component of delta-v y-component of delta-v z-component of delta-v delta-v magnitude total delta-v Major Minors 0 2 1 0 2 0 Step 3.677637 0.53599271890635e-001 rz (km) +6.97076711063437e+003 vy (kps) -1.2E-03 1.681319 49.8E+01 2 2.704854 -0.722328 9.2614431E-02 nS Penalty r n r 5.65787953981183e+002 vx (kps) +6.05868645918807e+000 initial delta-v vector and magnitude -----------------------------------x-component of delta-v y-component of delta-v z-component of delta-v delta-v magnitude 23.46280648235324e-001 inclination (deg) +2.098476 6.4374421E-03 1.279895 4.4E-03 -0.99999999999834e+003 vmag (kps) +7.00000409 56.31694485894215e+000 argper (deg) +0.444656 nCon Feasible 1 2.00000181 -0.85000000000013e+001 arglat (deg) +1.99999999999834e+003 vmag (kps) +7.00000000000015e+002 rx (km) +1.22908767828918e+000 eccentricity +4.4446559E-03 9.Orbital Mechanics with MATLAB rx (km) +1.116073 meters/second meters/second meters/second meters/second final position vector error components and magnitude ---------------------------------------------------x-component of delta-r y-component of delta-r z-component of delta-r delta-r magnitude transfer time 0.06182466181549e+000 orbital elements and state vector of the final orbit ---------------------------------------------------sma (km) +7.20553203642353e+000 ry (km) -7.18684693788357e+002 rmag (km) +7.99999999999668e+003 raan (deg) +1.97076711063437e+003 vy (kps) -1.93755404440110e-013 true anomaly (deg) +1. two-body guess for initial delta-v vector and magnitude x-component of delta-v y-component of delta-v z-component of delta-v delta-v magnitude 0.640625 -4.00000271 -0.5E-02 1.689171 1.65787953981183e+002 vx (kps) +6.

41785116756291e+000 argper (deg) +2.99999999999899e+003 vmag (kps) +7.07719077733820e+000 inclination (deg) +2.2963247E-02 5 1.8E-03 2.00000000000000e+000 period (min) +1.00000000000000e+002 rx (km) -1.00712131217702e+003 raan (deg) +9.85000000000000e+001 arglat (deg) +0.20553203642318e+000 eccentricity +9.9E-09)(1.97076711063496e+003 inclination (deg) +2.62861993419040e+002 vz (kps) -3.1097984E-01 6 (4.69999999999938e+002 ry (km) -7.98377817600496e+001 rx (km) +1.65787953978664e+002 vx (kps) +6.00000000000000e+000 vz (kps) +3.00000000000000e+000 vz (kps) +3.87846202409766e+003 vy (kps) -1.18684693788431e+002 rmag (km) +8.optimal solution found program lambert4 < perturbed Earth orbit Lambert problem > orbital elements and state vector of the initial orbit -----------------------------------------------------sma (km) +8.89460861700790e+001 arglat (deg) +0.00000000000000e+003 raan (deg) +1.00000000000000e+003 vmag (kps) +7.1097984E-01 1.47239268486159e+002 period (min) +1.5E-02 n 2.38918542133544e+003 vx (kps) -6.36811407992746e+000 argper (deg) +0.70142259361262e+002 rz (km) +6.1E-04 1.5E+00 9.00000000000000e+003 vmag (kps) +7.2E-03 3.05868645918807e+000 orbital elements and state vector of the transfer orbit after the initial delta-v --------------------------------------------------------------------------------sma (km) +8.97076711063496e+003 vy (kps) -1.68258357800548e-004 true anomaly (deg) +3.0E+00 1.00000000000000e+000 rz (km) +0.99999999998343e+003 raan (deg) +1.00000000000000e+002 rx (km) -1.18845731080629e+002 rmag (km) +8.10862530883712e+001 period (min) +1.65787953978664e+002 eccentricity +2.87846202409766e+003 vy (kps) -1.99999999999899e+003 page 21 .0E+00 1.29029908751038e+001 ry (km) -7.6E-02 7.53599271892845e-001 inclination (deg) +2.06187465926947e+000 orbital elements and state vector of the transfer orbit prior to the final delta-v ---------------------------------------------------------------------------------sma (km) +8.36706800685980e+000 argper (deg) +1.5E-09) 1.62861993419040e+002 argper (deg) +0.10905247177800e+000 eccentricity +0.9E-09)(1.00000000000000e+000 true anomaly (deg) +0.0E+00 4 1.00000000000000e+000 rz (km) +0.00000000000018e+002 rx (km) +1.1097988E-01 6 (4.Orbital Mechanics with MATLAB 3 4 5 5 0 0 0 0 1.5E-02 n c c EXIT -.07550945881474e+000 inclination (deg) +2.5E-09) 1.84999999999968e+001 arglat (deg) +1.38918542133544e+003 vx (kps) -6.06182466181505e+000 orbital elements and state vector of the final orbit ---------------------------------------------------sma (km) +7.89453516511178e+001 arglat (deg) +1.18684693788062e+002 rmag (km) +7.00000000000000e+000 ry (km) +7.69999999999938e+002 rz (km) +6.0E+00 1.5E-02 m 2.08536330087376e+000 eccentricity +9.18843202316582e+002 rmag (km) +7.00723489628503e+003 raan (deg) +1.05098256580623e-012 true anomaly (deg) +1.65333285316214e-004 true anomaly (deg) +2.00000000000000e+000 period (min) +1.38913746911629e+002 ry (km) +7.

689171 1. glambert.979840 meters/second final position vector error components and magnitude ---------------------------------------------------x-component of delta-r y-component of delta-r z-component of delta-r delta-r magnitude 0.000000 minutes Lambert Functions This section describes two MATLAB functions that solve the two-body form of Lambert’s boundary value problem.123148 55.00000123 -0.476379 meters/second meters/second meters/second meters/second final delta-v vector and magnitude ---------------------------------x-component of delta-v y-component of delta-v z-component of delta-v delta-v magnitude 23. vf] = glambert(cbmu. sv2.863767 meters/second meters/second meters/second meters/second total delta-v 110. nrev) % Gooding's solution of Lambert's problem % input page 22 .00000000 -0.00000390 0.Orbital Mechanics with MATLAB vx (kps) +6.555642 7.46280648236246e-001 vz (kps) -3. sv1.00000001 meters/second meters/second meters/second meters/second transfer time 56.00000523 -0.00000001 0.m – Gooding’s solution of Lambert’s problem This two-body Lambert function has the following syntax. function [vi.22908767828381e+000 vy (kps) -1.00000664 meters meters meters meters final velocity vector error components and magnitude ---------------------------------------------------x-component of delta-v y-component of delta-v z-component of delta-v delta-v magnitude 0.00000000 0.05868645918165e+000 initial delta-v vector and magnitude -----------------------------------x-component of delta-v y-component of delta-v z-component of delta-v delta-v magnitude 23. tof.318624 50.737088 78.31694485893854e+000 vmag (kps) +7.681319 49.

Vol. Volume 3. Number 1. 172-190. AIAA Journal. This iterative solution is valid for elliptic. Gedeon.Orbital Mechanics with MATLAB % % % % % % % cbmu sv1 sv2 tof nrev = = = = = central body gravitational constant initial 6-element state vector (position + velocity) final 6-element state vector (position + velocity) time of flight (+ posigrade.m – Geodeon’s solution of Lambert’s problem The algorithm used in this MATLAB function is based on the method described in “A Practical Note on the Use of Lambert’s Equation” by Geza Gedeon.retrograde) number of full revolutions (positive for long period orbit. parabolic and hyperbolic transfer orbits which may be either posigrade or retrograde. I. Proceedings of the 12th International Astronautical Congress. S. The Gedeon algorithm introduces the following variable z s 2a and solves the problem with a Newton-Raphson procedure. 1965. a is the semimajor axis of the transfer orbit and r r c s 1 2 2 This algorithm also makes use of the following constant: w   1 c s The function to be solved iteratively is given by: page 23 . In this equation. pages 149-150. The elliptic form of the general Lambert Theorem is t a3  1  k  m  k   sin       sin      where k may be either +1 (posigrade) or –1 (retrograde). and involve one or more revolutions about the central body. and m is the number of revolutions about the central body. Additional information can also be found in G. . “Lambertian Mechanics”. negative for short period orbit) % output % % vi = initial velocity vector of the transfer orbit vf = final velocity vector of the transfer orbit lambfunc.

sn) = right ascension of the ascending node statev(12. nsol] = lambfunc(ri. sn) = position vector z component statev(4. revmax) % solve Lambert's orbital two point boundary value problem % input % % % % % ri rf tof direct revmax = = = = = initial ECI position vector (kilometers) final ECI position vector (kilometers) time of flight (seconds) transfer direction (1 = posigrade. sn) = semimajor axis statev(8. rf. direct. page 24 . sn) = argument of perigee statev(11. tof. sn) = position vector y component statev(3. sn) = velocity vector y component statev(6.Orbital Mechanics with MATLAB N z  1 zz 1/ 2 21/ 2 1/ 2  1/ 2 1/ 2 1/ 2 1/ 2 1/ 2 1/ 2 1  k m  k  z  z 1  z     w z  w z  w z 1  w2 z            2 The Newton-Raphson algorithm also requires the derivative of this equation given by  3N  z   1  dN k w3  N  z      1/ 2  1/ 2 1/ 2 1/ 2  dz z 2  1  z  2 1  w2 z     The iteration for z is as follows: zn 1  zn  This Lambert function has the following syntax. sn) = orbital inclination statev(10. sn) = true anomaly where sn is the solution number Please note the value of the central body gravitational constant (mu) should be passed to this function with a global statement located in the main MATLAB script. sn) = velocity vector z component statev(7. N  zn  N   zn  function [statev. sn) = velocity vector x component statev(5. sn) = orbital eccentricity statev(9. sn) = position vector x component statev(2. -1 = retrograde) maximum number of complete orbits % output % % % % % % % % % % % % % % % % nsol = number of solutions statev = matrix of state vector solutions of the transfer trajectory after the initial delta-v statev(1.

scribd
/*********** DO NOT ALTER ANYTHING BELOW THIS LINE ! ************/ var s_code=s.t();if(s_code)document.write(s_code)//-->
{"extraction_info": {"found_math": false, "script_math_tex": 0, "script_math_asciimath": 0, "math_annotations": 0, "math_alttext": 0, "mathml": 0, "mathjax_tag": 0, "mathjax_inline_tex": 0, "mathjax_display_tex": 0, "mathjax_asciimath": 0, "img_math": 0, "codecogs_latex": 0, "wp_latex": 0, "mimetex.cgi": 0, "/images/math/codecogs": 0, "mathtex.cgi": 0, "katex": 0, "math-container": 0, "wp-katex-eq": 0, "align": 0, "equation": 0, "x-ck12": 0, "texerror": 0, "math_score": 0.8200035095214844, "perplexity": 27748.91444422919}, "config": {"markdown_headings": true, "markdown_code": true, "boilerplate_config": {"ratio_threshold": 0.18, "absolute_threshold": 10, "end_threshold": 15, "enable": true}, "remove_buttons": true, "remove_image_figures": true, "remove_link_clusters": true, "table_config": {"min_rows": 2, "min_cols": 3, "format": "plain"}, "remove_chinese": true, "remove_edit_buttons": true, "extract_latex": true}, "warc_path": "s3://commoncrawl/crawl-data/CC-MAIN-2016-07/segments/1454701166222.10/warc/CC-MAIN-20160205193926-00087-ip-10-236-182-209.ec2.internal.warc.gz"}