Jenab6

jenab6


Jenab's Fireside Chats


Previous Entry Share Next Entry
Elliptical Transfer Orbits
Jenab6
jenab6
Define the initial orbit of the transfer object and the orbit of the target.

I have shown how to calculate Hohmann transfer orbits. Now, I propose to show how to determine elliptical transfer orbits in a more general way. As an example case, I will find an elliptical transfer orbit beginning at the asteroid 2001-YB5 and ending at an interception of Earth.

2001-YB5, orbital elements
a = 2.349557177836 au
e = 0.8624274715129
i = 5.490700413641 degrees
L = 109.3451209415 degrees
w = 114.2474452629 degrees
T = JD 2453637.57768


For the purpose of this example problem, we will assume that Earth's orbit will be:

Earth, orbital elements
a = 1.0000001124 au
e = 0.0167102192
i = 0
L = 0
w = 103.078101 degrees
T = JD 2454468.667


The orbital elements are the semimajor axis (a), the eccentricity (e), the inclination (i), the longitude of the ascending node (L), the argument of the perihelion (w), and the time of perihelion passage (T). I won't bother using subscript tags for orbital elements; instead, I'll rely on context to keep them distinct.

Choose trial times for departure and for arrival.

The first task is to choose trial times for the departure from the initial orbit and for the arrival of the transfer object at the target. These are the endpoints of the intended trajectory in the transfer orbit. The significant event for the departure time is the rocket burn required to make the appropriate change in velocity (delta-vee) to move from the initial orbit into the transfer orbit. The significant event at the arrival time is either impact or a rendezvous delta-vee.


t1 = time of departure = 18h UT on 29 April 2018 = JD 2458238.25
t2 = time of arrival = 18:28:48 UT on 6 January 2020 = JD 2458855.27


Conversion from Calendar Date to Julian Date
Conversion from Julian Date to Calendar Date

Of course, it is necessary that t2 > t1, and it also helps if the difference t2-t1 is within bounds that might be deemed plausible in some seat-of-the-pants sense. I've already done the necessary hunting for suitable such times.

In our example, the required transit time, t2 - t1 = 617.02 days.

For simplicity, since this is only an illustration of technique, we will ignore the complication of Earth's gravitational acceleration near the arrival end of the transfer orbit.

The next task is to locate the positions of departure and arrival. These will define two points on the transfer orbit. The departure position is calculated from the elements of the initial orbit and the selected trial value for the departure time. The arrival position is calculated from the elements of the target's orbit and the selected trial value for the arrival time.

In our example, the position of departure will be the position of the asteroid at time t1, and the position of arrival will be the position Earth at time t2.

We should define some useful constants:

AU = 1.49597870691E+11 m au^-1
GM = 1.32712440018E+20 m^3 sec^-2
pi = 3.141592653589793238462643
k0 = 0.01745329251994329576923691
k1 = 365.256898326 days

Finding position and velocity in heliocentric ecliptic coordinates from a set of orbital elements and a time.

Find the period, P, of the orbit in days.
P = k1 a^1.5

Find the mean anomaly, m, of the orbit at time t.
m = 2 pi (t - T) / P

Adjust m to the interval [0, 2 pi).

Find the eccentric anomaly. Try Danby's method first.

U1 = m
REPEAT...
.   U0 = U1
.   F0 = U0 - e sin U0 - m
.   F1 = 1 - e cos U0
.   F2 = e sin U0
.   F3 = e cos U0
.   D1 = -F0 / F1
.   D2 = -F0 / [ F1 + D1 F2 / 2 ]
.   D3 = -F0 / [ F1 + D1 F2 / 2 + (D2)^2 F3 / 6 ]
.   U1 = U0 + D3
UNTIL |U1-U0| < 1E-15
u = U1

If Danby's method fails to converge, then try to find the eccentric anomaly with Father Wiggly's Famous Reverse Interpolation Bisection Method.

n = -1
REPEAT...
.   n = n+1
.   U0 = 10 n k0
.   U2 = 10 (n+1) k0
.   M0 = U0 - e sin U0
.   M2 = U2 - e sin U2
.   dM0 = m - M0
.   dM2 = m - M2
UNTIL  { (dM0)(dM2) < = 0 } or {n>35}
if n>35 then STOP {procedure failed}
REPEAT...
.   U1 = (U0+U2) / 2
.   M0 = U0 - e sin U0
.   M1 = U1 - e sin U1
.   dM0 = m - M0
.   dM1 = m - M1
.   if (dM0)(dM1) < 0 then U2=U1 else U0=U1
UNTIL (U2-U0) < 1E-13
M0 = U0 - e sin U0
M2 = U2 - e sin U2
u = U0 + (U2-U0)(m-M0) / (M2-M0)

To get the U's to converge to within 1E-15 or 1E-13 tolerances, your calculator or computer will need to be using variables with double precision or better. If you don't have double precision, you'll need to increase the error tolerance to something your computational means can handle.

Position.

Find the canonical (triple prime) heliocentric position vector.

x''' = a (cos u - e)
y''' = a sin u sqrt (1-e^2)
z''' = 0

Find the true anomaly, Q. We'll use it below when we find the velocity.

Q = Arctan( y''' , x''' )

Arctan(y,x) is the two-dimensional arctangent function, related to the principle value arctangent (atn) as follows:

Function Arctan( y , x )
.   if x = 0 and y > 0 then angle = +pi / 2
.   if x = 0 and y = 0 then angle = 0
.   if x = 0 and y < 0 then angle = -pi / 2
.   if x > 0 and y > 0 then angle = atn( y / x )
.   if x < 0 then angle = atn( y / x ) + pi
.   if x > 0 and y < 0 then angle = atn( y / x ) + 2 pi
Arctan = angle


Rotate the triple-prime position vector by the argument of the perihelion, w.

x'' = x''' cos w - y''' sin w
y'' = x''' sin w + y''' cos w
z'' = z''' = 0

Rotate the double-prime position vector by the inclination, i.

x' = x''
y' = y'' cos i
z' = y'' sin i

Rotate the single-prime position vector by the longitude of the ascending node, L.

x = x' cos L - y' sin L
y = x' sin L + y' cos L
z = z'

The unprimed position vector is the position in heliocentric ecliptic coordinates.

Velocity.

Find the canonical (triple prime) heliocentric velocity vector.

k2 = sqrt { GM / [ a AU (1 - e^2) ] }
k2 is a speed in meters per second.

Vx''' = -k2 sin Q
Vy''' = k2 (e + cos Q)
Vz''' = 0

Rotate the triple-prime velocity vector by the argument of the perihelion, w.

Vx'' = Vx''' cos w - Vy''' sin w
Vy'' = Vx''' sin w + Vy''' cos w
Vz'' = Vz''' = 0

Rotate the double-prime velocity vector by the inclination, i.

Vx' = Vx''
Vy' = Vy'' cos i
Vz' = Vy'' sin i

Rotate the single-prime velocity vector by the longitude of the ascending node, L.

Vx = Vx' cos L - Vy' sin L
Vy = Vx' sin L + Vy' cos L
Vz = Vz'

The unprimed velocity vector is the sun-relative velocity in ecliptic coordinates.

Here are the results for our example:

The position and velocity of 2001-YB5 at 
JD 2458238.25 are
x1 = +3.159148898997291 au
y1 = +3.003558117525086 au
z1 = -0.3821685497977586 au
Vxo = -3565.785981875893 m/s
Vyo = +3891.390270455813 m/s
Vzo = +199.4993435825594 m/s

The position and velocity of Earth at 
JD 2458855.27 are
x2 = -0.2819965365811233 au
y2 = +0.9420187015477031 au
z2 = 0
Vxd = -29022.48342622212 m/s
Vyd = -8655.470317741644 m/s
Vzd = 0


The time difference, t2-t1, is the required transit time for the transfer object in the transfer orbit. This is distinct from the calculated transit time, dt, which has yet to be found. In order for a valid transfer orbit to exist between r1 at t1 and r2 at t2, it is necessary that dt = t2-t1. In general, this will not be the case. However, as I said earlier, I've already done the preliminary work, and the required and calculated transit times will be very nearly equal for the example.

The eccentricity of a hypothetical transfer orbit.

We calculate the distances of the sides of the triangle whose vertices are the sun, the position of departure, and the position of arrival.

r1 = sqrt [ (x1)^2 + (y1)^2 + (z1)^2 ]
r2 = sqrt [ (x2)^2 + (y2)^2 + (z2)^2 ]
d = sqrt [ (x2-x1)^2 + (y2-y1)^2 + (z2-z1)^2 ]

For Hohmann transfers, both departure and arrival are marked by transfer orbit apsides. However, in this more general case, only one of those endpoints is so marked: it must be either, but it is never both.

We d
efine the integer variable K and permit it to have only the values 1 and 2.

If K=1, an apside (perihelion or aphelion) of the transfer orbit occurs at departure.
If K=2, an apside
(perihelion or aphelion) of the transfer orbit occurs at arrival.

K = either 1 or 2
J = 3 - K
N = (-1)^J

If K=1, then J=2 and N=1
If K=2, then J=1 and N=-1

The variables K and J will usually be subscripts. The variable N is a sign toggle factor.

In our example, K=1 and J=2 and N=1.

Now then, if the heliocentric distance of the endpoint having the apside is less than that of the endpoint which does not have the apside, then the apside in question must be the perihelion, and therefore the true anomaly of the apsidal endpoint, QK, must be zero.

On the other hand, if the heliocentric distance of the endpoint having the apside is greater than that of the endpoint which doesn't have the apside, then the apside in question must be the aphelion, and therefore the true anomaly of the apsidal endpoint, QK, must be pi radians.

If rK < rJ then QK = 0 else QK = pi radians

We immediately set the mean anomaly and the eccentric anomaly equal to whatever is the value of the true anomaly, in the case of the apsidal endpoint of the intended trajectory.

mK = uK = QK

In our example, m1 = u1 = Q1 = pi radians, meaning that the aphelion of a hypothetical transfer orbit occurs at the point of departure from the asteroid 2001-YB5 at the time t1
= JD 2458238.25.

The eccentricity of a conic section, having the sun at a focus, which includes the point of departure and the point of arrival, is found by solving, simultaneously, the polar equation which relates the heliocentric distance with the true anomaly and the law of cosines. After some algebra, we get

e = 2 cos(QK) rK ( rK - rJ ) / ( rJ^2 - rK^2 - d^2 )

For our example,
K=1, Q1=pi radians, cos Q1 = -1
r1 = 4.375801175995221 au
r2 = 0.9833215550925033 au
d = 4.029575594635826 au

e = 0.8626144800739287


Conditions for which no elliptical transfer orbit exists.

If any of the following conditions exist, then STOP and try the other value of K.

e = 1

Although parabolic transfer orbits are possible, they are outside the scope of this article

e > 1 and QK = 0

The above condition corresponds to a hyperbolic transfer orbit having its perihelion at either departure (if K=1) or at arrival (if K=2). However, hyperbolic transfer orbits are outside the scope of this article. I have a procedure for finding hyperbolic transfer orbits in another post.

e = 0

Although circular transfer orbits are possible, they are outside the scope of this article.

e < 0

No orbit whatever has an eccentricity less than zero.

e > 1 and QK <> 0

The above condition corresponds to a hyperbolic transfer orbit having its aphelion at either departure or at arrival. However, this requires rK=infinity, which contradicts the finite value calculated for rK. Therefore, this is an impossible condition, and obtaining it means that no such transfer orbit exists.

In any of those circumstances, there exists no elliptical transfer orbit for the value of K being used. It is possible that neither of the values permitted to K (namely 1 and 2) will permit an elliptical transfer orbit. Which means that the values of t1 and t2 that were inappropriately selected. It becomes necessary, in that case, to try other times for t1 or for t2, or for both, and work the math to this point once again.

Furthermore, even if none of the above conditions exists, any orbit corresponding to an otherwise acceptable eccentricity [0 < e < 1] must provide a calculated transit time, dt, that is equal to the required transit time, t2-t1. Until it is known that dt=t2-t1, the orbit is only a hypothetical elliptical transfer orbit.

The true anomaly of the non-apsidal endpoint of the intended trajectory in a hypothetical elliptical transfer orbit.

We find the cosine, c, of the angle subtended at the sun between the position of departure and the position of arrival. That angle is the arc of true anomaly traveled by the transfer object in the hypothetical elliptical transfer orbit between departure and arrival.

c = ( rK^2 + rJ^2 - d^2 ) / ( 2 rK rJ )

Because we'll assign to the arc of true anomaly the principle value of ArcCos(c), all of the transfer orbits found by this procedure will be elliptical transfer orbits of the short path, i.e., elliptical transfer orbits in which the transfer object moves less than 180 degrees of true anomaly from departure to arrival. This is important. There are ways to adjust this procedure for elliptical orbits of the long path, and for super-periodic transfer orbits in either direction, but we will confine ourselves in this article to elliptical transfer orbits of the short path only.

Recall that in our example K=1, J=2, and N=1.

We find the true anomaly of the non-apsidal end of the hypothetical transfer orbit.

QJ = QK + N ArcCos c

The principle value arctangent function can be substituted for the ArcCos, if desired.

QJ = QK + N { (pi/2) - atn[ c / sqrt( 1 - c^2 ) ] }

Adjust QJ to the interval [0, 2 pi), if necessary.

For our example,
Q1 = pi radians
c = 0.4505275708291625
Q2 = 4.245032787432119 radians


The semimajor axis of a hypothetical elliptical transfer orbit.

We find the semimajor axis of the hypothetical transfer orbit. If the eccentricity, e, of the hypothetical transfer orbit is strictly between zero and one [0 < e < 1], then the hypothetical transfer orbit is an ellipse, and the semimajor axis is found from

a = rK / ( 1 - e cos QK )

In our example, a = 2.349279049855524 au.

The eccentric and mean anomalies in the hypothetical elliptical transfer orbit at the non-apsidal endpoint.

The sine, s, and the cosine, c, of the eccentric anomaly in the transfer orbit at the non-apsidal endpoint of the intended trajectory is found from

s = ( rJ / a ) sin QJ / sqrt ( 1-e^2 )
c = e + ( rJ / a ) cos QJ

The eccentric anomaly, u, is found by using the two dimensional arctangent function.

uJ = Arctan ( s , c )

In our example, u2 = 5.452053669284399 radians.

The mean anomaly in the transfer orbit at the non-apsidal endpoint of the intended trajectory is found as the straightforward solution to Kepler's equation.

mJ = uJ - e sin uJ

In our example, m2 = 6.089262339183971 radians.

The calculated transit time in the hypothetical transfer orbit.

At last we will find out whether our departure time (t1) and arrival time (t2) were well chosen. Although we will not usually be continuing along a transfer orbit for a whole period, it is still necessary to determine what that period is.

P = k1 a^1.5

There P is the period of the transfer orbit in days. The mean motion, b, of the transfer object in the transfer orbit is found from

b = 2 pi / P

There b is the mean motion in the transfer orbit in radians per day. The calculated transit time in the transfer orbit, from departure to arrival, is

dt = (N/b) [ mJ - pi sin(QK/2) ]

In our example,
P = 1315.225848439035 days
b = 0.004777267200638389 radians/day
K=1, J=2, N=1
dt = 617.0200580784495 days


Recall that t2-t1 = 617.02 days. The difference between the calculated transit time and the required transit time is only about five seconds. That's close enough. We can drop the "hypothetical" from our language and acknowledge that there does indeed exist a valid transfer orbit from the asteroid 2001-YB5 to Earth having a departure at 18h UT on 29 April 2018 and an arrival at 18:28:48 UT on 6 January 2020.

The time of perihelion passage for the transfer object in the transfer orbit.


Perihelion passage will occur either at departure or at arrival, or else it will not occur at all during the transit of the transfer object along its intended trajectory. It is necessary to know the time of perihelion passage anyway, so that we can generate an ephemeris for the transfer object, and thereby make sure that the transfer object reaches the proper place at the proper time, intercepting the target.

T = tK - mK/b


In our example,
K=1
t1 =
JD 2458238.25
m1 = pi radians
b = 0.004777267200638389 radians/day
T = JD 2457580.637075781

The inclination of the transfer orbit.

We find the cross product r1 x r2 as follows:

Xp' = y1 z2 - z1 y2
Yp' = z1 x2 - x1 z2
Zp' = x1 y2 - y1 x2

The magnitude of the cross product is

Rp = sqrt [ (Xp)^2 + (Yp)^2 + (Zp)^2 ]

We obtain the unit vector normal to the orbit and codirectional with the orbit's angular momentum.

Xp = Xp' / Rp
Yp = Yp' / Rp
Zp = Zp' / Rp

The inclination of the transfer orbit, i, is found from

i = ArcCos (Zp)

The principle value arctangent function can be substituted for the ArcCos, if desired.

i = pi/2 - atn { Zp / sqrt [ 1 - (Zp)^2 ] }

In our example, i = 5.61408792389817 degrees.

Velocity in the elliptical transfer orbit at the apsidal endpoint.

A vector codirectional with the velocity in the transfer orbit at the apsidal endpoint can be found from the cross product Rp x rK.

VxK'' = Yp zK - Zp yK
VyK'' = Zp xK - Xp zK
VzK'' = Xp yK - Yp xK

The magnitude of the cross product is

VK'' = sqrt [ (VxK'')^2 + (VyK'')^2 + (VzK'')^2 ]

The codirectional unit vector is

VxK' = VxK'' / VK''
VyK' = VyK'' / VK''
VzK' = VzK'' / VK''

The actual magnitude of the velocity at the apsidal endpoint is found from the Vis Viva equation.

VK = sqrt[ (GM/AU) ( 2 / rK - 1/a ) ]

The velocity in the elliptical transfer orbit at the apsidal endpoint is

VxK = VK VxK'
VyK = VK VyK'
VzK = VK VzK'

In our example, K=1, meaning that departure is the apsidal endpoint. The position and velocity of the transfer object in the transfer orbit (i.e. immediately after delta-vee) at departure are

x1 = +3.159148898997291 au
y1 = +3.003558117525086 au
z1 = -0.3821685497977586 au
Vx1 = -3618.095915873970 m/s
Vy1 = +3835.117316284865 m/s
Vz1 = +232.6042211888594 m/s


Now that the heliocentric position and the sun-relative velocity in the transfer orbit are both known for the same moment of time, tK, we may be sure that we can fully determine the transfer orbit's elements.

[ xK, yK, zK, VxK, VyK, VzK ] & [ tK ] ---> [ a, e, i, L, w, T ]

Once we have the orbital elements, we can calculate the position and velocity in the transfer orbit at any specified time within a reasonable interval from the elements' epoch.

[ a, e, i, L, w, T ] & [ t ] ---> [ x, y, z, Vx, Vy, Vz ]

We have already found its eccentricity (e), its semimajor axis (a), its time of perihelion passage (T), and its inclination (i). What remains is the longitude of the ascending node (L) and the argument of the perihelion (w).

The longitude of the ascending node of the transfer orbit.

We begin by finding the transfer orbit's angular momentum (per unit mass) from the cross product rK x VK.

hx = AU [ yK VzK - zK VyK ]
hy = AU [ zK VxK - xK VzK ]
hz = AU [ xK VyK - yK VxK ]


In our example,
hx = 3.237748988940674E+14 m^2 sec^-1
hy = 9.692312898885492E+13 m^2 sec^-1
hz = 3.438188115976321E+15 m^2 sec^-1


The longitude of the ascending node, L, is found from

L = Arctan ( hx , -hy )

In our example, L = 106.6652516775637 degrees.

Be wary in the use of the above expression. The quantity hx is playing the role of sin L, while the quantity (-hy) is playing the role of cos L.

The argument of the perihelion of the transfer orbit.

c = ( xK cos L + yK sin L ) / rK

if [ sin i = 0 ] then s = ( yK cos L - xK sin L ) / rK
if [ sin i <> 0 ] then s = zK / ( rK sin i )

w' = Arctan ( s , c )

The quantity w' is the sum of the argument of the perihelion and the true anomaly at the point on the orbit that is being considered, which in this procedure is the apsidal endpoint of the intended trajectory. We already know the true anomaly at this point, so we find the argument of the perihelion from

w = w' - QK

Adjust w to the interval [0, 2 pi), if necessary.


Transfer orbit elements.
Departing from asteroid 2001-YB5 at JD 2458238.25
Arriving at Earth at JD 2458855.27
Transit time 617.02 days.

a = 2.349279049855524 au
e = 0.8626144800739287
i = 5.61408792389817 degrees
L = 106.6652516775637 degrees
w = 116.7775373854853 degrees
T = JD 2457580.637075781

Position and Velocity in the transfer orbit at the non-apsidal endpoint.

Refer to the section headed "Finding position and velocity in heliocentric ecliptic coordinates from a set of orbital elements and a time." We will do now for the transfer object in the transfer orbit exactly what we did when we were calculating the position and velocity of the target at time t2 or those of the transfer object in its initial orbit at t1.

Because dt might be only approximately equal to t2-t1, it will be necessary to recalculate the position of the transfer object in the transfer orbit for time t2, rather than merely assuming that it's position will be identical with that of the target at t2. The velocity of the transfer object must be calculated, of course, to determine the delta-vee at time t2.

In our example,
The non-apsidal endpoint is arrival, so J=2.
The target is Earth.
t2 = JD 2458855.27


Find the canonical (triple prime) heliocentric position vector.

xJ''' = a (cos uJ - e)
yJ''' = a sin uJ sqrt (1-e^2)
zJ''' = 0

Find the true anomaly, Q. We'll use it below when we find the velocity.

QJ (recalculated) = Arctan( yJ''' , xJ''' )

In our example, Q2 (recalculated) = 4.245031986300209 radians. Compare this to Q2 =
4.245032787432119 radians. The error in true anomaly at arrival [Q2 - Q2 (recalculated) <> 0] arises because dt is not exactly equal to t2-t1. If we had worked a bit harder to find values of t2 and t1 that made dt more nearly equal to t2-t1, the error would have been smaller. [Note: if t2 were changed to JD 2458855.26990126, the error in true anomaly would be reduced from about 8E-7 radians to about 3E-12 radians. I'm not going to change t2; I am only showing the dependence of precision in the output on precision in the input.]

Rotate the triple-prime position vector by the argument of the perihelion, w.

xJ'' = xJ''' cos w - yJ''' sin w
yJ'' = xJ''' sin w + yJ''' cos w
zJ'' = zJ''' = 0

Rotate the double-prime position vector by the inclination, i.

xJ' = xJ''
yJ' = yJ'' cos i
zJ' = yJ'' sin i

Rotate the single-prime position vector by the longitude of the ascending node, L.

xJ (recalculated) = xJ' cos L - yJ' sin L
yJ (recalculated) = xJ' sin L + yJ' cos L
zJ (recalculated) = zJ'

The unprimed position vector is the (recalculated) position in heliocentric ecliptic coordinates. It should be checked against the already known values for [xJ,  yJ,  zJ].

In our example,
x2 (recalculated) = -0.2819960700947116 au
y2 (recalculated) = +0.9420198770150876 au
z2 (recalculated) = -0.0000000770657545 au
The error in position, |r2 - r2 (recalculated)|, is about 189.5 kilometers. [Note: if t2 were changed to JD 2458855.26990126, the position error would be reduced to 70 centimeters!]


Find the canonical (triple prime) heliocentric velocity vector.

k2 = sqrt { GM / [ a AU (1 - e^2) ] }
k2 is a speed in meters per second.

VxJ''' = -k2 sin QJ
VyJ''' = k2 (e + cos QJ)
Vz'J'' = 0

Rotate the triple-prime velocity vector by the argument of the perihelion, w.

VxJ'' = VxJ''' cos w - VyJ''' sin w
VyJ'' = VxJ''' sin w + VyJ''' cos w
VzJ'' = VzJ''' = 0

Rotate the double-prime velocity vector by the inclination, i.

VxJ' = VxJ''
VyJ' = VyJ'' cos i
VzJ' = VyJ'' sin i

Rotate the single-prime velocity vector by the longitude of the ascending node, L.

VxJ = VxJ' cos L - VyJ' sin L
VyJ = VxJ' sin L + VyJ' cos L
VzJ = VzJ'

The unprimed velocity vector is the sun-relative velocity in ecliptic coordinates.


In our example,
Vx2 = -13907.07996471122 m/s
Vy2 = -35043.47505289391 m/s
Vz2 = +2297.514387170954 m/s


Checking the work by evolving the state vector at time t1 by time integration, and comparing the result to the state vector at time t2.

We begin with the state vector in the transfer orbit at the time of departure, t1.

t1 = JD 2458238.25
x1 = +3.159148898997291 au = +472601948485.8118 meters
y1 = +3.003558117525086 au = +449325898878.4212 meters
z1 = -0.3821685497977586 au = -57171601294.81209 meters
Vx1 = -3618.095915873970 m/s
Vy1 = +3835.117316284865 m/s
Vz1 = +232.6042211888594 m/s


We integrate this vector with some appropriate time increment from t1 to t2, using a numerical model of the restricted two-body problem.

In the example, t2-t1 = 617.02 days, so we integrate for 53310528 seconds.

And compare the resulting state vector with our already calculated state vector for time t2.

Here is the integrated (step=10 sec) state vector at (approximately) t2.

t2 =
JD 2458855.27
x2 (integrated) = -42186037845.24232 meters = -0.2819962453368013 au
y2 (integrated) = +140924297714.0508 meters = +0.942020745770742 au
z2 (integrated) = -12723.83233543315 meters = -0.000000085053566 au
Vx2 (integrated) = -13907.09059089995 m/s
Vy2 (integrated) = -35043.43295392387 m/s
Vz2 (integrated) = +2297.514201064887 m/s


Here, once again, is the previously calculated state vector at t2.

t2 = JD 2458855.27
x2 = -0.2819965365811233 au
y2 = +0.9420187015477031 au
z2 = 0

Vx2 = -13907.07996471122 m/s
Vy2 = -35043.47505289391 m/s
Vz2 = +2297.514387170954 m/s


It seems fairly close. Evidently, I haven't screwed up yet.


The delta-vee for departure.

We determine the sun-relative departure delta-vee, in ecliptic coordinates, from the vector difference V1-Vo.

dVx1 = Vx1 - Vxo
dVy1 = Vy1 - Vyo
dVz1 = Vz1 - Vzo

Recall that the initial orbit was the orbit of asteroid 2001-YB5 and that at departure our transfer object was sharing the asteroid's position and velocity. The sun-relative velocity, in ecliptic coordinates, of the transfer object in the initial orbit, immediately prior to the departure delta-vee, is

Vxo = -3565.785981875893 m/s
Vyo = +3891.390270455813 m/s
Vzo = +199.4993435825594 m/s

The sun-relative velocity, in ecliptic coordinates, of the transfer object in the transfer orbit, immediately after the departure delta-vee, is

Vx1 = -3618.095915873970 m/s
Vy1 = +3835.117316284865 m/s
Vz1 = +232.6042211888594 m/s


The vector difference V1 - Vo is the departure delta-vee in ecliptic coordinates.

dVx1 = -52.309933998077 m/s
dVy1 = -56.272954170948 m/s
dVz1 = +33.104877606300 m/s


We may rotate the delta-vee from ecliptic coordinates to celestial coordinates. Let k3 be the obliquity of the Earth at the time of departure, t1, expressed in full Julian date. The obliquity in degrees is approximately found from

k3 = 23.439282 - 3.563E-7 (t1 - 2451543.5)

In our example, t1 =
JD 2458238.25, and therefore k3 = 23.436896660575 degrees.

dVx1' = dVx1
dVy1' = dVy1 cos k3 - dVz1 sin k3
dVz1' = dVy1 sin k3 + dVz1 cos k3

The primed delta-vee is in rectangular celestial coordinates. We prefer them in spherical coordinates, so

dV1 = sqrt [ (dVx1')^2 + (dVy1')^2 + (dVz1')^2 ]
RAdv1 = Arctan { dVy1' / dV1 , dVx1' / dV1 }
DECdv1 = ArcSin { dVz1' / dV1 }

In our example, the departure delta-vee has a magnitude of
83.659473 m/s. The direction of thrust is right ascension 15h 24m 20.7902s and declination +5.4816562 degrees. The rockets are aimed toward a direction in the constellation Serpens Caput, about five degrees west and one degree south of Unukalhai (Alpha Serpentis), which is near the galaxy NGC 5921.

[Note: changing t2 to
JD 2458855.26990126 would result in dv1=83.660071 m/s, RAdv1=15h 24m 21.8469s, DECdv1=5.4807962 degrees.]

The delta-vee for arrival.

We determine the sun-relative departure delta-vee, in ecliptic coordinates, from the vector difference Vd-V2.

dVx2 = Vxd - Vx2
dVy2 = Vyd - Vy2
dVz2 = Vzd - Vz2


Vx2 = -13907.07996471122 m/s
Vy2 = -35043.47505289391 m/s
Vz2 = +2297.514387170954 m/s

Vxd = -29022.48342622212 m/s
Vyd = -8655.470317741644 m/s
Vzd = 0

dVx2 = -15115.40346151090 m/s
dVy2 = +26388.00473515226 m/s
dVz2 = -2297.514387170954 m/s

We may rotate the delta-vee from ecliptic coordinates to celestial coordinates. Let k3 (this time) be the obliquity of the Earth at the time of arrival, t2, expressed in full Julian date. The obliquity in degrees is approximately found from

k3 = 23.439282 - 3.563E-7 (t2 - 2451543.5)

In our example, t2 =
JD 2458855.27, and therefore k3 = 23.43667682 degrees.

dVx2' = dVx2
dVy2' = dVy2 cos k3 - dVz2 sin k3
dVz2' = dVy2 sin k3 + dVz2 cos k3


dV2 = sqrt [ (dVx2')^2 + (dVy2')^2 + (dVz2')^2 ]
RAdv2 = Arctan { dVy2' / dV2 , dVx2' / dV2 }
DECdv2 = ArcSin { dVz2' / dV2 }


In our example, the arrival delta-vee (if we are minded to apply one) has a magnitude of 30497.225908 m/s. The direction of thrust is right ascension 8h 4m 7.6051s and declination +15.9636363 degrees. Again, this procedure ignores the transfer object's acceleration by Earth's gravity near arrival. Earth's gravity will accelerate the transfer object so that it has a speed of 32548.45 m/s upon impact with Earth's surface, if no delta-vee is applied.

The main purpose of this procedure is to develop an initial approximation for the departure delta-vee, so that it can be used as an initial guess in a different sort of procedure. That next procedure would be a many-body simulation with kinematic integration over small increments of time, which would account for the motions of the transfer object, the sun, the target object, and every other object of gravitational importance. The time increment would be reduced in size until the resulting position at time t2 converged. Then the initial delta-vee would be very slightly modified in magnitude and in direction, with the resulting changes in position error at t2 being noted, until it became clear which delta-vee was correct.

In Memoriam: Victor G. Szebehely (1921-1997), Hungarian (White), celestial mechanic & engineer, NASA Project Apollo, University of Texas.

Jerry Abbott

i do not know why i'm writing this, but thank you for the this knowledge you are spreading. i don't think i had the capacity to grasp this read, ever, but i know there must be someone out there who benefited. i' ve only been seeing LJs / blogs for art & lit uses but never this scientific things. so yeah, thanks again!

Some background in celestial mechanics is very helpful in understanding my essay on transfer orbits. With that background, you'd be able mentally to "see" what I was doing and, to an extent, to anticipate where I was going. This is the math behind spaceship piloting, and, if mankind ever gets to the point where spaceflight is commonplace, there will be job opportunities for people with this knowledge. By putting the transfer tutorial here, I'm hoping to popularize celestial mechanics and make it at least a little easier to learn than it was before.

You are viewing jenab6