Transfer orbits for dummies A hillbilly tutorial.

[Jenab]

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]

Restrictions inherent in this procedure.

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).

Jerry Abbott

 

[Jenab]

Initial data.

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]

Preburn and Rendezvous State Vectors.

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]

Inclination of the transfer orbit.

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]

Eccentricity and semimajor axis of the transfer orbit.

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]

Velocity in transfer orbit at apsidal end of intended trajectory.

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]

Longitude of the ascending mode of the transfer orbit.

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]

Argument of the perihelion of the transfer orbit.

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]

True anomaly in the transfer orbit of the non-apsidal end of the intended trajectory.

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]

Eccentric and mean anomalies in the transfer orbit at the non-apsidal endpoint.

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 transit time in the transfer orbit.

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:

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.

Obviously the proposed hyperbolic orbit fails the transit time test and is therefore not a valid transfer orbit. A spaceship leaving Vesta at the time of departure will reach the point of arrival long before Earth does, which would be rather foolish. The spaceship pilot would lose his job.

However, the proposed elliptical orbit is a near-match (six hours difference). Some adjustment in the time of arrival will probably result in a perfect match.

The only thing left to do is calculate the delta-vees.

Jerry Abbott

 

[Jenab]

The delta-vees for departure and for arrival.

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

Since the spaceship's transit time wasn't exactly a match for the Earth's, I won't plug in the numbers on this occasion. These equations do, however, show the way it is done.

UPDATE: See a later post for the delta-vees for departure and arrival for the elliptical transfer orbit.

Jerry Abbott

 

[Jenab]

Long-path and superperiodic elliptical transfer orbits.

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]

Comments? Applause? Honorary master's degree?

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

Jerry Abbott

 

[Jenab]

Matching up transit times.

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

 

[Janitor]

Call me impressed!

One topic that came up here some time ago was whether the advantage to burning as much of a rocket's propellant while low in a gravitational field applies only when the field is considered as varying in strength (as when a probe flies by a planet), or if it also applies to modest rockets going straight up to an apogee and then falling straight back down in a field that is considered constant, e.g. a flight of a couple of miles up from Earth's surface. A key assumption was made that the beefing-up of the structure to allow for higher acceleration would be ignored. So you've got two rockets to compare, both with same structure weight and same propellant weight and same total impulse. Your thoughts?

 

[Jenab]

It sounds like the first rocket is taking advantage of the gravity well of a planet by holding on to its fuel's mass until most of its potential energy has become kinetic energy, then blowing the fuel away contrary to the spaceship's velocity, which transfers some (how much depends on the exhaust speed) of the fuel mass's kinetic energy into the spaceship. The ship then rises out of the gravity well still having all of the energy from the fall, the energy from the thrust, plus the gravitational potential energy stolen from the fuel.

For the second spaceship, there's also an advantage to losing the weight of the propellant as fast as possible. The least efficient (zero efficiency) launching rocket is one that can do no more than hover by making thrust equal only to its own weight. When launching from a planet with an atmosphere, though, there may be a tradeoff between thrust and aerodynamic pressure, so that you'd want to get above the soupy part of the atmosphere before flooring the pedal.

I'm really not much of an expert in astrodynamics. I really am a hillbilly. I live in the Alleghenies near Hillsboro, West Virginia, where I raise goats and grow potatoes and apple and walnut trees. Now and then I like to figure things out, like I did with transfer orbits a few years ago.

Jerry Abbott

 

[remcook]

Depends on what you want I guess.

If you want to gain speed using planetary a flyby, it's more advantageous to fire at pericentre. This is because speed is largest there. The speed that you add by burning fuel will always be the same (it depends on the fuel, engine, etc). So, because kinetic energy goes with the square of the velocity, adding that same velocity to a large speed will give you a larger increase in kinetic energy (and because potential energy stays constant during the burn, the total energy increase is the same as the kinetic energy increase for both cases.)

 

[Janitor]

Thanks to both of you for your responses.

For the case of the modest rocket going straight up and straight down in a constant-g field... If we ignore the effect of an atmosphere, we ought to be able to graph apogee vs. thrust and get some sort of curve, with fixed parameters being structure mass, propellant mass, and total impulse, and the constraint that thrust times burn time equals the total impulse. As Jenab noted, "The least efficient (zero efficiency) launching rocket is one that can do no more than hover by making thrust equal only to its own weight." To my intuition, that means the curve will not have its peak over toward the left side of the graph where the thrust is low. In fact, my intuition suggests the curve is monotonically increasing with increasing thrust. (In the real world, structure mass would have to increase when thrust gets bigger to keep the rocket from collapsing or blowing up, so that its non-constancy would have to be taken into account, so the real-world version of the graph would actually have some maximum apogee for some finite value of thrust, beyond which increasing thrust lowers the apogee; "the point of diminishing returns" you might call it.)

Remcook, when you say, "because potential energy stays constant during the burn," are you making the approximation that the burn time is short compared to the time it takes for the gravitational potential to change appreciably when the probe is near its minimum point in the potential field?

And another question for Jenab: can a West Virginian still grab onto a slow-moving coal hopper and go places in your state? Or are those days long gone?

As it happens, I recently read an article on M&K Junction in West Virginia, across the Cheat River from Rowlesburg.

 

[remcook]

Remcook, when you say, "because potential energy stays constant during the burn," are you making the approximation that the burn time is short compared to the time it takes for the gravitational potential to change appreciably when the probe is near its minimum point in the potential field?

Yes

Both in apogee and in perigee the velocity is perpendicular to the gravity acceleration anyway (no gravity loss). But the satellite moves quicker at perigee of course. But if the burn time is small...

 

[Jenab]

And another question for Jenab: can a West Virginian still grab onto a slow-moving coal hopper and go places in your state? Or are those days long gone?

As it happens, I recently read an article on M&K Junction in West Virginia, across the Cheat River from Rowlesburg.

I never tried jumping on a coal hopper. I'm in Pocahontas County, about 4 miles NW of Hillsboro, or 2 miles NW of Mill Point, or 8 miles SW of Marlinton, or 33 miles north of Lewisburg. The nearest river is the Greenbriar, which runs N-S a little to the east of Hwy 219. I moved up here from Alabama in 1998.

Jerry Abbott

 

[Janitor]

Thanks for confirming, Remcook.

Jenab, M&K Junction is three counties north of you in Preston. Ever been to Harper's Ferry or Spruce Knob? I have not, but those are places on my to-do list, to get to one of these years.

 

[Jenab]

I don't travel much, except locally to shop for groceries, animal feed, to buy building supplies/tools, or to send and receive mail. I came originally to West Virginia to work for Dr. William Pierce as an editor for National Vanguard Books. When he died, the National Alliance came under new management, and I soon discovered that I had irreconcilable differences with the new chief, so I quit.

Now I live very pastorally on my ranch/orchard, and circumstances don't give me a lot of time for sightseeing. Once I chased a presumptuous young black bear away from my trash can, off my porch and into the woods. Another time, I shot a bobcat as he was stalking my chicken, Mr. Bronfman; but unfortunately the same chicken met an untimely demise two weeks later under the talons of a hawk.

My previous cat, Father Wiggly, had a run-in with the bobcat a week before I had my own show-down with it. ("Bobcat, this ranch isn't big enough for the two of us.") Wiggles managed to get away, but he was scarred up some. But he disappeared later that year and never showed up again. So now I have two new kittens, Spooky and Goblin.

Most of my trouble with wildlife, though, comes from deer entering the property to nibble on my young apple trees. I have each tree caged inside a circle of fencing 3 feet wide and six feet high. That works, most of the time, but I lost three of my 30 trees last winter to rabbits nibbling off the bark at the bottom. So now I've added another circumferential layer of chicken wire down there.

Jerry Abbott

 

[selfAdjoint]

Sounds like a beautiful life. Are you anywhere near Spencer? That's where my father grew up. He used to spend summers on his uncle's farm and had a bunch of stories just like yours that he used to tell me.

 

[Jenab]

Sounds like a beautiful life. Are you anywhere near Spencer? That's where my father grew up. He used to spend summers on his uncle's farm and had a bunch of stories just like yours that he used to tell me.

It's a good place to live, especially if you like being close to nature. I like to sit on my front porch and read. My house is on a rise, with a road coming up from a valley to the south and halfway circling my house on a spiral upward. The hilly nature of the area is such that my front porch faces an even higher rise, beyond the road, to the NE. Yet higher (and more distant) hilltops are visible to the north and west. And there's woods covering most of the hills.

On a windy day, the leaves make a very nice susurration, which is a kind of leafy "white noise," during which the trees seem to oscillate between dark green (when you see the tops of the leaves) and very light greenish white (when the wind shows you the lighter flip side of the leaves). And on top of that, there's cloud-play alternately shadowing the foreground, then the background, etc.

But the best thing about this area is that it is smack in the middle of a demographic safe-zone. When fossil fuels become, in 20 years or so, too expensive for government subsidies to keep mechanized agriculture going, there's going to be famine along the coast and in urban areas generally, which means that the same people who will rob you now for the money in your wallet will be breaking into homes looking for cans of beans or sacks of rice.

These "safe-zones" I refer to are areas with low population density, especially in regard to demographic groups that have an elevated statistical propensity for causing crime. Pardon the circumlocution. We'll have a food shortage here, but we have a chance at maintaining ourselves with our own efforts, and there will be fewer bandits to contend with.

Spencer is west of me. It's in western West Virginia, over by Ohio. I'm in eastern West Virginia, over by Virginia.

Jerry Abbott

 

[Janitor]

Our federal government in the form of the I.R.S. and welfare agencies has a bad track record of rewarding people for having lots of kids. My opinion is that coming shortages of fuel plus coming shortages (out West anyway) of water are aggravated by population growth. I doubt whether more than a tiny handful of congresspersons would agree with me on that, though.

 

[Jenab]

I'm not a biological scientist, but I have a halfway decent understanding for how evolution works to achieve favorable adaptations and how nature's usual methodology can be frustrated.

Births per se aren't the problem. It's the lack of natural rigor to life that leads to an early death for the biologically unfit, the poorly adapted, the defective. Fossil fuels enabled mankind to remove this rigor from his existence for a time, but at the price of accumulating a load of bad genes that increasingly require mechanical aid to compensate.

The depletion of fossil fuels will bring the bill for these accumulated genetic costs due for payment. No longer will mankind be able to evade natural rigor, which is essentially as rigorous as it ever was, but we are much less adapted than we once were for meeting its challenges.

In fact, the demographic groups that maintain a high birthrate, despite the disaster that will sweep the world, will have an advantage over those that attempt to redress the food shortage by reducing their family size. Conflict over resources is inevitable, and in the usual course of things people will tend to sort themselves by biological similiarity first and foremost. Other things being equal, the side with the biggest army wins.

It is time for us to consider how we will maintain our existence, and that of those who are important to us, in the tempestuous struggles that lie ahead. I hadn't thought to introduce politics and biology into a discussion of celestial mechanics, but things that have been on my mind have a way of coming up.

For those of you with the money and inclination to change your residences, here is a short and partial list of safe-zones in the United States:

Virginia/West Virginia.
VA: Highland, Bath.
WV: Pocahontas, Pendleton, Tucker, Randolph, Grant.

Georgia/North Carolina.
NC: Clay, Cherokee.
GA: Fannin, Union, Towns.

Arkansas/Missouri.
AR: Fulton, Baxter, Marion, Izard.
MO: Howell, Oregon, Ozark, Shannon, Douglas, Wright, Texas.

Wisconsin/Iowa.
WI: Grant, Juneau, Vernon, Richland, Crawford.
IA: Winnishiek, Allamakee, Howard, Clayton.

Nebraska/Kansas.
KS: Cheyenne, Rawlins, Decatur, Norton, Phillips, Smith, Jewell, Republic, Cloud, Mitchell, Osborne, Rooks, Grayham, Sheridan, Thomas, Sherman.
NE: Dundy, Hitchcock, Red Willow, Furnas, Harlan, Chase, Hayes.

Central Nebraska.
NE: Grant, Hooker, Thomas, Blaine, Loup, Garfield, Wheeler, Holt, Custer, Valley, Greeley, Logan, McPherson, Arthur, Brown, Rock.

No doubt there are others that I have not found.

You can find other relatively safe areas by going to the Census Bureau's "American Factfinder" pages, which begin on

http://factfinder.census.gov/home/s...n.html?_lang=en

From there, click the link "Data Sets" under the header "Getting Detailed Data." It's in the line: "Expert User? Go directly to Data Sets."

From there, click the link "List all maps" under the header "Select from the following options."

From there, choose from the list a demographic statistic that you want to display on the map. You'll find this statistic especially important:

* Persons per Square Mile: 2000

If you are interested in racial data, the Census Bureau has that as well, and maps to the same scale and with the same center can be overlaid to provide whatever cumulative measure of threat-assessment you consider to be relevant. You might also try convolving the color-coded county-level output with a Gaussian probability distribution to account for some degree of threat from mobile predators.

Jerry Abbott

 

[Jenab]

It looks like I allowed roundoff error to accumulate as I lazily truncated decimals to save button pushes on the calculator. I just now programmed the transfer orbit procedure and got a transit time of 224.85 days for the spaceship, as compared with 225.1 days for the Earth. So my elliptical transfer orbit was better than I thought.

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

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

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 order to enter the transfer orbit.

Jerry Abbott

 

[Jenab]

Rotation of sun-relative delta-vee into local celestial coordinates.

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

 

[remcook]

Jenab: if I may ask: why did you post this stuff? I mean..how did you come up with this? Why would you want to calculate this?
Ehm..can't find the right question to ask.
Hope I am not rude. you're doing a great job

just wondering...

 

[Nereid]

*SNIP

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

The departure thrust will be roughly toward [arrgh! I don't have a star atlas!]

The arrival thrust will be roughly toward [arrgh! I don't have a star atlas!]

I can't vouch for (or against) any of these, but some are free, and all can work on your PC (so they claim):
http://www.seds.org/billa/astrosoftware.html

 

[Jenab]

Jenab: if I may ask: why did you post this stuff? I mean..how did you come up with this? Why would you want to calculate this?
Ehm..can't find the right question to ask.
Hope I am not rude. you're doing a great job
just wondering...

There's really two questions here. First, what provoked a goatherd into learning celestial mechanics? Second, what prompted me to post a demonstration on calculating transfer orbits here?

It would be poetic for me to point out that herdsmen have been watching the sky for thousands of years and must somehow acquire a natural curiosity about how things work up there. But really what got me started on celestial mechanics were two books by Robert A. Heinlein: The Moon is a Harsh Mistress and The Rolling Stones. In both books, celestial mechanics plays a role in the story.

I decided that I'd see how tough it was to learn how to fly a spaceship by the seat of my pants, even if our ridiculous Government would never let me get close to one. It turns out that basic astronavigation isn't as hard as it's usually cracked up to be. The concepts are fairly simple. It's the presentation that's usually wrong: too bent on throwing the history and derivation behind every physical law discovered since Ptolemy into unprepared faces.

The way to teach celestial mechanics is the way I've done it here. Put all the important equations out in the open. Minimize on the derivations - those can come later as advanced, further reading topics. Keep the notation simple, at the high school algebra level as much as possible. Follow the essential equations with a fully worked-out example. Dish the work out in bite-sized chunks.

Voila, now anybody can follow the procedure. Or almost.

I can't vouch for (or against) any of these, but some are free, and all can work on your PC (so they claim):
http://www.seds.org/billa/astrosoftware.html

Thanks! I downloaded one of the free ones. And I bought a 20th edition Norton's Star Atlas from Amazon.com.

Jerry Abbott

 

[Jenab]

Determining a state vector in a hyperbolic orbit

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]

Transfer orbit program available, free download

My first program was in GWbasic (old DOS interperter programming software), which some of you might not have. So I rewrote it in FREE PASCAL and compiled it. You can download (34 kB) a zip file containing the source code (transit.pas), a compiled executable (transit.exe), and the necessary external data file (transit.dat) from

http://www.jabpage.org/features/transit.zip

NOTICE: I've caught a couple of control-flow bugs myself since uploading this program originally. And I fixed a failure to restore the calculated argument of perihelion to the interval [0,2 pi). So you might want to download the latest version.

The data file initially will contain the orbital elements of Vesta and of Earth and the departure and arrival dates that I used in this thread.

Free download. No restrictions on copying, modifying or sharing, except don't charge the fellow you share it with anything: it must stay freeware.

Bug alerts and suggestions for improvement can be posted here or emailed to goatlyones@moonshinehollow.com.

Jerry Abbott