
Elliptical TransferIntercept Orbits
NonHohmann Transfer Orbits that Take You Somewhere
Most of the time, when someone speaks of transfer orbits, he's referring to a special case known as a Hohmann transfer orbit. Hohmann transfer orbits have a departure occurring at one of its apsides (perihelion or aphelion) and an arrival occurring at the other apside. Thus, I could describe Hohmann transfer orbits as transfer orbits having two anchored apsides.
In practice, very few interplanetary transfer orbits will be Hohmann transfers. The orbits of solar system bodies are ellipses, not circles, and the planes in which the orbits exist are tilted with respect to each other. The direction in which a rocket must apply thrust, in order to enter a transfer orbit bound for an asteroid that is almost in the ecliptic plane might be, itself, in an angular sense, quite far from the ecliptic plane.
In this essay, I will treat a more general case of transfer orbits that have only one anchored apside, which turns out to be enough to close the equation set and permit the Keplerian elements of the transfer orbit to be found, as well as the changes of velocity required for transfer orbit insertion and, later, for matching velocity with the destination object.
In what follows, the following example problem will be used for illustration:
A spaceship is initially in Earth's orbit, but is on the opposite side of the sun from Earth. Its captain wants to enter a transfer orbit, bound for Vesta, at 12h UT on 26 June 2017. The navigator does some trial runs on a computer and discovers an elliptical transfer orbit having its aphelion at Vesta upon arrival at 4h 45m 36.036s UT on 12 June 2018. Check the navigator's work to ensure that an elliptical transfer orbit does exist for these times for departure and arrival. Show the elements of the transfer orbit and the deltavees required for transfer orbit insertion (departure) and for matching velocity with Vesta at arrival.
Spaceship initial orbit.
a = 1.000002 AU
e = 0.016711
i = 0.0°
Ω = 0.0°
ω = 103.095°
T = JD 2454285.96
Vesta's orbital elements.
a = 2.36126914 AU
e = 0.089054753
i = 7.13518389°
Ω = 103.91484282°
ω = 149.85540185°
T = JD 2454267.1969204
Departure time,
t₁ = 12h UTC, 26 June 2017
Arrival time,
t₂ = 4h 45m 36.036s UTC, 12 June 2018
It is convenient to convert t₁ and t₂ from calendar date format to Julian date format.
Converting the time of departure, t₁, from calendar date to Julian date
t₁ = 12h UTC, 26 June 2017
Y = 2017
M = 6
D = 26
Q = 12
A = 0
B = 2489909
C = 122
E = 69
F = 51
t₁ = JD 2457931.0
Converting the time of arrival, t₂, from calendar date to Julian date
t₂ = 4h 45m 36.036s UTC, 12 June 2018
Y = 2018
M = 6
D = 12
Q = 4.76001
A = 0
B = 2490274
C = 122
E = 69
F = 51
t₂ = JD 2458281.69833375
Instead of having the initial position vectors given to us, we must calculate them by reducing the elements of the spaceship's initial orbit (around the sun) and the time of departure therefrom, t₁, in order to obtain the position vector r₁, and by reducing the elements of Vesta's orbit and the time of arrival thereto, t₂, in order to obtain the position vector r₂.
For what passes below, the Sun's gravitational parameter,
GM = 1.32712440018ᴇ20 m³ sec⁻²
The ratio of the astronomical unit to the meter,
AU = 1.495978707ᴇ11 m au⁻¹
And the
Unless otherwise indicated, the coordinate system to which all unprimed vectors in this essay refer is ecliptic coordinates — heliocentric for position, and sunrelative for velocity.
Calculate the position and velocity of the spaceship in its initial orbit at the time of departure
Calculate the position and velocity of Vesta at the time of arrival
We will refer to a "hypothetical" transfer orbit until we have assured ourselves that it satisfies the condition that the calculated transit time be equal, or very nearly equal, to the required transit time.
The required transit time is the amount of time that the destination object (Vesta, in our example) takes to go from where it is at t₁ to where it is at t₂. This time difference is, of course, t₂−t₁.
The calculated transit time is the amount of time, Δt, that the spaceship takes to travel, along the hypothetical transfer orbit, from where it is at t₁ to the intersection of the hypothetical transfer orbit with the orbit of the destination object.
In general, Δt will differ substantially from t₂−t₁. It is necessary that t₁ and t₂ be chosen such that Δt is nearly equal to t₂−t₁. Once we know that to be the case, we can drop the word "hypothetical," for we will have determined that the transfer orbit does, indeed, exist.
The two known points on the hypothetical transfer orbit are
x₁ = −0.092732158 AU
y₁ = +0.979054316 AU
z₁ = 0
x₂ = −0.13298229 AU
y₂ = −2.14957848 AU
z₂ = +0.080867606 AU
The sides of the sundeparturearrival triangle are
r₁ = 0.98343612 AU
r₂ = 2.15520567 AU
d = 3.129936551 AU
We are given that the apoapsis of the hypothetical transfer orbit occurs at the arrival position, so
β = 2
φ = 1
N = −1
m₂ = u₂ = θ₂ = π radians
The eccentricity and the semimajor axis of the hypothetical transfer orbit are
e = 0.37484849
a = 1.56759505 AU
The true, eccentric, and mean anomalies in the hypothetical transfer orbit at the nonapsidal endpoint of the intended trajectory (i.e. the departure position) are
θ₁ = 0.16062918 radians
u₁ = 0.10844236 radians
m₁ = 0.067872532 radians
The period and mean motion of the hypothetical transfer orbit are
P = 716.884602 days
μ = 0.00876457005 radians/day
The calculated and required transit times in the hypothetical transfer orbit, from departure to arrival, are
Δt = 350.698335 days
t₂−t₁ = 350.69833375 days
Subject to roundoff error, the difference is about onetenth of a second. That's close enough. The transfer orbit exists.
The time of perihelion passage in the transfer orbit (although the spaceship is never there) is
T = JD 2457923.256033 = 18h 8m 41s UTC on 18 June 2017
Find the unit normal vector to the plane of the transfer orbit
Xn' = +0.079173779
Yn' = +0.0074990276
Zn' = +0.329531936
Rn' = 0.338992654
Xn = +0.233556030
Yn = +0.0221215048
Zn = +0.972091673
The transfer orbit's inclination to the ecliptic
i = 13.56812324°
Find the velocity in the transfer orbit at the apsidal endpoint of the intended trajectory
Vx₂'' = +2.091376254
Vy₂'' = −0.148158093
Vz₂'' = −0.499105248
V₂'' = 2.155205676
Vx₂' = +0.970383605
Vy₂' = −0.068744295
Vz₂' = −0.231581261
V₂ = 16041.367805 m/s
Vx₂ = +15566.280326 m/s
Vy₂ = −1102.752513 m/s
Vz₂ = −3714.880181 m/s
The angular momentum per unit mass
hx = +1.207943483ᴇ15 m²/s
hy = +1.144116363ᴇ14 m²/s
hz = +5.027623568ᴇ15 m²/s
The transfer orbit's longitude of the ascending node
Ω = 95.41068849°
Find the transfer orbit's argument of the perihelion
cos ω'' = −0.987126840
sin ω'' = +0.159939380
ω'' = 2.980963411 radians
ω' = −0.160629243 radians
ω = 350.79662233°
Keplerian elements of the transfer orbit
a = 1.56759505 AU
e = 0.37484849
i = 13.56812324°
Ω = 95.41068849°
ω = 350.79662233°
T = JD 2457923.256033
Now that you have the elements of the transfer orbit, you can calculate the changesofvelocity needed for transfer orbit insertion (departure) and for matching velocity with the target asteroid at arrival.
When we reduce the elements of the transfer orbit with the time of departure, t₁, we find that the velocity of the spaceship in the transfer orbit is
Vx₁ = −34166.4329 m/s
Vy₁ = −1690.83202 m/s
Vz₁ = +8247.34992 m/s
The velocity of the spaceship in its initial orbit at t₁ was
Vxi = −30140.9504 m/s
Vyi = −2921.69307 m/s
Vzi = 0.0 m/s
So the change of velocity required at departure for transfer orbit insertion is
ΔVx₁ = −4025.4825 m/s
ΔVy₁ = +1230.8611 m/s
ΔVz₁ = +8247.3499 m/s
ΔV₁ = 9259.4983 m/s
The velocity of Vesta, when it is intercepted by the spaceship, is
Vxf = +20933.6861 m/s
Vyf = −1766.64767 m/s
Vzf = −2490.40168 m/s
When we reduce the elements of the transfer orbit with the time of arrival, t₂, we find that the velocity of the spaceship in the transfer orbit is
Vx₂ = +15566.2801 m/s
Vy₂ = −1102.75259 m/s
Vz₂ = −3714.88014 m/s
So the change of velocity required of the spaceship at arrival to match velocity with Vesta is
ΔVx₂ = +5367.4060 m/s
ΔVy₂ = −663.8951 m/s
ΔVz₂ = +1224.4785 m/s
ΔV₂ = 5545.1917 m/s
Remember that all of the unprimed vectors in this tutorial are referred to ecliptic coordinates. If you want them in celestial coordinates (so that you can use a star chart to show you the right ascension and declination in which to point the nose of your spaceship when you apply thrust), you'll still have that to do.
Converting a velocity vector from rectangular ecliptic coordinates to spherical celestial coordinates
Let's convert the departure deltavee in our example. Here it is in ecliptic coordinates:
ΔVx₁ = −4025.4825 m/s
ΔVy₁ = +1230.8611 m/s
ΔVz₁ = +8247.3499 m/s
The obliquity of the ecliptic, ε, can be estimated from the Laskar equation.
T = (t−2451545) / 3652500
Ɛ = 84381.448″ − 4680.93″ T − 1.55″ T² + 1999.25″ T³ − 51.38″ T⁴
− 249.67″ T⁵ − 39.05″ T⁶ + 7.12″ T⁷ + 27.87″ T⁸ + 5.79″ T⁹ + 2.45″ T¹⁰
ε = Ɛ / 206264.806247
t = t₁ = JD 2457931.0
T = 0.00174839151266
Ɛ = 84373.263907"
ε = 0.409053126623 radians
The magnitude of the vector does not change by the rotation. It remains
ΔV₁' = ΔV₁ = √[ (ΔVx₁)² + (ΔVy₁)² + (ΔVz₁)² ]
ΔV₁' = 9259.4983 m/s
The velocity vector in celestial coordinates is found from
ΔVx₁' = ΔVx₁
ΔVy₁' = ΔVy₁ cos ε − ΔVz₁ sin ε
ΔVz₁' = ΔVy₁ sin ε + ΔVz₁ cos ε
ΔVx₁' = −4025.4825 m/s
ΔVy₁' = −2150.9947 m/s
ΔVz₁' = +8056.4894 m/s
The right ascension of the direction in which the velocity vector points is
α₁ = (12 hours/π) arctan( ΔVy₁' , ΔVx₁' )
α₁ = 13.8745051 hours
The declination of the direction in which the velocity vector points is
δ₁ = (180°/π) arcsin( ΔVz₁' / ΔV₁' )
δ₁ = +60.467750°
A check on the accuracy of the method by a numerical evolution of the state vector
As calculated from the Keplerian elements of the transfer orbit, at time t₁ = JD 2457931.0, the spaceship's heliocentric state vector is
x₁ = −0.092732158 AU
y₁ = +0.979054316 AU
z₁ = 0
Vx₁ = −34166.4329 m/s
Vy₁ = −1690.83202 m/s
Vz₁ = +8247.34992 m/s
As calculated from the Keplerian elements of the transfer orbit, at time t₂ = JD 2458281.69833375, the spaceship's heliocentric state vector is
x₂ = −0.13298229 AU
y₂ = −2.14957848 AU
z₂ = +0.080867606 AU
Vx₂ = +15566.2801 m/s
Vy₂ = −1102.75259 m/s
Vz₂ = −3714.88014 m/s
The transit time of the spaceship in the transfer orbit is
t₂−t₁ = 30300336.036 sec
If we take the state vector at t₁ and numerically walk it forward in time by 1 second intervals for 30300336 seconds, we get this result.
x = −0.132983124 AU
y = −2.149579643 AU
z = +0.080867833 AU
Vx = +15566.27354 m/s
Vy = −1102.76889 m/s
Vz = −3714.87818 m/s
Most of the time, when someone speaks of transfer orbits, he's referring to a special case known as a Hohmann transfer orbit. Hohmann transfer orbits have a departure occurring at one of its apsides (perihelion or aphelion) and an arrival occurring at the other apside. Thus, I could describe Hohmann transfer orbits as transfer orbits having two anchored apsides.
In practice, very few interplanetary transfer orbits will be Hohmann transfers. The orbits of solar system bodies are ellipses, not circles, and the planes in which the orbits exist are tilted with respect to each other. The direction in which a rocket must apply thrust, in order to enter a transfer orbit bound for an asteroid that is almost in the ecliptic plane might be, itself, in an angular sense, quite far from the ecliptic plane.
In this essay, I will treat a more general case of transfer orbits that have only one anchored apside, which turns out to be enough to close the equation set and permit the Keplerian elements of the transfer orbit to be found, as well as the changes of velocity required for transfer orbit insertion and, later, for matching velocity with the destination object.
In what follows, the following example problem will be used for illustration:
A spaceship is initially in Earth's orbit, but is on the opposite side of the sun from Earth. Its captain wants to enter a transfer orbit, bound for Vesta, at 12h UT on 26 June 2017. The navigator does some trial runs on a computer and discovers an elliptical transfer orbit having its aphelion at Vesta upon arrival at 4h 45m 36.036s UT on 12 June 2018. Check the navigator's work to ensure that an elliptical transfer orbit does exist for these times for departure and arrival. Show the elements of the transfer orbit and the deltavees required for transfer orbit insertion (departure) and for matching velocity with Vesta at arrival.
Spaceship initial orbit.
a = 1.000002 AU
e = 0.016711
i = 0.0°
Ω = 0.0°
ω = 103.095°
T = JD 2454285.96
Vesta's orbital elements.
a = 2.36126914 AU
e = 0.089054753
i = 7.13518389°
Ω = 103.91484282°
ω = 149.85540185°
T = JD 2454267.1969204
Departure time,
t₁ = 12h UTC, 26 June 2017
Arrival time,
t₂ = 4h 45m 36.036s UTC, 12 June 2018
It is convenient to convert t₁ and t₂ from calendar date format to Julian date format.
Converting from Calendar Date to Julian Date
After Fliegel and van Flandern (1968).
The time zone must be Greenwich, Zulu, UT, UTC (all the same zone)
Y = the fourdigit year
M = the month of the year (1=January... 12=December)
D = the day of the month
Q = the time of the day in decimal hours
A = integer [ (M−14) / 12 ]
B = integer { [ 1461 (Y + 4800 + A) ] / 4 }
C = integer { [ 367 (M − 2 − 12A) ] / 12 }
E = integer [ (Y + 4900 + A) / 100 ]
F = integer [ (3E) / 4 ]
t = B + C − F + D − 32075.5 + Q/24
Converting the time of departure, t₁, from calendar date to Julian date
t₁ = 12h UTC, 26 June 2017
Y = 2017
M = 6
D = 26
Q = 12
A = 0
B = 2489909
C = 122
E = 69
F = 51
t₁ = JD 2457931.0
Converting the time of arrival, t₂, from calendar date to Julian date
t₂ = 4h 45m 36.036s UTC, 12 June 2018
Y = 2018
M = 6
D = 12
Q = 4.76001
A = 0
B = 2490274
C = 122
E = 69
F = 51
t₂ = JD 2458281.69833375
Instead of having the initial position vectors given to us, we must calculate them by reducing the elements of the spaceship's initial orbit (around the sun) and the time of departure therefrom, t₁, in order to obtain the position vector r₁, and by reducing the elements of Vesta's orbit and the time of arrival thereto, t₂, in order to obtain the position vector r₂.
For what passes below, the Sun's gravitational parameter,
GM = 1.32712440018ᴇ20 m³ sec⁻²
The ratio of the astronomical unit to the meter,
AU = 1.495978707ᴇ11 m au⁻¹
And the
Definition of the twodimensional arctangent function.
atn(z) = single argument arctangent function of the argument z.
Function arctan( y , x )
. if x = 0 and y greater than 0 then angle = +π/2
. if x = 0 and y = 0 then angle = 0
. if x = 0 and y less than 0 then angle = −π/2
. if x greater than 0 and y greater than 0 then angle = atn(y/x)
. if x less than 0 then angle = atn(y/x) + π
. if x greater than 0 and y less than 0 then angle = atn(y/x) + 2π
arctan = angle
Unless otherwise indicated, the coordinate system to which all unprimed vectors in this essay refer is ecliptic coordinates — heliocentric for position, and sunrelative for velocity.
Reducing Keplerian orbital elements and a time to position and velocity in heliocentric ecliptic coordinates
Find the period, P, in days.
P = (365.256898326 days) a¹·⁵
Find the mean anomaly, m, in radians.
m₀ = (t − T) / P
m = 2π [ m₀ − integer(m₀) ]
Find the eccentric anomaly, u, in radians.
The Danby first approximation for the eccentric anomaly, u, in radians.
u' = m
+ (e − e³/8 + e⁵/192) sin(m)
+ (e²/2 − e⁴/6) sin(2m)
+ (3e³/8 − 27e⁵/128) sin(3m)
+ (e⁴/3) sin(4m)
The Danby's method refinement for the eccentric anomaly.
u = u'
REPEAT
U = u
F₀ = U − e sin U − m
F₁ = 1 − e cos U
F₂ = e sin U
F₃ = e cos U
D₁ = −F₀ / F₁
D₂ = −F₀ / [ F₁ + D₁F₂/2 ]
D₃ = −F₀ / [ F₁ + D₁F₂/2 + D₂²F₃/6 ]
u = U + D₃
UNTIL u−U is less than 1ᴇ14
The loop, just above, converges u to the correct value of the eccentric anomaly. Usually. However, when e is near one and the orbiting object is near the periapsis of its orbit, there is a chance that this loop will fail to converge. In such cases, a different rootfinding method will be needed.
Find the canonical position vector of the object in its orbit at time t.
x''' = a (cos u − e)
y''' = a sin u √(1−e²)
z''' = 0
Find the true anomaly, θ. We'll use it below when we find the velocity.
θ = arctan( y''' , x''' )
Rotate the tripleprime position vector by the argument of the perihelion, ω.
x'' = x''' cos ω − y''' sin ω
y'' = x''' sin ω + y''' cos ω
z'' = z''' = 0
Rotate the doubleprime position vector by the inclination, i.
x' = x''
y' = y'' cos i
z' = y'' sin i
Rotate the singleprime position vector by the longitude of the ascending node, Ω.
x = x' cos Ω − y' sin Ω
y = x' sin Ω + y' cos Ω
z = z'
The unprimed position vector [x,y,z] is the position in heliocentric ecliptic coordinates.
Find the canonical (tripleprime) heliocentric velocity vector.
k = √{ GM / [ a AU (1 − e²) ] }
k is a speed in meters per second.
Vx''' = −k sin θ
Vy''' = k (e + cos θ)
Vz''' = 0
Rotate the tripleprime velocity vector by the argument of the perihelion, ω.
Vx'' = Vx''' cos ω − Vy''' sin ω
Vy'' = Vx''' sin ω + Vy''' cos ω
Vz'' = Vz''' = 0
Rotate the doubleprime velocity vector by the inclination, i.
Vx' = Vx''
Vy' = Vy'' cos i
Vz' = Vy'' sin i
Rotate the singleprime velocity vector by the longitude of the ascending node, Ω.
Vx = Vx' cos Ω − Vy' sin Ω
Vy = Vx' sin Ω + Vy' cos Ω
Vz = Vz'
The unprimed velocity vector [Vx,Vy,Vz] is the sunrelative velocity in ecliptic coordinates.
Calculate the position and velocity of the spaceship in its initial orbit at the time of departure
P = 365.257994  y'' = +0.979054316  Vy''' = +30019.1146 
m₀ = 9.97935722  x' = −0.092732158  Vx'' = −30140.9504 
m = 6.15348288  y' = +0.979054316  Vy'' = −2921.69307 
u' = 6.15128508  z' = 0  Vx' = −30140.9504 
u = 6.15128508  xi = −0.092732158  Vy' = −2921.69307 
x''' = +0.974604719  yi = +0.979054316  Vz' = 0 
y''' = −0.131499998  zi = 0  Vxi = −30140.9504 
θ = 6.14906877  k = 29788.8217  Vyi = −2921.69307 
x'' = −0.092732158  Vx''' = +3983.20734  Vzi = 0 
Calculate the position and velocity of Vesta at the time of arrival
P = 1325.30752  y'' = +0.651051227  Vy''' = +20727.481 
m₀ = 3.02910935  x' = −2.0545179  Vx'' = −6748.92645 
m = 0.182899417  y' = +0.646009389  Vy'' = −20049.7967 
u' = 0.200646945  z' = +0.080867606  Vx' = −6748.92645 
u = 0.200648459  xf = −0.13298229  Vy' = −19894.5281 
x''' = +2.10361404  yf = −2.14957848  Vz' = −2490.40168 
y''' = +0.468742457  zf = +0.080867606  Vxf = +20933.6861 
θ = 0.219245394  k = 19460.2928  Vyf = −1766.64767 
x'' = −2.0545179  Vx''' = −4232.48025  Vzf = −2490.40168 
We will refer to a "hypothetical" transfer orbit until we have assured ourselves that it satisfies the condition that the calculated transit time be equal, or very nearly equal, to the required transit time.
The required transit time is the amount of time that the destination object (Vesta, in our example) takes to go from where it is at t₁ to where it is at t₂. This time difference is, of course, t₂−t₁.
The calculated transit time is the amount of time, Δt, that the spaceship takes to travel, along the hypothetical transfer orbit, from where it is at t₁ to the intersection of the hypothetical transfer orbit with the orbit of the destination object.
In general, Δt will differ substantially from t₂−t₁. It is necessary that t₁ and t₂ be chosen such that Δt is nearly equal to t₂−t₁. Once we know that to be the case, we can drop the word "hypothetical," for we will have determined that the transfer orbit does, indeed, exist.
The determination of an elliptical transferintercept orbit from a position and time of departure and from a position and time of arrival
At time t₁ a spaceship in free orbit around the sun (i.e. there is no planet nearby) has this state vector:
xi , yi , zi , Vxi , Vyi , Vzi
At time t₂ (such that t₂>t₁) an asteroid in free orbit around the sun (i.e. there is no significant perturbing third mass) has this state vector:
xf , yf , zf , Vxf , Vyf , Vzf
We want to find out whether or not there exists a transfer orbit between the position elements of those two state vectors, such that
x₁ = xi
y₁ = yi
z₁ = zi
x₂ = xf
y₂ = yf
z₂ = zf
The subscript 1 denotes "pertaining to the transfer orbit at transfer orbit insertion" or "departure."
The subscript 2 denotes "pertaining to the transfer orbit at its intersection with the destination object's orbit" or "arrival." (Whether the destination object is actually there at that same time is a question we will answer presently.)
r₁ = √[ x₁² + y₁² + z₁² ]
r₂ = √[ x₂² + y₂² + z₂² ]
d = √[ (x₂−x₁)² + (y₂−y₁)² + (z₂−z₁)² ]
We define the integer variable β and permit it to have only the values 1 and 2.
If β=1, an apside (perihelion or aphelion) of the transfer orbit occurs at departure.
If β=2, an apside (perihelion or aphelion) of the transfer orbit occurs at arrival.
β = either 1 or 2
φ = 3 − β
N = (−1)ᵠ
The variables β and φ will usually be subscripts. The variable N is a sign toggle factor.
m : mean anomaly
u : eccentric anomaly
θ : true anomaly
If the apside at the apsidal endpoint of the intended trajectory is the perihelion, then
mᵦ = uᵦ = θᵦ = 0
If the apside at the apsidal endpoint of the intended trajectory is the aphelion, then
mᵦ = uᵦ = θᵦ = π radians
The eccentricity of a conic section, having the sun at a focus, which includes the point of departure and the point of arrival, is found by solving, simultaneously,
The polar equations which relates the heliocentric distances with the true anomalies,
r₁ = a (1−e²) / (1 + e cos θ₁)
r₂ = a (1−e²) / (1 + e cos θ₂)
The law of cosines,
d² = r₁² + r₂² − 2 r₁ r₂ cos(θ₂−θ₁)
In those equations, the known quantities are d, r₁, r₂, and whichever of the true anomalies, θᵦ, has the hypothetical transfer orbit's apside. The cosine of θᵦ is +1 if that apside is the perihelion, and it is −1 if it is the aphelion. The condition that one of the endpoints of the intended trajectory, either departure or arrival, occurs at one of the apsides of the hypothetical transfer orbit is why you don't need a third point on the transfer orbit to determine its elements. The three equations have three unknown quantities, namely a, e, and θᵩ.
After some algebra, we get the eccentricity of the hypothetical transfer orbit.
e = 2 (cos θᵦ) rᵦ (rᵦ−rᵩ) / (rᵩ² − rᵦ² − d²)
The semimajor axis of the hypothetical transfer orbit is found from
a = rᵦ / (1 − e cos θᵦ)
The true anomaly in the hypothetical transfer orbit at the nonapsidal endpoint of the intended trajectory is found as follows:
θᵩ = θᵦ + N arccos{(rᵦ² + rᵩ² − d²) / (2rᵦrᵩ)}
The eccentric anomaly in the hypothetical transfer orbit at the nonapsidal endpoint of the intended trajectory is found as follows:
sin uᵩ = (rᵩ/a) sin θᵩ / √(1−e²)
cos uᵩ = (rᵩ/a) cos θᵩ + e
uᵩ = arctan(sin uᵩ , cos uᵩ)
The mean anomaly in the hypothetical transfer orbit at the nonapsidal endpoint of the intended trajectory is found as follows:
mᵩ = uᵩ − e sin uᵩ
The period of the hypothetical transfer orbit is
P = (365.256898326 days) a¹·⁵
The mean motion in the hypothetical transfer orbit is
μ = 2π/P
For short path trajectories (for which the arc of true anomaly going from departure to arrival is less than π radians), the calculated transit time in the hypothetical transfer orbit is
Δt = (N/μ) [mᵩ − π sin(θᵦ/2)]
Here's the test that determines whether a transfer orbit exists between heliocentric position r₁ at time t₁ and heliocentric position r₂ at time t₂. It is necessary that
Δt ≈ t₂ − t₁
And the match should be a close one, ideally a small fraction of a second. In general, this will not be the case. If the difference in the required and the calculated transit times is unacceptably large, then the spaceship pilot will have to choose either a different departure time, or a different arrival time, or both, and try again.
The procedure being demonstrated here finds elliptical transfer orbits of the short path, by which it is meant that the arc of true anomaly along the intended trajectory, from departure to arrival, is strictly less than π radians. One of the transfer orbit's apsides will occur at either the position of departure or at the position of arrival, but the other apside will not occur at all within the intended trajectory. To be complete about things, we will calculate the time of perihelion passage in the transfer orbit, whether or not the spaceship is ever there.
T = tᵦ − mᵦ/μ
The inclination of the transfer orbit is found from the cross product of the heliocentric position vectors of departure and arrival, r₁ x r₂.
Xn' = y₁ z₂ − z₁ y₂
Yn' = z₁ x₂ − x₁ z₂
Zn' = x₁ y₂ − y₁ x₂
Rn' = √[ (Xn')² + (Yn')² + (Zn')² ]
Xn = Xn' / Rn'
Yn = Yn' / Rn'
Zn = Zn' / Rn'
The vector Rn is a unit normal to the transfer orbit in the direction of the orbit's angular momentum.
i = arccos(Zn)
Having found the components of the vector normal to the transfer orbit (in the direction of the angular momentum), we now use it to find the velocity in the transfer orbit at the apsidal endpoint of the intended trajectory.
Vxᵦ'' = Yn zᵦ − Zn yᵦ
Vyᵦ'' = Zn xᵦ − Xn zᵦ
Vzᵦ'' = Xn yᵦ − Yn xᵦ
Vᵦ'' = √[ (Vxᵦ'')² + (Vyᵦ'')² + (Vzᵦ'')² ]
Vxᵦ' = Vxᵦ'' / Vᵦ''
Vyᵦ' = Vyᵦ'' / Vᵦ''
Vzᵦ' = Vzᵦ'' / Vᵦ''
Vᵦ = √[ (GM/AU) (2/rᵦ − 1/a) ]
Vxᵦ = Vᵦ Vxᵦ'
Vyᵦ = Vᵦ Vyᵦ'
Vzᵦ = Vᵦ Vzᵦ'
Now we find the angular momentum per unit mass in the transfer orbit.
hx = AU (yᵦ Vzᵦ − zᵦ Vyᵦ)
hy = AU (zᵦ Vxᵦ − xᵦ Vzᵦ)
hz = AU (xᵦ Vyᵦ − yᵦ Vxᵦ)
From here, we find the longitude of the ascending node of the transfer orbit.
Ω = arctan( hx , −hy )
Notice that hx is proportional to sin Ω, while −hy is proportional to cos Ω.
We find the argument of the perihelion of the transfer orbit as follows:
cos ω'' = (xᵦ cos Ω + yᵦ sin Ω) / rᵦ
If sin i = 0 then sin ω'' = (yᵦ cos Ω − xᵦ sin Ω) / rᵦ
If sin i ≠ 0 then sin ω'' = zᵦ / (rᵦ sin i)
ω'' = arctan( sin ω'' , cos ω'' )
ω' = ω'' − θᵦ
If ω' ≥ 0 then ω = ω'
If ω' < 0 then ω = ω' + 2π
The two known points on the hypothetical transfer orbit are
x₁ = −0.092732158 AU
y₁ = +0.979054316 AU
z₁ = 0
x₂ = −0.13298229 AU
y₂ = −2.14957848 AU
z₂ = +0.080867606 AU
The sides of the sundeparturearrival triangle are
r₁ = 0.98343612 AU
r₂ = 2.15520567 AU
d = 3.129936551 AU
We are given that the apoapsis of the hypothetical transfer orbit occurs at the arrival position, so
β = 2
φ = 1
N = −1
m₂ = u₂ = θ₂ = π radians
The eccentricity and the semimajor axis of the hypothetical transfer orbit are
e = 0.37484849
a = 1.56759505 AU
The true, eccentric, and mean anomalies in the hypothetical transfer orbit at the nonapsidal endpoint of the intended trajectory (i.e. the departure position) are
θ₁ = 0.16062918 radians
u₁ = 0.10844236 radians
m₁ = 0.067872532 radians
The period and mean motion of the hypothetical transfer orbit are
P = 716.884602 days
μ = 0.00876457005 radians/day
The calculated and required transit times in the hypothetical transfer orbit, from departure to arrival, are
Δt = 350.698335 days
t₂−t₁ = 350.69833375 days
Subject to roundoff error, the difference is about onetenth of a second. That's close enough. The transfer orbit exists.
The time of perihelion passage in the transfer orbit (although the spaceship is never there) is
T = JD 2457923.256033 = 18h 8m 41s UTC on 18 June 2017
Find the unit normal vector to the plane of the transfer orbit
Xn' = +0.079173779
Yn' = +0.0074990276
Zn' = +0.329531936
Rn' = 0.338992654
Xn = +0.233556030
Yn = +0.0221215048
Zn = +0.972091673
The transfer orbit's inclination to the ecliptic
i = 13.56812324°
Find the velocity in the transfer orbit at the apsidal endpoint of the intended trajectory
Vx₂'' = +2.091376254
Vy₂'' = −0.148158093
Vz₂'' = −0.499105248
V₂'' = 2.155205676
Vx₂' = +0.970383605
Vy₂' = −0.068744295
Vz₂' = −0.231581261
V₂ = 16041.367805 m/s
Vx₂ = +15566.280326 m/s
Vy₂ = −1102.752513 m/s
Vz₂ = −3714.880181 m/s
The angular momentum per unit mass
hx = +1.207943483ᴇ15 m²/s
hy = +1.144116363ᴇ14 m²/s
hz = +5.027623568ᴇ15 m²/s
The transfer orbit's longitude of the ascending node
Ω = 95.41068849°
Find the transfer orbit's argument of the perihelion
cos ω'' = −0.987126840
sin ω'' = +0.159939380
ω'' = 2.980963411 radians
ω' = −0.160629243 radians
ω = 350.79662233°
Keplerian elements of the transfer orbit
a = 1.56759505 AU
e = 0.37484849
i = 13.56812324°
Ω = 95.41068849°
ω = 350.79662233°
T = JD 2457923.256033
Now that you have the elements of the transfer orbit, you can calculate the changesofvelocity needed for transfer orbit insertion (departure) and for matching velocity with the target asteroid at arrival.
When we reduce the elements of the transfer orbit with the time of departure, t₁, we find that the velocity of the spaceship in the transfer orbit is
Vx₁ = −34166.4329 m/s
Vy₁ = −1690.83202 m/s
Vz₁ = +8247.34992 m/s
The velocity of the spaceship in its initial orbit at t₁ was
Vxi = −30140.9504 m/s
Vyi = −2921.69307 m/s
Vzi = 0.0 m/s
So the change of velocity required at departure for transfer orbit insertion is
ΔVx₁ = −4025.4825 m/s
ΔVy₁ = +1230.8611 m/s
ΔVz₁ = +8247.3499 m/s
ΔV₁ = 9259.4983 m/s
The velocity of Vesta, when it is intercepted by the spaceship, is
Vxf = +20933.6861 m/s
Vyf = −1766.64767 m/s
Vzf = −2490.40168 m/s
When we reduce the elements of the transfer orbit with the time of arrival, t₂, we find that the velocity of the spaceship in the transfer orbit is
Vx₂ = +15566.2801 m/s
Vy₂ = −1102.75259 m/s
Vz₂ = −3714.88014 m/s
So the change of velocity required of the spaceship at arrival to match velocity with Vesta is
ΔVx₂ = +5367.4060 m/s
ΔVy₂ = −663.8951 m/s
ΔVz₂ = +1224.4785 m/s
ΔV₂ = 5545.1917 m/s
Remember that all of the unprimed vectors in this tutorial are referred to ecliptic coordinates. If you want them in celestial coordinates (so that you can use a star chart to show you the right ascension and declination in which to point the nose of your spaceship when you apply thrust), you'll still have that to do.
Converting a velocity vector from rectangular ecliptic coordinates to spherical celestial coordinates
Let's convert the departure deltavee in our example. Here it is in ecliptic coordinates:
ΔVx₁ = −4025.4825 m/s
ΔVy₁ = +1230.8611 m/s
ΔVz₁ = +8247.3499 m/s
The obliquity of the ecliptic, ε, can be estimated from the Laskar equation.
T = (t−2451545) / 3652500
Ɛ = 84381.448″ − 4680.93″ T − 1.55″ T² + 1999.25″ T³ − 51.38″ T⁴
− 249.67″ T⁵ − 39.05″ T⁶ + 7.12″ T⁷ + 27.87″ T⁸ + 5.79″ T⁹ + 2.45″ T¹⁰
ε = Ɛ / 206264.806247
t = t₁ = JD 2457931.0
T = 0.00174839151266
Ɛ = 84373.263907"
ε = 0.409053126623 radians
The magnitude of the vector does not change by the rotation. It remains
ΔV₁' = ΔV₁ = √[ (ΔVx₁)² + (ΔVy₁)² + (ΔVz₁)² ]
ΔV₁' = 9259.4983 m/s
The velocity vector in celestial coordinates is found from
ΔVx₁' = ΔVx₁
ΔVy₁' = ΔVy₁ cos ε − ΔVz₁ sin ε
ΔVz₁' = ΔVy₁ sin ε + ΔVz₁ cos ε
ΔVx₁' = −4025.4825 m/s
ΔVy₁' = −2150.9947 m/s
ΔVz₁' = +8056.4894 m/s
The right ascension of the direction in which the velocity vector points is
α₁ = (12 hours/π) arctan( ΔVy₁' , ΔVx₁' )
α₁ = 13.8745051 hours
The declination of the direction in which the velocity vector points is
δ₁ = (180°/π) arcsin( ΔVz₁' / ΔV₁' )
δ₁ = +60.467750°
A check on the accuracy of the method by a numerical evolution of the state vector
As calculated from the Keplerian elements of the transfer orbit, at time t₁ = JD 2457931.0, the spaceship's heliocentric state vector is
x₁ = −0.092732158 AU
y₁ = +0.979054316 AU
z₁ = 0
Vx₁ = −34166.4329 m/s
Vy₁ = −1690.83202 m/s
Vz₁ = +8247.34992 m/s
As calculated from the Keplerian elements of the transfer orbit, at time t₂ = JD 2458281.69833375, the spaceship's heliocentric state vector is
x₂ = −0.13298229 AU
y₂ = −2.14957848 AU
z₂ = +0.080867606 AU
Vx₂ = +15566.2801 m/s
Vy₂ = −1102.75259 m/s
Vz₂ = −3714.88014 m/s
The transit time of the spaceship in the transfer orbit is
t₂−t₁ = 30300336.036 sec
If we take the state vector at t₁ and numerically walk it forward in time by 1 second intervals for 30300336 seconds, we get this result.
x = −0.132983124 AU
y = −2.149579643 AU
z = +0.080867833 AU
Vx = +15566.27354 m/s
Vy = −1102.76889 m/s
Vz = −3714.87818 m/s