[Jenab]

I'm putting this transfer orbit tutorial here because someone suggested the idea of getting methane from Jupiter or Titan.

The procedure given in the posts to follow will permit the user to calculate the elements of a transfer orbit between a known orbit of departure and a known destination orbit, with a departure at a specified moment of time, in a specified interval of transit time. It will also inform the user of the necessary changes of velocity - or delta-vees - to enter the transfer orbit at departure and to enter the rendezvous orbit at arrival.

This method for finding transfer orbits will treat only ellipses and hyperbolas. Circular and parabolic transfer orbits will not be considered.

The default coordinate system will be heliocentric ecliptic coordinates in which the XY plane (the ecliptic) shall be the current plane of Earth's orbit, with positions being relative to the center of the sun, and with velocities being sun-relative, unless otherwise noted.

All transfer orbits sought by this method will have an apside at one, but not at both, endpoints of the intended trajectory. That is, either departure or arrival will occur at either aphelion or perihelion of the transfer orbit (and both of those instances of "or" are really "xor" - exclusive or).

HELPFUL DESCRIPTIONS.

a : semimajor axis
For elliptical orbits, the semimajor axis is half the length of the longest line segment that can be drawn inside the ellipse. For hyperbolic orbits, the semimajor axis is the distance from the vertex to the nearest point on either branch.

e : eccentricity
The orbit is a circle when e=0. The orbit is an ellipse when 0<e<1. The orbit is a parabola when e=1. The orbit is a hyperbola when e>1.

i : inclination to ecliptic
I prefer to use the principal value of the arc-cosine for my range of inclination [0 , pi radians]. I take the ecliptic to be the current plane of Earth's orbit. (This is the hillbilly tutorial, remember?)

L : longitude of ascending node
The angle, subtended at the sun, measured in the ecliptic, from the Vernal Equinox to the ascending node. The ascending node is the point at which the orbit crosses the ecliptic having the Z component of its velocity in the direction of the North Ecliptic Pole. (Most people use a capital omega for this quantity.)

w : argument of perihelion
The angle, subtended at the sun, measured in the plane of the orbit in the direction of motion, from the ascending node to the orbit's perihelion. (Most people use a lower-case omega for this quantity.)

T : time of perihelion passage
The time (e.g. a decimal calendar or Julian date) at which the orbiting object passes through the perihelion of its orbit.

Q : true anomaly
The angle, subtended at the sun, measured in the plane of the orbit in the direction of motion, from the perihelion to the current position of the orbiting object. No Greek letters can be rendered here, sorry.

M : mean anomaly
The angle, subtended at the sun, measured in the plane of the orbit in the direction of motion, from the perihelion to the point at which the orbiting object would be found, if it moved constantly by its average angular velocity.

u : eccentric anomaly
The angle, subtended at the geometric center of the orbit, measured in the plane of the orbit in the direction of motion, from the perihelion to a point on the circle circumscribing the orbit, which point is found by extending a line perpendicular to the major axis through the current position of the orbiting object until it intersects the circumscribing circle. No pictures, sorry.

WHERE TO FIND ORBITAL ELEMENTS.

Planetary orbital elements (NASA/JPL):
http://ssd.jpl.nasa.gov/elem_planets.html

Except don't believe the values for Earth's inclination and longitude of ascending node. Set them both equal to zero. The other elements in that table are OK. Why they didn't just go ahead and give the times of perihelion passage (or a mean anomaly at epoch), I've no idea.

Planetary orbital elements (some other guy):
http://www.astro.uu.nl/~strous/AA/en/reken/hemelpositie.html

Don't believe his value for Earth's longitude of ascending node, either. Set it to zero.

Asteroid orbital elements:
http://ssd.jpl.nasa.gov/sb_elem.html

Convert a calendar date to a Julian date:
http://wwwmacho.mcmaster.ca/JAVA/JD.html

Convert a Julian date to a calendar date:
http://wwwmacho.mcmaster.ca/JAVA/CD.html

Jerry Abbott

[Jenab]

The unfortunate fact to someone wanting to figure a transfer orbit is being unable to know, ahead of time, where the arrival position is. What you can actually see is where the spaceship and the rendezvous object are at the moment of departure. And without knowing in advance where the arrival will be, you must do some trial and error work in order to ensure that the spaceship and the rendezvous object will reach the arrival position in the same amount of transit time.

But while developing the theory of transfer orbit calculation, it is getting ahead of things to worry about transit time matching. We will do the trial-and-error work later, once the theory is in place for determining the transfer orbit's elements and the transit time in the transfer orbit. For now, we will assume that we know where the arrival will occur.

Let these be the orbital elements of Earth.

a = 1.00000011 AU
e = 0.01671022
i = 0
L = 0
w = 102.94719 degrees
T = JD 2453009.3

Let these be the orbital elements of Vesta.

a = 2.3626478 AU
e = 0.08887781
i = 7.13485 degrees
L = 103.94712 degrees
w = 149.67895 degrees
T = JD 2452941.1

We assert that a spaceship will depart from Vesta at

t1 = 4 February 2004, 19h 12m UT
t1 = JD 2453040.3

We assert that the spaceship will arrive at Earth at

t2 = 16 September 2004, 21h 36m UT
t2 = JD 2453265.4

We thus have established a required transit time for the spaceship, going from Vesta to Earth, of t2-t1, or 225.1 days.

Our job is to find out whether these assertions can possibly be true, astrodynamically speaking. That is, does there exist a valid transfer orbit from the departure position at t1 to the arrival position at t2?

Jerry Abbott

[Jenab]

Given its orbital elements and a specified moment of time, you can find the heliocentric position and sun-relative velocity at that moment.

Find the orbital period.

P = (365.256898326 days) a^1.5

Find the mean motion.

m = 2 pi radians / P

Find the mean anomaly.

M = m (t-T)

Find the eccentric anomaly. (Danby's Method used.)

u(0) = M

Repeat over subscript j,

f0 = u(j) - e sin u(j) - M
f1 = 1 - e cos u(j)
f2 = e sin u(j)
f3 = e cos u(j)
d1 = -f0/f1
d2 = -f0/(f1 + d1 f2 / 2)
d3 = -f0/(f1 + d1 f2 / 2 + d2^2 f3 / 6)
u(j+1) = u(j) + d3

Until | u(j+1) - u(j) | < 1E-12

The eccentric anomaly, u, is the converged value from the loop.

Finding the canonical position vector.

x''' = a [ (cos u) - e ]

y''' = a (1 - e^2)^0.5 sin u

z''' = 0

Rotation by the argument of the perihelion.

x'' = x''' cos w - y''' sin w

y'' = x''' sin w + y''' cos w

z'' = z''' = 0

Rotation by the inclination.

x' = x''

y' = y'' cos i

z' = y'' sin i

Rotation by the longitude of ascending node.

x = x' cos L - y' sin L

y = x' sin L + y' cos L

z = z'

This vector [x,y,z] is the heliocentric position vector at time [t] of a planet having orbital elements [a , e , i , L , w , T].

Finding the true anomaly.

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

The arctan2 function is the two-argument arctangent. It is related to the principal value arctan function as follows:

Function arctan2( y , x );
begin
if x=0 then
begin
if y>0 then Q:=+pi/2
if y=0 then Q:=0
if y<0 then Q:=-pi/2
end else begin
Q:=arctan(y/x)
if x<0 then Q:=Q+pi
if x>0 and y<0 then Q:=Q+2 pi
end
arctan2:=Q
end

Velocity in an elliptical orbit.

You will want to convert the semimajor axis from astronomical units to meters.

1 astronomical unit = 1.49597870691E+11 meters

GMsun = 1.32712440018E+20 m^3 sec^-2

VX’’’ = -sin Q { GMsun / [ a (1-e^2) ] }^0.5

VY’’’ = (e + cos Q) { GMsun / [ a (1-e^2) ] }^0.5

VZ’’’ = 0

The coordinate rotations by w, i, and L proceed as with the position vector.

Putting in the numbers...

I will from now on usually present angles in radians in the interval [0, 2 pi.)

Finding the HEC state vector for Vesta on JD 2453040.3

a = 2.3626478 AU
e = 0.08887781
i = 0.1245266 radians
L = 1.81422 radians
w = 2.612391 radians
T = JD 2452941.1

P = 1326.468 days
m = 0.004736778 radians per day
t-T = 99.2 days
M = 0.4698881 radians
u = 0.5135515 radians

x''' = +1.847896
y''' = +1.156106
z''' = 0

x'' = -2.178777
y'' = -0.06506175
z'' = 0

x' = -2.178777
y' = -0.06455795
z' = -0.008080998

x = +0.5877971
y = -2.098983
z = -0.008080998

Q = 0.5590546 radians

Vx''' = -10318.3 m/s
Vy''' = +18221.6 m/s
Vz''' = 0

...intermediate steps skipped...

Vx = +20234.0 m/s
Vy = +4724.0 m/s
Vz = -2600.6 m/s

Finding the HEC state vector for Earth on JD 2453265.4

a = 1.00000011 AU
e = 0.01671022
i = 0
L = 0
w = 1.796767 radians
T = JD 2453009.3

P = 365.2569 days
m = 0.0172021 radians per day
t-T = 256.1 days
M = 4.405455 radians
u = 4.389608 radians

x''' = -0.3339517
y''' = -0.9482125
z''' = 0

x'' = +0.9989289
y'' = -0.1130118
z'' = 0

x' = +0.9989289
y' = -0.1130118
z' = 0

x = +0.9989289
y = -0.1130118
z = 0

Q = 4.373764 radians

Vx''' = +28097.2 m/s
Vy''' = -9397.8 m/s
Vz''' = 0

...intermediate steps skipped...

Vx = +2863.6 m/s
Vy = +29488.5 m/s
Vz = 0

Summary: Preburn and rendezvous state vectors.

The HEC position and velocity of Vesta at the moment of departure are:

x1 = +0.5877971
y1 = -2.098983
z1 = -0.008080998
Vx(preburn) = +20234.0 m/s
Vy(preburn) = +4724.0 m/s
Vz(preburn) = -2600.6 m/s

The HEC position and velocity of Earth at the moment of arrival are:

x2 = +0.9989289
y2 = -0.1130118
z2 = 0
Vx(rendezvous) = +2863.6 m/s
Vy(rendezvous) = +29488.5 m/s
Vz(rendezvous) = 0

I could have designated the positions with (preburn) or with (rendezvous), too, but they are also equal to the positions of the transit object at departure and at arrival, respectively. For quantities pertaining to a transfer orbit, departure will usually be denoted by the subscript 1, while arrival will usually be denoted by the subscript 2.

Jerry Abbott

[Jenab]

We have two HEC position vectors, R1 and R2, which we suppose to be in the plane of the transfer orbit. The vectors must be distinct and not exactly opposite in direction.

We find the cross product, R1xR2.

Xn' = y1 z2 - z1 y2

Yn' = z1 x2 - x1 z2

Zn' = x1 y2 - y1 x2

We find the magnitude of the primed normal vector, Rn'.

Rn' = [ (Xn')^2 + (Yn')^2 + (Zn')^2 ]^0.5

We divide the primed vector's components each by the magnitude of the primed vector.

Xn = Xn' / Rn'

Yn = Yn' / Rn'

Zn = Zn' / Rn'

The vector [Xn, Yn, Zn] is a unit normal vector with respect to the transfer orbit.

i = pi/2 - arcsin(Zn)

The variable i is the inclination of the transfer orbit plane with respect to the ecliptic.

Plugging in the numbers...

We suppose that these two HEC position vectors are in the transfer orbit plane.

x1 = +0.5877971
y1 = -2.098983
z1 = -0.008080998

x2 = +0.9989289
y2 = -0.1130118
z2 = 0

The cross product gives

Xn' = -9.132484E-4
Yn' = -8.072343E-3
Zn' = +2.030307

The magnitude, Rn', is 2.030324 AU

The unit normal to the transfer orbit is

Xn = -4.498044E-4
Yn = -3.975890E-3
Zn = +0.9999920

The inclination of the transfer orbit,

i = 0.003996730 radians

Jerry Abbott

[Jenab]

We find the lengths of the three sides of a triangle having the sun, the point of departure, and the point of arrival as their vertices.

r1 = [ (x1)^2 + (y1)^2 + (z1)^2 ]^0.5

r1 is the distance from the sun to the point of departure.

r2 = [ (x2)^2 + (y2)^2 + (z2)^2 ]^0.5

r2 is the distance from the sun to the point of arrival.

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

d is the line-of-light separation between the points of departure and arrival.

In order to get the eccentricity and semimajor axis of the transfer orbit from these three values, we must restrict our allowed transfer orbits to those having an apside at either endpoint of the intended trajectory.

That is, perihelion or aphelion must occur at departure or at arrival. (Both instances of "or" in the previous sentence are exclusive. If one endpoint of the intended trajectory has an apside, the other endpoint may not have the other one, or this procedure will fail.)

We proceed by considering, first, the Law of Cosines, understanding that the angle between r1 and r2 is the difference in true anomaly Q2-Q1.

cos(Q2-Q1) = ( r1^2 + r2^2 - d^2 ) / (2 r1 r2)

Secondly, we solve the polar equation of the orbit, with the pole at the solar focus, for the cosine of the true anomaly.

cos Q = [ a (1 - e^2) - r ] / (e r)

There are four general cases for a subsequent solution:

1. Departure occurs from the perihelion of the transfer orbit.
Outward bound trajectory, Q1 = 0.

2. Arrival occurs at the aphelion of the transfer orbit.
Outward bound trajectory, Q2 = pi radians.

3. Arrival occurs at the perihelion of the transfer orbit.
Inward bound trajectory, Q2 = 0.

4. Departure occurs from the aphelion of the transfer orbit.
Inward bound trajectory, Q1 = pi radians.

Now, I'm not about to solve all four cases. I'm just going to solve case #1 explicitly, to provide proof of principle. You can solve the other three cases yourself, if you want to check my work. I'm only going to provide the results for cases #2, #3, and #4.

INSIST Q1=0. Therefore, r1 = a(1-e).

Why? Because Q1=0 means that departure occurs at the transfer orbit's perihelion, and the heliocentric distance at perihelion equals the product of the semimajor axis and of one minus the eccentricity.

We thus have two equations for cos Q2, namely,

cos Q2 = ( r1^2 + r2^2 - d^2 ) / (2 r1 r2)

cos Q2 = [ r1 (1-e^2) / (1-e) - r2 ] / (e r2)

Meaning that...

e r2 ( r1^2 + r2^2 - d^2 ) = 2 r1^2 r2 (1+e) - 2 r1 r2^2

...solve, solve, solve...

e = 2 r1 (r1 - r2) / (r2^2 - r1^2 - d^2)

a = r1 / (1-e)

How easy that was.


Case 1. Perihelion at Departure: Q1 = 0.

e = 2 r1 (r1 - r2) / ( r2^2 - r1^2 - d^2 )

If e is in [0,1) then an elliptical transfer orbit exists, and a = r1 / (1-e)

If e>1 then a hyperbolic transfer orbit exists, and a = r1 / (e-1)

If e<0 then there is no transfer orbit having perihelion at departure.


Case 2. Aphelion at Arrival: Q2 = pi radians.

e = 2 r2 (r1 - r2) / ( r1^2 - r2^2 - d^2 )

If e is in [0,1) then an elliptical transfer orbit exists, and a = r2 / (1+e)

If e is not in [0,1) then there is no transfer orbit having aphelion at arrival.


Case 3. Perihelion at Arrival: Q2 = 0.

e = 2 r2 (r2 - r1) / ( r1^2 - r2^2 - d^2 )

If e is in [0,1) then an elliptical transfer orbit exists, and a = r2 / (1-e)

If e>1 then a hyperbolic transfer orbit exists, and a = r2 / (e-1)

If e<0 then there is no transfer orbit having perihelion at arrival.


Case 4. Aphelion at Departure: Q1 = pi radians.

e = 2 r1 (r2 - r1) / ( r2^2 - r1^2 - d^2 )

If e is in [0,1) then an elliptical transfer orbit exists, and a = r1 / (1+e)

If e is not in [0,1) then no transfer orbit exists having aphelion at departure.


Since our example problem involves a transfer from Vesta to Earth, the trajectory will be inward bound and the transfer orbit we seek will be found in either Case #3 or Case #4.

Plugging in the numbers...

r1 = 2.179748
r2 = 1.005301

d = 2.028097

Case 1: e = -0.6519092 (unacceptable)
Case 2: e = -6.339075 (unacceptable)
Case 3: e = 6.339075 (hyperbolic orbit), a = 0.1882913
Case 4: e = 0.6519092 (elliptical orbit), a = 1.319533

So we shall propose that two orbits be investigated as possible transfer orbits between r1 at the time of departure and r2 at the time of arrival.

Proposed hyperbolic transfer orbit.

a = 0.1882913
e = 6.339075
i = 0.003996730 radians
L = ?
w = ?
T = ?

Proposed elliptical transfer orbit.

a = 1.319533
e = 0.6519092
i = 0.003996730 radians
L = ?
w = ?
T = ?

Jerry Abbott

[Jenab]

I define a subscript variable K to designate the apsidal endpoint of the trajectory. When K=1, the apside is at departure. When K=2, the apside is at arrival.

We must determine the sun-relative velocity of the spaceship at the apsidal end of its intended trajectory. Although it would be possible to use the general method (shown for ellipses in an earlier post), there is another way to proceed in this special circumstance.

We've found the heliocentric position vector to the apsidal endpoint of the intended trajectory, rK = [xK , yK, zK] and a unit vector normal to the transfer orbit Rn = [Xn , Yn, Zn].

We will obtain the cross product Rn x rK.

VxK'' = Yn zK - Zn yK

VyK'' = Zn xK - Xn zK

VzK'' = Xn yK - Yn xK

We find the magnitude of this double-primed vector.

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

We divide the components of the double-primed vector by the magnitude of the double-primed vector, obtaining a unit vector in the direction of the velocity in the transfer orbit at the apsidal endpoint.

VxK' = VxK'' / VK''

VyK' = VyK'' / VK''

VzK' = VzK'' / VK''

We refer to the Vis Viva equation to get the sun-relative speed in the transfer orbit at the apsidal endpoint. To keep the units consistent, you may want to convert [a] and [rK] from astronomical units to meters.

1 astronomical unit = 1.49597870691E+11 meters

GMsun = 1.32712440018E+20 m^3 sec^-2

Elliptical orbits: VK = [ GMsun ( 2/rK - 1/a ) ]^0.5

Hyperbolic orbits: VK = [ GMsun ( 2/rK + 1/a ) ]^0.5

VxK = VK VxK'

VyK = VK VyK'

VzK = VK VzK'

The vector VK = [VxK , VyK , VzK] is the velocity in the transfer orbit at the apsidal endpoint of the intended trajectory, referred to heliocentric ecliptic coordinates.

This method only works at apsides, not in general. It works at apsides because the heliocentric position vector and the sun-relative velocity are perpendicular to each other there.

Plugging in the numbers...

Xn = -4.498044E-4
Yn = -3.975890E-3
Zn = +0.9999920

The proposed hyperbolic transfer orbit.

x2 = +0.9989326
y2 = -0.1129745
z2 = 0

Vx2'' = +0.1129736
Vy2'' = +0.9989246
Vz2'' = +4.022925E-3

V2'' = 1.005301

Vx2' = +0.1123779
Vy2' = +0.9936575
Vz2' = +4.001713E-3

a = 0.1883722 AU
r2 = 1.005301 AU

V2 = 80463.3 m/s

Vx2 = +9042.3 m/s
Vy2 = +79953.0 m/s
Vz2 = +322.0 m/s

The proposed elliptical transfer orbit.

x1 = +0.5877971
y1 = -2.098983
z1 = -0.008080998

Vx1'' = +2.098999
Vy1'' = +0.5877888
Vz1'' = +3.281149E-3

V1'' = 2.179748

Vx1' = +0.9629546
Vy1' = +0.2696590
Vz1' = +1.505288E-3

a = 1.319533 AU
r1 = 2.179748 AU

V1 = 11902.5 m/s

Vx1 = +11461.5 m/s
Vy1 = +3209.6 m/s
Vz1 = +17.9 m/s

The significance of the work so far.

We have solved for three of the orbital elements of the transfer orbit (i, e, a), and we have determined a state vector for the transfer orbit at the apsidal endpoint of its trajectory.

State vector for the proposed hyperbolic transfer orbit.

x2 = +0.9989289
y2 = -0.1130118
z2 = 0
Vx2 = +9042.3 m/s
Vy2 = +79953.0 m/s
Vz2 = +322.0 m/s

State vector for the proposed elliptical transfer orbit.

x1 = +0.5877971
y1 = -2.098983
z1 = -0.008080998
Vx1 = +11461.5 m/s
Vy1 = +3209.6 m/s
Vz1 = +17.9 m/s

As yet, we know nothing about the velocity at the non-apsidal endpoint of the intended trajectory. But we shall soon fix that!

Jerry Abbott

[Jenab]

The angular momentum per unit mass, h, is equal to the cross product of the heliocentric radius vector and the velocity vector at the apsidal endpoint: rK x VK. Angular momentum is a conserved quantity, and the components of h are constant for the entire transfer orbit.

You will want to convert xK, yK, zK from astronomical units to meters.

hX = yK VzK - zK VyK

hY = zK VxK - xK VzK

hZ = xK VyK - yK VxK

h = [ (hX)^2 + (hY)^2 + (hZ)^2 ]^0.5

The longitude of the ascending node, L, of the transfer orbit can be calculated as follows:

L = arctan2(hX , -hY)

Plugging in the numbers...

The proposed hyperbolic orbit.

hX = -5.441887E+12 m^2 sec^-1

hY = -4.811776E+13 m^2 sec^-1

hZ = +1.210085E+16 m^2 sec^-1

h = 1.2100947E+16 m^2 sec^-1

L = 6.170569 radians

The proposed elliptical orbit.

hX = -1.745789+12 m^2 sec^-1

hY = -1.543129E+13 m^2 sec^-1

hZ = +3.881186E+15 m^2 sec^-1

h = 3.881217E+15 m^2 sec^-1

L = 6.170532 radians

The longitude of the ascending node had better be the same for both the proposed hyperbolic orbit and the proposed elliptical orbit! (That's because they share the heliocentric position vectors for arrival and departure and therefore occur in the same plane.)

Jerry Abbott

[Jenab]

cos(w+QK) = (xK cos L + yK sin L) / rK

if (i = 0) or (i = pi radians), then
sin(w+QK) = (yK cos L - xK sin L) / rK

if (i<>0) and (i<>pi radians), then
sin(w+QK) = zK / (rK sin i)

w = arctan2[ sin(w+QK) , cos(w+QK) ] - QK

If necessary, adjust w to the interval [0, 2 pi).

Plugging in the numbers...

The proposed hyperbolic transfer orbit.

x2 = +0.9989326
y2 = -0.1129745
z2 = 0

r2 = 1.005301 AU

i = 0.003996730 radians
L = 6.170569 radians

Q2 = 0

(This orbit's perihelion occurs at arrival at a point on Earth's orbit.)

w = 0

Yes, the argument of the perihelion of the proposed hyperbolic transfer orbit is zero. Why? Because the arrival occurs at a point which is both perihelion and on the ecliptic plane, where the Earth is to be found. This node is the ascending node because Vz2 is positive. Hence, the argument of the perihelion, for this orbit, is zero.

The proposed elliptical transfer orbit.

x1 = +0.5877971
y1 = -2.098983
z1 = -0.008080998

r1 = 2.179748 AU

i = 0.003996730 radians
L = 6.170532 radians

Q1 = pi radians

(This orbit's aphelion is at departure from Vesta.)

w = 1.9560970 radians

Jerry Abbott

[Jenab]

I define a subscript variable J to designate the non-apsidal endpoint of the intended trajectory. If the apside is at departure (K=1), then J=2. If the apside is at arrival (K=2), then J=1.

We want to find the heliocentric vector to the non-apsidal endpoint in the coordinate system in which the XY plane contains the transfer orbit, with the +X axis extended through the transfer orbit's perihelion. Once we have this canonical vector, we can use it to get the true anomaly at the non-apsidal endpoint of the intended trajectory.

xJ’ = yJ sin L + xJ cos L

yJ’ = yJ cos L - xJ sin L

zJ’ = zJ

xJ’’ = xJ’

yJ’’ = zJ’ sin i + yJ’ cos i

zJ’’ = zJ’ cos i - yJ’ sin i

xJ’’’ = yJ’’ sin w + xJ’’ cos w

yJ’’’ = yJ’’ cos w - xJ’’ sin w

zJ’’’ = zJ’’

This triple-primed position vector is the one we were looking for. Important check: Within a reasonable allowance for roundoff error, the value of zJ’’’ should be zero.

QJ = arctan2 (yJ’’’ , xJ’’’)

Finding the true anomaly of the non-apsidal endpoint of the intended trajectory is a milestone in solving the transfer orbit problem.

Plugging in the numbers...

The proposed hyperbolic transfer orbit.

x1 = +0.5878052
y1 = -2.098982
z1 = -8.082041E-3

L = 6.170569 radians
i = 3.99673E-3 radians
w = 0

...intermediate steps skipped...

x1''' = +0.8199621
y1''' = -2.019646
z1''' = -1.008436E-5

Q1 = 5.098051 radians

The proposed elliptical transfer orbit.

x2 = +0.9989289
y2 = -0.1130118
z2 = 0

L = 6.170532 radians
i = 3.99673E-3 radians
w = 1.9560970 radians

...intermediate steps skipped...

x2''' = -0.3778306
y2''' = -0.9315980
z2''' = +8.933364E-11

Q2 = 4.327089 radians

Jerry Abbott

[Jenab]

For hyperbolic orbits, we proceed as follows:

B = cosh uJ = (1/e) (1 + rJ /a)

uJ’ = ln [ B + (B^2 - 1)^0.5 ] = arc-cosh(B)

(Depending on whether your calculator has an arc-cosh button on it.)

If sin QJ < 0 then uJ = -uJ’ else uJ = uJ’

MJ = e sinh uJ - uJ


For elliptical orbits, we proceed as follows:

sin uJ = (rJ /a) sin QJ / (1-e^2)^0.5

cos uJ = e + (rJ /a) cos QJ

uJ = arctan2(sin uJ , cos uJ)

MJ = uJ - e sin uJ


The proposed hyperbolic transfer orbit.

e = 6.33678
a = 0.1883722 AU
r1 = 2.179749 AU
Q1 = 5.098051 radians

u1' = 1.307609 radians
sin Q1 is negative.
u1 = -1.307609 radians

M1 = -9.550008 radians

The proposed elliptical transfer orbit.

e = 0.6519092
a = 1.319533
r2 = 1.005301
Q2 = 4.327089 radians

sin u2 = -0.9310417
cos u2 = +0.3655728

u2 = 5.086543 radians
M2 = 5.693351 radians

Jerry Abbott

[Jenab]

The proposed hyperbolic transfer orbit.

GMsun = 1.32712440018E+20 m^3 sec^-2

1 astronomical unit = 1.49597870691E+11 meters

You will want to convert the semimajor axis [a] from astronomical units to meters.

m = (86400 sec/day) ( GMsun / a^3 )^0.5

Where m is the hyperbolic mean motion in radians per day.

Perihelion at Departure: K=1, J=2.
dt = M2 / m, if Q1=0 (thus M1=0)

Perihelion at Arrival: K=2, J=1.
dt = -M1 / m, if Q2=0 (thus M2=0)

In our example, we were investigating a hyperbolic orbit connecting the positions of Vesta at the time of departure and of Earth at the time of arrival. Since the orbit is inward bound (arrival is closer to the sun than departure), the perihelion was made to coincide with arrival. Thus Q2=0.

M1 = -9.550008 radians
m = 0.2104051 radians per day
dt = -M1/m = 45.4 days

The proposed elliptical transfer orbit.

P = (365.256898326 days) a^1.5

m = 2 pi radians / P

Where P is the orbital period in days, and m is the mean motion in radians per day.

Apside at Departure: K=1, J=2.
dt = M2 / m, if Q1=0 (thus M1=0)
dt = (M2-pi) / m, if Q1=pi (thus M1=pi)

Apside at Arrival: K=2, J=1.
dt = (2 pi - M1) / m, if Q2=0 (thus M2=0)
dt = (pi - M1) / m , if Q2=pi (thus M2=pi)

In our example, we were investigating an elliptical orbit connecting the positions of Vesta at the time of departure and of Earth at the time of arrival. The orbit had its aphelion at the point of departure. Thus Q1=pi radians.

P = 553.6415 days
m = 0.01134884 radians per day
M2 = 5.693351 radians
dt = (M2-pi) / m = 224.85 days

Significance of our work so far.

In an earlier post, I wrote:

$this->bbcode_second_pass_quote('Jenab', 'W')e assert that a spaceship will depart from Vesta at

t1 = 4 February 2004, 19h 12m UT
t1 = JD 2453040.3

We assert that the spaceship will arrive at Earth at

t2 = 16 September 2004, 21h 36m UT
t2 = JD 2453265.4

We thus have established a required transit time for the spaceship, going from Vesta to Earth, of t2-t1, or 225.1 days.

[Jenab]

The canonical velocity in a hyperbolic orbit, in general, is found from

Vx’’’ = -(a/r) { GMsun / a }^0.5 sinh u

Vy’’’ = +(a/r) { GMsun / a }^0.5 (e^2 - 1)^0.5 cosh u

Vz’’’ = 0

The canonical velocity in an elliptical orbit, in general, is found from

Vx’’’ = -sin Q { GMsun / [ a (1-e^2) ] }^0.5

Vy’’’ = (e + cos Q) { GMsun / [ a (1-e^2) ] }^0.5

Vz’’’ = 0

These triple-primed vectors would be rotated (negatively) by the angular elements of the orbit (w,i,L) to heliocentric ecliptic coordinates.

Rotation by the argument of the perihelion, w.

Vx'' = Vx''' cos w - Vy''' sin w

Vy'' = Vx''' sin w + Vy''' cos w

Vz'' = Vz''' = 0

Rotation by the inclination, i.

Vx' = Vx''

Vy' = Vy'' cos i

Vz' = Vy'' sin i

Rotation by the longitude of ascending node, L.

Vx = Vx' cos L - Vy' sin L

Vy = Vx' sin L + Vy' cos L

Vz = Vz'

By comparing transit times for Earth and for the spaceship from their respective starting points (at the time of departure) to the arrival location, we determined that the elliptical orbit in our example problem was a valid (or almost so) transfer orbit from Vesta to Earth. That orbit had aphelion at the departure from Vesta.

Using a special method (involving a cross product and the Vis Viva equation) that works only at apsides, we determined the HEC velocity in the transfer orbit at departure: [Vx1, Vy1, Vz1].

Using the general method for finding velocity in an elliptical orbit, we obtained Vesta's own HEC velocity at departure:

[ Vx(preburn) , Vy(preburn) , Vz(preburn) ]

The delta-vee for departure is

dV1x = Vx1 - Vx(preburn)
dV1y = Vy1 - Vy(preburn)
dV1z = Vz1 - Vz(preburn)

We may now use the general method for finding velocity in an elliptical orbit to determine the velocity in the transfer orbit at the non-apsidal endpoint of the intended trajectory (i.e., at Earth). The result will be: [Vx2, Vy2, Vz2].

We may use the general method, also, to determine the HEC velocity of Earth in its own orbit at the moment of arrival, namely:

[ Vx(rendezvous) , Vy(rendezvous) , Vz(rendezvous) ]

The delta-vee required of the spaceship to match speed with Earth is

dV2x = Vx(rendezvous) - Vx2
dV2y = Vy(rendezvous) - Vy2
dV2z = Vz(rendezvous) - Vz2

Plugging in the numbers...

Transfer orbit from Vesta (JD 2453040.3) to Earth (JD 2453265.4).

Ellipse with aphelion at departure.

a = 1.319533 AU
e = 0.6519092
i = 0.2289958 degrees
L = 353.5454 degrees
w = 112.0761 degrees

True anomaly of arrival: 247.9239 degrees.

Delta-vee at departure (HEC)
dV1x = -8.772 km/sec
dV1y = -1.514 km/sec
dV1z = +2.619 km/sec

Delta-vee at arrival (HEC)
dV2x = +20.487 km/sec
dV2y = +1.515 km/sec
dV2z = -0.103 km/sec

The above delta-vees are in ecliptic coordinates. A coordinate rotation will get the velocity vector into celestial coordinates, which will permit the spaceship pilot to orient his thrust by observing the starfield into which he must accelerate, in comparison with a conventional star atlas, in order to enter the transfer orbit.

Jerry Abbott

[Jenab]

If the spaceship pilot has a star atlas referred to ecliptic coordinates, he won't need to do this step. But since most star atlases use celestial coordinates, I thought it best to include the rotation from ecliptic to celestial.

dVx' = dVx

dVy' = dVy cos q - dVz sin q

dVz' = dVy sin q + dVz cos q

Where [q] is Earth's obliquity at the moment of thrust. The J2000 value of the obliquity is

q = 23.439281 degrees = 0.40909263 radians

The magnitude of the delta-vee is independent of the rotation, of course,

dV' = dV = { dVx^2 + dVy^2 + dVz^2 }^0.5

The right ascension of the delta-vee (hence also the thrust) is

dVRA = arctan2( dVy' , dVx' )

The declination of the delta-vee (hence also the thrust) is

dVDEC = arcsin( dVz' / dV' )

Putting in the numbers...

Departure.

dV1x = -8.772 km/sec
dV1y = -1.514 km/sec
dV1z = +2.619 km/sec

dV1 = 9.279 km/sec
dVRA1 = 13h 1m 57.4s
dVDEC1 = +11d 11m 8s

The departure thrust will be roughly toward Vindemiatrix, Virgo.

Arrival.

dV2x = +20.487 km/sec
dV2y = +1.515 km/sec
dV2z = -0.103 km/sec

dV2 = 20.544 km/sec
dVRA2 = 0h 4m 11.0s
dVDEC2 = 0d 1m 29s

The arrival thrust will be toward a point in Pisces.

Jerry Abbott

[Jenab]

Once you have the above procedure coded, you can use trial-and-error to find valid transfer orbits.

You'd choose a time for departure (t1) and a time for arrival (t2), thereby fixing the positions of the rendezvous planet at those times, as well as fixing the transit time for the rendezvous planet (t2-t1).

Then you crank through the procedure to determine whether (or not) the proposed transfer orbit connecting the location of the spaceship at t1 with the rendesvous planet's location at t2 actually does occur in the same time difference.

If the transit time along the proposed transfer orbit isn't even close, discard that proposed transfer orbit: it's no good.

If the transit time along the proposed transfer orbit is close, you'd play around with different departure or arrival times until you have an exact transit time match with both the rendezvous planet and the spaceship.

Now, there's another way to proceed, if you don't mind the work. For any two objects - Vesta and Earth in our example - it would be possible to construct a "Big Book of Object 1 to Object 2 Transfer Orbits" comprised mostly of a table indexed by the heliocentric longitudes of the two objects. The index interval might be tenths of a degree, or milliradians, or something convenient to spaceship pilots.

This table is the result of some reverse interpolation. To construct it to begin with, you need to use the mean anomalies of Object 1 and Object 2 as indices going forward, with the transfer orbit(s) elements, transit times, and the heliocentric longitudes of both planets at the time of departure being outputs.

Finding the heliocentric longitude of the preburn position is relatively simple: one use of the arctan2 function on y1,x1.

The rendezvous planet must be backtracked in its own orbit by the transit time in the transfer orbit. To do this, a reduction of the rendezvous planet's orbital elements is conducted with the "time of interest" being t2-dt. The resulting HEC position vector might be labeled [xbt, ybt, zbt]. The heliocentric longitude of the rendezvous planet at the time of departure is arctan2(ybt,xbt).

The output is only made in the event that t2-dt = t1, to within some appropriately chosen error tolerance. Otherwise, the proposed transfer orbit in that case would be found to be invalid and accordingly of no interest.

Having found the heliocentric longitudes of both planets at the moment of departure, one now has the indices needed to construct the reverse interpolation table from which the sampled set of valid transfer orbits may be read off as a first rough approximation for spaceship pilots wanting to go from a point in one orbit to rendezvous with a planet in a different orbit.

Jerry Abbott

[Jenab]

The elliptical transfer orbits treated in the preceding posts are what I've come to call "elliptical transfer orbits of the short path," meaning that the transit object moves through an arc of true anomaly less than pi radians.

For each elliptical transfer orbit of the short path, there exists an elliptical orbit of the long path, in which the transit object moves along the same orbit in the opposite direction.

Long path elliptical orbits differ from short path elliptical orbits in the following ways.

At each point where the HEC position vectors are identical, the HEC velocity vectors differ (only) in sign.

The sum of the transit times for the long and short paths is equal to the period of the transfer orbit.

The inclinations sum to pi radians.

The angular momentum vector of the long path orbit is the reverse of that of the short path.

The longitude of the ascending node on the long path will be that of the short path plus or minus pi radians, whichever will keep it inside the interval [0 , 2 pi).

The sum of the short-path and long-path arguments of the perihelion should add up to pi radians +/- some whole multiple of 2 pi radians.

The sum of the non-apsidal endpoint true anomalies for the short path elliptical orbit and for the long-path elliptical orbit is 2 pi.

And the like.

There is yet another class of intercept orbit, which might be called a super-periodic transfer orbit. They may be either prograde or retrograde, as we normally reckon such, and involve transit times exceeding the transfer orbit's period.

Such extra possible orbits gives a spaceship pilot more opportunities to find a valid transfer orbit, however at the cost of increasing his transit time, often greatly.

Jerry Abbott

[Jenab]

In the example used for calculating a transfer orbit, we had a spaceship departing from Vesta at a certain time. We used Vesta's orbital elements and the departure time to calculate the preburn state vector for Vesta (hence also for the rocket) immediately before the rocket fires its engines to enter the transfer orbit.

Vesta is in an elliptical orbit, and the method shown for obtaining the preburn state vector was the method appropriate for elliptical orbits.

But suppose the rocket had been in a hyperbolic orbit, relative to the sun, instead? The calculation must proceed somewhat differently, in that case.

There is, of course, no period associated with a hyperbolic orbit. But we can determine an equivalent to the mean motion:

m = ( GMsun / a^3 )^0.5

Where

GMsun = 1.32712440018E+20 m^3 sec^-2

and remember to enter the semimajor axis of the hyperbolic orbit in meters. One astronomical unit equals 1.49597870691E+11 meters.

Likewise, the mean anomaly has no special geometric meaning for a hyperbolic orbit, but it nonetheless remains mathematically convenient as an intermediate quantity. The mean anomaly is zero at perihelion, negative prior to perihelion, and positive after perihelion.

M = m (t - T)

where [t] is the moment of interest (e.g. the time of departure) and [T] is the time of perihelion passage. This difference of time is entered in seconds, and M will result in radians. Writing the equation fully:

M = {GMsun / a^3)^0.5 (t - T)

If you'd rather input astronomical units for [a] and days for (t-T), then

M = 0.01720209895 (t-T) a^-1.5 AU^1.5 day^-1

and, again, M will be in radians.

It is important to remember that we do not correct the mean anomaly of hyperbolic orbits to the interval [0,2 pi). If it comes out negative, leave it that way.

Kepler's equation for hyperbolic orbits is

M = e sinh u - u

Where (u) is the hyperbolic eccentric anomaly, which, along with (M), must be in radians. As it was in the elliptical case, the equation is transcendental in the variable that we are trying to find, and we must use a differential calculus method for solving it.

Danby's Method for finding the eccentric anomalies of hyperbolic orbits.

u(0) = 0

Repeat over index j

f0 = e sinh u(j) - u(j) - M
f1 = e cosh u(j) - 1
f2 = e sinh u(j)
f3 = e cosh u(j)
d1 = -f0 / f1
d2 = -f0 / [ f1 + (d1)(f2)/2 ]
d3 = -f0 / [ f1 + (d1)(f2)/2 + d2^2 (f3)/6 ]
u(j+1) = u(j) + d3

Until |u(j+1)-u(j)| < 1E-12

The converged value for (u) from this loop is the eccentric anomaly for the hyperbolic orbit. We don't correct (u) to the interval [0,2 pi) either; if it comes out negative, we leave it that way.

Finding the canonical position vector.

The true anomaly is found from

Q' = arccos { (e - cosh u) / (e cosh u - 1) }

if u>0 then Q = Q'
if u=0 then Q = 0
if u<0 then Q = 2 pi - Q'

The heliocentric distance is

r = a (e cosh u - 1)

The canonical position vector is

x''' = r cos Q
y''' = r sin Q
z''' = 0

The canonical velocity vector is

Vx''' = -(a/r) { GMsun / a }^0.5 sinh u

Vy''' = +(a/r) { GMsun / a }^0.5 (e^2 - 1)^0.5 cosh u

Vz''' = 0

The triple-primed position and velocity, although relative to the sun, are not yet presented in the heliocentric ecliptic coordinate system. They are each rotated (negatively) by the orbital elements w (about the z''' axis), i (about the x'' axis), and L (about the z' axis) in order to appear in the (unprimed) HEC system.

A check on the magnitude of the velocity (i.e., the sun-relative speed) is available:

Vx^2 + Vy^2 + Vz^2 = GMsun { 2 / (x^2 + y^2 + z^2)^0.5 + 1/a }

The state vector of an object in a hyperbolic orbit of elements [ a , e , i , L , w , T ] at the moment of interest [t] is

[ x , y , z , Vx , Vy , Vz ]

Jerry Abbott

[Jenab]

We have now reached the end of this dissertation. I will entertain questions and engage in a certain amount of pedagogical banter.

Jerry Abbott

[Geingreen]

Dude, you forgot the Di-Lithium Crystals for your reactor core.

I also agree Afghani hash is the best, but have you tried some of our more ethnogenic "silly" psilocybin vessels?

http://www.erowid.org/plants/mushrooms/mushrooms.shtml

I think a mushroom cap would make a great vesel design, and the stem a great booster rocket (solid fuel) for when the warp drive is off line.
))
==)))
))

woah.... kinda looks like a date palm, but ya get the idea...

~ ~
-
<>

Hey look a self portrait, hope it dosent cause the server to crash with the out of this world grapix.

[jato]

I fly with Orbiter, a space simulator. You can get a free copy here:

http://www.medphys.ucl.ac.uk/~martins/orbit/orbit.html




There are fictional spacecraft as well. I thought it was cool when I performed my first transfer orbit from the Earth to the Moon.

I wanted to be an astronaut since I was a little guy.

Too bad space travel and exploration is coming to an end.

[Giengreen]

Ok dudes and dudets, that stinks , my way cool grapixs got royaly screwed.


=0 there a simplified shroom-ship.

I hope that clears up the techinical difficuties.



-------------------
Ghienges Khan has nothing on me, cause in my mind I rape rob pillaige everything I see...