CELESTIAL MECHANICS CHALLENGE: Find the Transfer Orbit

[Jenab]

Okay fellows. Here's my challenge.

There existed a transfer orbit from asteroid Vesta (IAU number 4) to Earth with a departure date of 4 February 2004.

(Go look up the orbital elements of Earth and Vesta; that's part of the work!)

Find: The classical orbital elements of the transfer orbit, the transit time, and the change-of-velocity vectors at departure and at arrival.

The classical orbital elements are: semimajor axis, eccentricity, inclination to the ecliptic, longitude of the ascending node, argument of perihelion, and Julian date of perihelion passage.

The delta-vee vectors should be referred to heliocentric ecliptic coordinates.

Please show your enough of your work to prove that you know what you're doing.

High praise will reward the poster who provides the correct answers first! Good luck to all who try.

Jerry Abbott

[Jenab]

Departure from Vesta:
4 February 2004, 19h 12m UT
JD 2453040.3

Arrival at Earth.
16 September 2004, 21h 36m UT
JD 2453265.4

The transit time is therefore 225.1 days.

The point of departure in heliocentric ecliptic coordinates can be calculated from a reduction of Vesta's orbital elements and the time of departure.

The point of arrival in heliocentric ecliptic coordinates can be calculated from a reduction of Earth's orbital elements and the time of arrival.

With these points (the endpoints of the intended trajectory) both known, it becomes much simpler to solve for the elements of the transfer orbit.

Had only the departure time and position been known, it would have been necessary to construct a large table of TRIAL transfer orbits from Vesta to various points on Earth's orbit, and checking each of them to see whether Earth and the spaceship had the same transit time. In general, they would not; i.e., most of the trial transfer orbits would not be actual transfer orbits between the two planets. Only when you discovered a match between transit times would you have found a valid transfer orbit.

Now that you have the arrival time, you can straightforwardly proceed to solve for the elements of the transfer orbit.

If nobody ventures a serious attempt at a solution in a few days, I'll begin solving the problem myself.

Jerry Abbott

[Jenab]

My guess is that not many people are interested in transfer orbits, and those who are intrigued by this challenge are daunted by what they think is a difficult task. But, really, if you've had vector algebra in high school and differential calculus and analytic geometry in college, you should be able to figure this out. And if you've had an astronomy course in celestial mechanics, the problem should be a cakewalk.

But, okay, for those who need a bit of jumpstarting...

The first thing to do is finding the position vector of departure and the position vector of arrival. That is, you want to calculate:

1. The position vector of Vesta, R1, at the time of departure, which is 4 February 2004, 19h 12m UT (or JD 2453040.3).

2. The position vector of Earth, R2, at the time of arrival, which is 16 September 2004, 21h 36m UT (or JD 2453265.4).

Both position vectors should be referred to heliocentric ecliptic coordinates.

Once these vectors are known, a cross product (R1xR2) will give you the vector normal to the transfer orbit, which should be divided by its own magnitude to result in a unit normal to the transfer orbit.

The inclination of the transfer orbit is found from the Z component of this unit normal vector.

The next elements to tackle are the eccentricity and the semimajor axis. Use the Pythagorean Theorem to get the line-of-light distance, d, between the point of departure and the point of arrival.

Now comes a neat trick. To find the eccentricity of the transfer orbit, solve the polar equation of the orbit (with the pole at the solar focus) simultaneously with the law of cosines. After a bit of algebra, a set of equations for the eccentricity, e, pops right out. After which, solving for the semimajor axis, a, is a trivial one-step process.

Jerry Abbott

[Tesseract]

You posted three times in a row, i would think that nobody wants to answer the question.Try a more interesting thread next time.

[jordan]

But, really, if you've had vector algebra in high school and differential calculus and analytic geometry in college, you should be able to figure this out. And if you've had an astronomy course in celestial mechanics, the problem should be a cakewalk.

Although I enjoy a good challenge, you just about summed up my problem right there. Sorry, but I have no idea what to do on this one. And don't mind Tesserect's cordial way with words.

[Jenab]

You posted three times in a row, i would think that nobody wants to answer the question.Try a more interesting thread next time.
By the time I'm finished in this thread, it will contain a mini-lesson in finding transfer orbits. Somebody might think it worthwhile, even if that somebody has not yet found it.

I've had a chance to look around this forum, and I've noticed that a great deal of the discussions are mostly cursory treatments of popular science.

I see here many of the same questions about relativity and quantum mechanics that I've seen in many other places, with the same hand-waving and lack of mathematical rigor. It's gee-whiz to the max, but the educational load in the average thread is really fairly low.

For once, let's change that.

Jerry Abbott

[Jenab]

The orbital elements that I'll use for Vesta, except the mean anomaly, came from

http://arnold.usno.navy.mil/murison/Asteroids/OrbitalElements.html

Vesta's mean anomaly at epoch was found at

http://cfa-www.harvard.edu/iau/Ephemerides/Bright/2004/00004.html

mean anomaly at epoch, M, 15.64391 degrees
epoch is JD 2453000.5 (27 Dec 2003)

Vesta's period is found by raising the semimajor axis to the power of 1.5 and multiplying the result by 365.256898326 days.

period, P, 1325.593765 days.

mean daily motion, m = 0.27157641 degrees per day

Vesta was 57.6040821 days past its perihelion at epoch, therefore...

VESTA ORBITAL ELEMENTS
time of perihelion passage, T, JD 2452942.9
semimajor axis, a, 2.36160912 astronomical units
eccentricity, e, 0.08871313
inclination, i, 7.133087 degrees
longitude of ascending node, L, 103.937391 degrees
argument of perihelion, w, 150.232154 degrees

EARTH ORBITAL ELEMENTS
time of perihelion passage, T, JD 2453009.3
semimajor axis, a, 1.00000011 astronomical units
eccentricity, e, 0.01671022
inclination, i, defined zero (ecliptic taken as plane of Earth's orbit)
longitude of ascending node, L, defined zero
argument of perihelion, w, 102.94719 degrees

If you want to convert from calendar date to Julian date, go to...
http://wwwmacho.mcmaster.ca/JAVA/JD.html

If you want to convert from Julian date to calendar date, go to...
http://wwwmacho.mcmaster.ca/JAVA/CD.html

Dates and times for Earth's perihelions, aphelions, solstices and equinoxes can be found at...
http://aa.usno.navy.mil/data/docs/EarthSeasons.html

Jerry Abbott

[Sayonara³]

You posted three times in a row, i would think that nobody wants to answer the question.Try a more interesting thread next time.
SILENCE! It's an interesting read.

[Tesseract]

SILENCE! It's an interesting read.
The kind of thing that will save the world.If you dont know this, how can you expect to be a valuable member of society? :-p :-p :-p

[Sayonara³]

What are you talking about? In fact no, don't tell me... just stop polluting this thread.

[Jenab]

Given its orbital elements and a specified moment of time, you can find where a planet was at that moment. The procedure for doing that is called a reduction of the orbital elements, and I present two examples below on how to do that. These particular examples were chosen, of course, because they are both part of the larger problem of finding a transfer orbit between the two points in space that will be the results.


THE PROCEDURE

Find the orbital period.

P = (365.256898326 days) a^1.5

Find the mean motion.

m = 2 pi radians / P

Find the time difference from perihelion.

dt = t - T

Find the mean anomaly.

M = m dt

Find the eccentric anomaly.

The equation relating the mean anomaly, the orbital eccentricity and the eccentric anomaly is known as Kepler's equation. It is transcendental in the eccentric anomaly, which is why it is necessary to use radian measure for this step. I'll use Danby's method to solve for the eccentric anomaly. Most of the time, the simpler Newton's method will also work, but Danby's method is faster and has fewer convergence failures.

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


Solving for VESTA (t = JD 2453040.3)

a = 2.36160912 astronomical units

P = 1325.594 days

m = 0.004739903 radians per day

dt = 2453040.3 - 2452942.9 = 97.4 days

M = 0.4616672 radians

e = 0.08871313

u = 0.5045526 radians

x''' = +1.857826
y''' = +1.137138
z''' = 0

w = 2.622046 radians

x'' = -2.177249
y'' = -0.06469965
z'' = 0

i = 0.1244959 radians

x' = -2.177249
y' = -0.06419891
z' = -0.008034048

L = 1.81405 radians

x = +0.5867243
y = -2.097687
z = -0.008034048


Solving for EARTH (t = JD 2453265.4)

a = 1.00000011 astronomical units

P = 365.2569 days

m = 0.0172021 radians per day

dt = 2453265.4 - 2453009.3 = 256.1 days

M = 4.405455 radians

e = 0.01671022

u = 4.389608 radians

x''' = -0.3339160
y''' = -0.9482244
z''' = 0

w = 1.796767 radians

x'' = +0.9989326
y'' = -0.1129745
z'' = 0

i = 0 radians

x' = +0.9989326
y' = -0.1129745
z' = 0

L = 0 radians

x = +0.9989326
y = -0.1129745
z = 0

Summary.

These two points in spacetime exist on the transfer orbit.

Departure from Vesta
x1 = +0.5867243
y1 = -2.097687
z1 = -0.008034048
t1 = JD 2453040.3

Arrival at Earth
x2 = +0.9989326
y2 = -0.1129745
z2 = 0
t2 = JD 2453265.4

The spatial components of these vectors are referred to heliocentric ecliptic coordinates, and the values are in astronomical units. The time component for each vector is referred to Julian Date (JD).

The accuracy of the vectors depends on the accuracy with which the orbital elements were known. Possibly, improved values could be found somewhere, but we'll go with these.

Jerry Abbott

[Jenab]

Refer to my previous post for nomenclature.

We will eventually want to know the velocity of Vesta at departure and the velocity of Earth at arrival, in order to calculate the necessary delta-vee at each end of the intended trajectory.

Let Q be the true anomaly of an object in its orbit around the sun. (I don't have access to the greek letters here.) To determine its value:

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

Where x''' and y''' are two of the components of the canonical position vector (see previous post).

The arctan2 function is the two-argument arctangent function. It is related to the single-argument 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

GMsun is the solar gravitational constant.

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

When using the equations below, make sure to keep the units consistant. You will probably want to enter the semimajor axis in meters, rather than attempt to adjust GMsun.

1 astronomical unit = 1.49597870691E+11 meters.

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

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

Vz''' = 0

In order to adjust this canonical velocity vector to heliocentric ecliptic coordinates, you'd rotate it in the same way that a position vector would be rotated. Again, see my preceding post.

Plugging in the numbers...


Orbital velocity of Vesta at JD 2453040.3

x''' = +1.857826
y''' = +1.137138

Q = 0.5492546 radians

a = 2.36160912
e = 0.08871313

Canonical velocity
Vx''' = -10158.3 m/s
Vy''' = +18322.5 m/s
Vz''' = 0

...intermediate steps skipped...

HEC velocity
Vx = +20241.4 m/s
Vy = +4735.7 m/s
Vz = -2601.2 m/s


Orbital velocity of Earth at JD 2453265.4

x''' = -0.3339160
y''' = -0.9482244

Q = 4.37380123 radians

a = 1.00000011
e = 0.01671022

Canonical velocity
Vx''' = +28097.7 m/s
Vy''' = -9396.8 m/s
Vz''' = 0

...intermediate steps skipped...

HEC velocity
Vx = +2862.5 m/s
Vy = +29488.7 m/s
Vz = 0


Jerry Abbott

[Jenab]

Given two heliocentric position vectors to points on an orbit, excepting a pair differing in true anomaly by pi radians, it is possible to determine the inclination of the orbit to the ecliptic.

Position vector 1, R1, has these components.
x1 = +0.5867243
y1 = -2.097687
z1 = -0.008034048

Position vector 2, R2, has these components.
x2 = +0.9989326
y2 = -0.1129745
z2 = 0

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)

Where [i] is the inclination of the transfer orbit to the ecliptic plane.

Plugging in the numbers...

Xn' = -0.0009076426

Yn' = -0.0080254725

Zn' = +2.0291630445

Rn' = 2.02917911803

Xn = -0.0004472955

Yn = -0.0039550340

Zn = +0.9999920788

i = 0.00398025345 radians

Jerry Abbott

[Tesseract]

You should be a professor Jenab you clearly show your talents of slowly explaining things. :-) ;)

[Jenab]

You should be a professor Jenab you clearly show your talents of slowly explaining things. :-) ;)

That's what the eMode test said after declaring me a Visionary Philosopher with a 142 IQ. :D

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

r2 = 1.005300740

d = 2.027082617

Case #3.
e = 6.287121148
a = 0.190141423
A hyperbolic orbit, with perihelion at arrival point.

Case #4.
e = 0.651493744
a = 1.318933507
An elliptical orbit, with aphelion at departure point.

It is not at all likely that both of these orbits are valid transfer orbits. One of them might be (ok, is...I've already checked), but the other will turn out to have an inappropriate transit time. A spaceship using the wrong orbit will reach the "arrival point" either much sooner or much later than the Earth does, which would be kind of foolish. The spaceship pilot will have to choose the right orbit, if he wants to remain employed.

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 in Post #12, there is another way to proceed in this special circumstance.

We found a unit vector normal to the transfer orbit Rn = [Xn , Yn, Zn] in Post #13. We also have the heliocentric position vector to the apsidal endpoint, rK = [xK , yK, zK], which was calculated in Post #11.

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.

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

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

Again, remember to keep the units consistent. You may want to convert [a] and [rK] from astronomical units to meters, for example.

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

Possible hyperbolic transfer orbit

Perihelion at arrival at Earth on JD 2453265.4
K=2

Xn = -0.0004472955
Yn = -0.0039550340
Zn = +0.9999920788

xK = +0.9989326
yK = -0.1129745
zK = 0

VxK'' = +0.11297361
VyK'' = +0.99892469
VzK'' = +0.00400135

VK'' = 1.00530074

VxK' = +0.11237792
VyK' = +0.99365757
VzK' = +0.00398025

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

1 astronomical unit = 1.49597870691E+11 meters

a = 2.8444752E+10 meters

rK = 1.50390850E+11 meters

VK = 80190.5 m/s

VxK = +9011.6 m/s
VyK = +79681.9 m/s
VzK = +319.2 m/s

If this hyperbolic orbit turns out to be an actual transfer orbit - which we don't yet know for certain - we would determine the necessary delta-vee for matching velocity with Earth.

Possible elliptical transfer orbit

Aphelion at Departure at Vesta on JD 2453040.3
K=1

Xn = -0.0004472955
Yn = -0.0039550340
Zn = +0.9999920788

xK = +0.5867243
yK = -2.097687
zK = -0.008034048

VxK'' = +2.09770216
VyK'' = +0.58671606
VzK'' = +0.00325880

VK'' = 2.17821044

VxK' = +0.96303926
VyK' = +0.26935692
VzK' = +0.00149609

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

1 astronomical unit = 1.49597870691E+11 meters

a = 1.97309644E+11 meters

rK = 3.25855643E+11 meters

VK = 11913.7 m/s

VxK = +11473.4 m/s
VyK = +3209.1 m/s
VzK = +17.8 m/s

If this elliptical orbit turns out to be an actual transfer orbit - which we don't yet know for certain - we would determine the necessary delta-vee for entering the transfer orbit from initially being at rest with respect to Vesta.

Significance of our progress so far.

We now have complete state vectors for two proposed transfer orbits at whichever endpoint of their intended trajectories occurs at an apside:

[xK, yK, zK, VxK, VyK, VzK]

The proposed hyperbolic transfer orbit has this state vector:

x2 = +0.9989326 AU
y2 = -0.1129745 AU
z2 = 0
Vx2 = +9011.6 m/s
Vy2 = +79681.9 m/s
Vz2 = +319.2 m/s

The proposed elliptical transfer orbit has this state vector:

x1 = +0.5867243 AU
y1 = -2.097687 AU
z1 = -0.008034048 AU
Vx1 = +11473.4 m/s
Vy1 = +3209.1 m/s
Vz1 = +17.8 m/s

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

Proposed hyperbolic transfer orbit

hX = -5.39471769E+12 m^2 sec^-2
hY = -4.77006702E+13 m^2 sec^-2
hZ = +1.19075189E+16 m^2 sec^-2

h = 1.19076157E+16 m^2 sec^-2

L = 6.17056860 radians

Proposed elliptical transfer orbit

hX = -1.72886746E+12 m^2 sec^-2
hY = -1.53519637E+13 m^2 sec^-2
hZ = +3.88213341E+15 m^2 sec^-2

h = 3.88216415E+15 m^2 sec^-2

L = 6.17104239 radians

Comment.

Strictly, the longitude of the ascending node for both of these proposed transfer orbits should be equal, as they both contain the same two heliocentric radius vectors and therefore occur in the same plane. The most likely reason that they are not quite equal is roundoff error. Their closeness provides a check on my work so far.

Jerry Abbott

[Jenab]

The inclination [symbolized as i] of an orbit is the angle between the plane of the orbit and the plane of the ecliptic. Herein the ecliptic plane is assumed to be the plane of Earth's orbit.

These two planes intersect in a line called "the line of nodes."

The orbit itself crosses the ecliptic at two points, called nodes. At one node, called the ascending node, the planet (or whatever the orbiting object happens to be) crosses the ecliptic with the Z component of its velocity in the direction of the North Ecliptic Pole. At the other node (the descending node), the planet crosses the ecliptic with the Z component of its velocity toward the South Ecliptic Pole.

The heliocentric ecliptic coordinate system is defined such that the sun is at the origin, with +X axis extending toward the Vernal Equinox, and +Z axis normal to the Earth's orbit. It is a right-handed system, in which the +Y axis is 90 degrees counterclockwise from the +X axis to an observer positioned somewhere up the +Z axis looking back toward the origin.

The longitude of the ascending node [L] is the angle, subtended at the sun (measured in the ecliptic counterclockwise by the aforementioned observer) from the Vernal Equinox to the ascending node.

The argument of the perihelion [w] is 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.

The true anomaly [Q] is the angle, subtended at the sun, measured in the plane of the orbit in the direction of motion, from the orbit's perihelion to the current position of the planet.

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

i = 0.00398025345 radians
L = 6.17056860 radians

cos(w+Q2) = 1.00000000 (to the limit of justified precision)

sin(w+Q2) = 0

w + Q2 = 0

This is the Case #3, having perihelion at arrival, so...

Q2 = 0

w = 0

The zero result for the argument of the perihelion should not surprise anyone. The arrival occurs at the perihelion of the transfer orbit and also in the plane of Earth's orbit - i.e., in the ecliptic. The ascending node and the perihelion are in the same spot, and thus the argument of the perihelion is zero.

The proposed elliptical transfer orbit

x1 = +0.5867243
y1 = -2.097687
z1 = -0.008034048

r1 = 2.17821044

i = 0.00398025345 radians
L = 6.17104239 radians

cos(w+Q1) = +0.37543977

sin(w+Q1) = -0.92666979

w + Q1 = 5.09732665 radians

This is the Case #4, having aphelion at departure, so...

Q1 = pi radians

w = 1.95573400 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 desire a canonical position vector for the non-apsidal endpoint of the intended trajectory. That is, 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.

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.5867243
y1 = -2.097687
z1 = -0.008034048

L = 6.17056860 radians
i = 0.00398025345 radians
w = 0

x1' = +0.81874324
y1' = -2.01846369
z1' = -0.008034048

x1'' = +0.81874324
y1'' = -2.01847968
z1'' = -8.50786288E-9

Notice how small z1'' is. Actually, it would be zero, but for roundoff error.

x1''' = +0.81874324
y1''' = -2.01847968
z1''' = -8.50786288E-9

Since w=0 for this orbit, the triple-primed vector is equal to the double-primed vector.

Q1 = 5.09773397 radians

The proposed elliptical transfer orbit

x2 = +0.9989326
y2 = -0.1129745
z2 = 0

L = 6.17104239 radians
i = 0.00398025345 radians
w = 1.95573400 radians

x2' = +1.00530063
y2' = -4.76296468E-4
z2' = 0

x2'' = +1.00530063
y2'' = -4.76292695E-4
z2'' = +1.89577565E-6

x2''' = -0.37793320
y2''' = -0.93155573
z2''' = +1.89577565E-6

Q2 = 4.32697753 radians

Jerry Abbott

[Jenab]

We have reached a point where we must begin treating hyperbolic and elliptical orbits somewhat differently.

The proposed hyperbolic transfer orbit

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

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

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

MJ = e sinh uJ - uJ

Here, the uJ is the hyperbolic eccentric anomaly in the transfer orbit of the non-apsidal endpoint of the intended trajectory, and MJ is the corresponding hyperbolic mean anomaly. When calculating MJ it is necessary that uJ be in radians.

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

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

1 astronomical unit = 1.49597870691E+11 meters

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

The hyperbolic mean motion will result to units of radians per day.

Calculation of the transit time, dt.

Once we have the transit time, we can compare it to the time difference required by the known times of departure and of arrival. If they match, or very nearly match, then we have found our transfer orbit. If they do not match, we are disgusted at having done so much work for nothing.

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

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

Note that hyperbolic orbits do not have aphelia.

The proposed elliptical transfer orbit

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 uJ is the eccentric anomaly in the transfer orbit of the non-apsidal endpoint of the intended trajectory, and MJ is the corresponding mean anomaly. When calculating MJ it is necessary that uJ be in radians.

Although you generally don’t intend to follow an elliptical transfer orbit around its complete circuit, its orbital period can be found from

P = (365.256898326 days) a^1.5

The mean motion (radians per day) of a spaceship in the transfer orbit would be

m = 2 pi radians / P

Calculation of the transit time, dt.

Once we have the transit time, we can compare it to the time difference required by the known times of departure and of arrival. If they match, or very nearly match, then we have found our transfer orbit. If they do not match, we are disgusted at having done so much work for nothing.

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

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

If necessary, correct dt to the interval [0,P) by adding or subtracting the appropriate multiple of P.

Plugging in the numbers...

Proposed hyperbolic transfer orbit

r1 = 2.17821044
e = 6.287121148
a = 0.190141423 AU
Q1 = 5.09773397 radians

u1 = -1.30600661 radians

M1 = -9.44655271 radians

m = 0.20747527 radians per day

dt = 45.5 days

Whoops! This transit time does not equal the required time difference of 225.1 days! So, although the hyperbolic orbit does connect r1 with r2, a spaceship travelling along it would go from r1 to r2 too quickly, and the destination planet would not have arrived yet. So this orbit won't work.

Proposed elliptical transfer orbit

r2 = 1.00530074
e = 0.651493744
a = 1.318933507 AU
Q2 = 4.32697753 radians

sin u2 = -0.53583329
cos u2 = +0.36494920

u2 = 5.31030870 radians
M2 = 5.84877378 radians

P = 553.26446728 days
m = 0.0113565676 radians per day

dt = 238.4 days

Hm, that's what I get for using orbital elements different from those I used in my last episode of working out this problem. Ok, the spaceship pilot sees that his trip will take a little too long, and so he applies a course correction.

Jerry Abbott

[Tesseract]

How many posts do you think are left Jerry? :-)

[Jenab]

We now find the velocity vector, VJ, in the transfer orbit at the non-apsidal endpoint of the intended trajectory.

Canonical velocity at non-apsidal endpoint for hyperbolic orbits

VxJ’’’ = -(a/rJ) { GMsun / a }^0.5 sinh uJ

VyJ’’’ = +(a/rJ) { GMsun / a }^0.5 (e^2 - 1)^0.5 cosh uJ

VzJ’’’ = 0

Canonical velocity at non-apsidal endpoint for elliptical orbits

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

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

VzJ’’’ = 0

Rotation to heliocentric ecliptic coordinates

VxJ’’ = VxJ’’’ cos w - VyJ’’’ sin w

VyJ’’ = VxJ’’’ sin w + VyJ’’’ cos w

VzJ’’ = VzJ’’’

VxJ’ = VxJ’’

VyJ’ = VyJ’’ cos i

VzJ’ = VyJ’’ sin i

VxJ = VxJ’ cos L - VyJ’ sin L

VyJ = VxJ’ sin L + VyJ’ cos L

VzJ = VzJ’

By replacing J with K in all the above, you can find the sun-relative velocity in the transfer orbit at the apsidal endpoint, also.

Finding the delta-vees on each end of the intended trajectory goes like this:

dVx1 = Vx1 - Vx(preburn)
dVy1 = Vy1 - Vy(preburn)
dVz1 = Vz1 - Vz(preburn)

dVx2 = Vx(rendezvous) - Vx2
dVy2 = Vy(rendezvous) - Vy2
dVz2 = Vz(rendezvous) - Vz2

In my example problem, V(preburn) is the orbital velocity of Vesta at the moment of departure, and V(rendezvous) is the orbital velocity of Earth at the moment of arrival.

You can work it out explicitly, if you want to. Their magnitudes will be, roughly:

dV1 = 9.2 kilometers per second

dV2 = 20.4 kilometers per second

THE END

Jerry Abbott