This is a short follow-up from earlier posts in this thread.
In post #8, I outlined a method for calculating the radial and transverse components of the initial and final velocity vectors of the Keplerian arc connecting the initial and final points, [MATH]\mathbf{R}_i[/MATH] and [MATH]\mathbf{R}_f[/MATH]. from this, I said that is is easy to calculate the fuel 3-vector form of the velocity vectors using standard coordinate transformations. While this is all true, it occurred to me that ought to be a more efficient way of carrying out these calculations. And indeed there is!
Such a method uses quaternion algebra - so, for completeness, I thought I would set out an updated method:
Step 1
As before, I'm going to assume that we know the initial and final positions along the Keplerian arc, [MATH]\mathbf{R}_i = (x_i, y_i, z_i)[/MATH] and [MATH]\mathbf{R}_f = (x_f, y_f, z_f)[/MATH], and as before we also know the semi-major axis, [math]a[/math] of the Keplerian arc; and the gravitational parameter, [math]\mu[/math].
(N.B., in the above posts, the semi-major axis was calculated from a time-of-flight equation.)
Then, calculate the following quantities:
[math]
\begin{aligned}
r_i &= \left|\mathbf{R}_i\right| \\
r_f &= \left|\mathbf{R}_f\right| \\
\\
A &= r_i + r_f \\
\ell &= \left|\mathbf{R}_f - \mathbf{R}_i\right| \\
B &= \pm\sqrt{A^2 - \ell^2}\\
\end{aligned}
[/math]
Here, [math]r_i[/math] and [math]r_f[/math] are just the distances of the initial and final points from the centre of the gravitating body; [math]A[/math] is just the sum of these two distances; and [math]\ell[/math] is the length of the chord connecting the initial and final point. The quantity [math]B[/math] doesn't have a simple geometric interpretation, but it is just a simple function of [math]A[/math] and [math]\ell[/math].
Note, though, the [math]\pm[/math]in front of the square root in the expression for [math]B[/math]. Here, we have to make a choice about the sign of [math]B[/math] that we want to use in our calculations and the way to do that is as follows: imagine that in moving from the initial point [math]\mathbf{R}_i[/math] to the final point [math]\mathbf{R}_f[/math], as seen from the gravitating body, the satellite traces out a trajectory across the sky that spans some angle [math]\theta[/math]. If that angle is between [math]0[/math] and [math]\pi[/math] (a 'short' arc), we use the positive sign to calculate [math]B[/math]. And if the angle traced out is between [math]\pi[/math] and [math]2\,\pi[/math] (a 'long' arc), we use the negative sign to calculate [math]B[/math]. In other words, by choosing the sign of [math]B[/math] in the above, we elect to examine either the 'short' arc or the 'long' arc. Without further information, both arcs are equally valid solutions to the two-point boundary value problem and we need to supply additional information about our problem if we are to choose between them.
Step 2
Before going any further, we need to check that if the given semi-major axis is positive (
i.e., an elliptical orbit) that:
[math]
a \ge \frac{1}{4}\,\left(r_i + r_f + \ell\right)
[/math]
In other words, if the orbit is elliptical, there is a minimum permissible semi-major axis. If the semi-major axis provided is positive and less then this value, this this value is not consistent with the other data provided. For hyperbolic trajectories, we require only that [math]a<0[/math].
(N.B. If the semi-major axis equals the lowest permissible value -
i.e., [math]a = (r_i + r_f + \ell)/4[/math] then the resulting trajectories correspond to the least energy transfer between the initial and final points.
Step 3
Then, carry out the following simple calculations:
[math]
\begin{aligned}
\phi_i &= \sqrt{(z_i + r_i) / 2} \\
\phi_f &= \sqrt{(z_f + r_f) / 2} \\
\end{aligned}
[/math]
[math]
\begin{aligned}
\xi_i &= \frac{x_i}{2\,\phi_i} \\
\xi_f &= \frac{x_f}{2\,\phi_f} \\
\zeta_i &= \frac{y_i}{2\,\phi_i} \\
\zeta_f &= \frac{y_f}{2\,\phi_f} \\
\end{aligned}
[/math]
[math]
\begin{aligned}
c &= \frac{2}{B}\,\left(+\xi_i\,\xi_f + \zeta_i\,\zeta_f + \phi_i\,\phi_f\right) \\
s &= \frac{2}{B}\,\left(-\xi_i\,\zeta_f + \zeta_i\,\xi_f\right) \\
\end{aligned}
[/math]
Admittedly, it's hard to see a geometric motivation for these calculations and I'm not going to try to justify or explain these expressions in this note - although I have a proof available for those sufficiently motivated to go through it. However, these computations are straightforward with the most complex arithmetic operation used being the calculation of two square roots. Note in particular that no trigonometric functions are used; and there are no coordinate transformations. All expressions involve normal real numbers only.
Step 4
Having calculated all of the above as a precursor to the main calculation, we now want to compute two quaternions as follows:
[math]
\begin{aligned}
\mathcal{Q}_i &= \xi_i\,\hat{\mathbf{i}} + \zeta_i\,\hat{\mathbf{j}} +\phi_i\,\hat{\mathbf{k}}\\
\mathcal{Q}_f &= s\,\phi_f + \left(c\,\xi_f -s\,\zeta_f\right)\,\hat{\mathbf{i}} + \left(c\,\zeta_f + s\,\xi_f\right)\,\hat{\mathbf{j}} + c\,\phi_f\,\hat{\mathbf{k}}\\
\end{aligned}
[/math]
(N.B. An arbitrary quaternion [math]\mathcal{Q}[/math] can be written in the form: [math]\mathcal{Q} = q_0\, + q_1\,\hat{\mathbf{i}} + q_2\,\hat{\mathbf{j}} + q_3\,\hat {\mathbf{k}}[/math] and it's conjugate as [math]\mathcal{Q}^\dagger = q_0\, - q_1\,\hat{\mathbf{i}} - q_2\,\hat{\mathbf{j}} - q_3\,\hat{\mathbf{k}}[/math])
Step 5
As in the earlier posts, calculate the two roots, [math]X_\pm[/math], of the quadratic, [math]A - B\,X = 2\,a\,(1 - X^2)[/math].
These two roots encode all the information that we need about the two possible Keplerian arcs that connect the initial and final points with semi-major axis [math]a[/math] - given our choice of the sign of [math]B[/math].
For parabolic orbits where we would expect a divergence in the semi-major axis [math]a[/math], we can side-step solving for the roots of the quadratic and simply set [math]X=\pm 1[/math]. Also note that for elliptical orbits [math]|X|<1[/math]; and for hyperbolic orbits [math]|X|>1[/math].
Step 6
For each root, [math]X_\pm[/math], calculate the quaternions [math]\mathcal{V}_i[/math] and [math]\mathcal{V}_f[/math]:
[math]
\begin{aligned}
\mathcal{V}_i &= \frac{1}{r_i}\,\sqrt{\frac{2\,\mu}{A - B\,X_\pm}}\,\left(\mathcal{Q}_f - X_\pm\,\mathcal{Q}_i\right)\,\hat{\mathbf{k}}\,\mathcal{Q}_i^\dagger\\
\mathcal{V}_f &= \frac{1}{r_f}\,\sqrt{\frac{2\,\mu}{A - B\,X_\pm}}\,\left(X_\pm\,\mathcal{Q}_f - \mathcal{Q}_i\right)\,\hat{\mathbf{k}}\,\mathcal{Q}_f^\dagger\\
\end{aligned}
[/math]
The resulting quaternions are 'pure' quaternions. To extract the velocity vector at the initial point, one takes the coefficients of [math]\hat{\mathbf{i}}[/math], [math]\hat{\mathbf{j}}[/math] and [math]\hat{\mathbf{k}}[/math] of [math]\mathcal{V}_i[/math] to get the x, y and z coefficients of the velocity 3-vector; and similarly for the velocity at the final point.
So, problem solved! - the full 3-vector form for the initial and final velocity vectors (and without making use of a single trigonometric function).
Some python code
This may all seem like computational gobbled-gook but, in actual fact, it is a fast and efficient way of solving the problem. To illustrate, here is some Python code that solves this problem:
Code:
import numpy as np
import quaternion
def calcs(X, Y):
aa = np.sqrt(2 * mu / (A - B * X))
V1 = aa * (Q2 - X*Q1) * K * np.conj(Q1) / ri
V2 = aa * (X*Q2 - Q1) * K * np.conj(Q2) / rf
v1 = np.array([V1.x, V1.y, V1.z])
v2 = np.array([V2.x, V2.y, V2.z])
e = np.sqrt(Y * Y + (1 - Y * Y) * (rf - ri)*(rf - ri)/(A*A - B*B))
dt = 2*np.sqrt(a*a*a/mu)*(np.arccos(X)+Y*np.sqrt(1-X*X))
print("The initial velocity vector is (m/s): ", v1)
print("The final velocity vector is (m/s): ", v2)
print("The orbital eccentricity is: ", e)
print("The transfer time is (s): ", dt)
print()
return
# input the inital and final 3-vector coordinates
Ri = np.array([-178932619.,-28063045.,-17448755.])
Rf = np.array([-178337419.,-30878143.,-18771071.])
a = 181997386.85
mu = 1.26687e17
# from these inputs, calculate some basic things
[xi, yi, zi] = Ri
[xf, yf, zf] = Rf
ri = np.sqrt(np.dot(Ri, Ri))
rf = np.sqrt(np.dot(Rf, Rf))
A = ri + rf
B = np.sqrt(2 * (np.dot(Ri, Rf) + ri * rf))
f1 = np.sqrt((zi + ri) / 2)
f2 = np.sqrt((zf + rf) / 2)
xi1 = xi / f1 / 2
xi2 = xf / f2 / 2
ze1 = yi / f1 / 2
ze2 = yf / f2 / 2
c = 2 * ( +xi1 * xi2 + ze1 * ze2 + f1 * f2) / B
s = 2 * ( -xi1 * ze2 + xi2 * ze1) / B
# define some quaternions
K = np.quaternion( 0, 0, 0, 1)
Q1 = np.quaternion( 0, xi1, ze1, f1)
Q2 = np.quaternion( s*f2, c*xi2 - s*ze2, c*ze2 + s*xi2, c*f2)
# Find the two roots of the equation A - B X = 2 a (1 - X^2)
Xp = (B + np.sqrt(16*a*a - 8*a*A + B*B)) / 4 / a
Xm = (B - np.sqrt(16*a*a - 8*a*A + B*B)) / 4 / a
# Finally, solve the problem for the first root
calcs(Xp, Xm)
When this code snippet is run, it produces the following output:
Code:
The initial velocity vector is (m/s): [ 4734.49198839 -23495.77898041 -11041.88918176]
The final velocity vector is (m/s): [ 5185.26925794 -23421.41132194 -10996.18966314]
The orbital eccentricity is: 0.0010376147570018118
The transfer time is (s): 119.9998118041033
which matches the values calculated earlier.
