Orbiter-Forum  

Go Back   Orbiter-Forum > Far Side of the Moon > Math & Physics
Register Blogs Orbinauts List Social Groups FAQ Projects Mark Forums Read

Math & Physics Mathematical and physical problems of space flight and astronomy.

Reply
 
Thread Tools
Old 06-02-2019, 01:35 AM   #1
ncc1701d
Orbinaut
Default Can I get the Vector components xyz given only this information?

hello,
I have these 3 dates
and I know the position X Y Z of moon Amalthea relative to Jupiter
and then I know the Velocity Vectors V for each one.
I need to get the VX VY VZ components for each of the 3 date-times.
Do I have enough information to get those from what you see?
I believe I need the instantaneous values.
I am aware of VX = dx/dt, VY = dy/dt, VZ = dz/dt
and I know VX + VY + VZ = V
I am interested in the math behind finding those components given what I have.
If you know the procedure can you tell me what it is? or shed some light on how I might solve the problem if its possible?
thank you.

1979 MAR 05 00:00:00
-178932.619 km -28063.045 km -17448.755 km V = -29.804114 km/sec

1979 MAR 05 00:02:00
-178337.419 km -30878.143 km -18771.071 km V = -29.231286 km/sec

1979 MAR 05 00:04:00
-177688.033 km -33683.859 km -20087.682 km V = 28.649604 km/sec

thank you.
I do have the answers if someone knew how to do it. I can check to see if your right.

Last edited by ncc1701d; 06-02-2019 at 01:45 AM.
ncc1701d is offline   Reply With Quote
Old 06-02-2019, 01:48 AM   #2
Quick_Nick
Passed the Turing Test
 
Quick_Nick's Avatar
Default

If I understand correctly, you have XYZ position information and a speed, and you want XYZ velocity information.

This can get really complicated if you care about a perfect answer.

However I would say simply knowing that the orbit is nearly circular and nearly equatorial is enough information to have a pretty good idea of the velocity vector.
You need to know what coordinate frame those XYZ values are in though.

Edit: I'm thinking you may actually care about a more accurate answer given your three times are less than an orbit apart.
I suspect you can in fact evaluate those three points together to get some orbital parameters out.

Last edited by Quick_Nick; 06-02-2019 at 01:55 AM.
Quick_Nick is offline   Reply With Quote
Old 06-03-2019, 01:16 PM   #3
Marijn
Orbinaut
 
Marijn's Avatar
Default

I think this article describes the required math. It's not easy. https://science.larouchepac.com/gaus...erProblem.html
Marijn is offline   Reply With Quote
Thanked by:
Old 06-06-2019, 02:32 PM   #4
BrianJ
Addon Developer
Default

Hi,
thanks for the interesting problem to think about :-)

I can get the semimajor axis (a) and the standard gravitational parameter (mu) from the vis-viva equation;

v^2 = mu*(2/r - 1/a)

But then I get stuck!

Then I thought maybe it was an energy-conservation problem (Potential + Kinetic = Constant) but I'm not getting anywhere with that.

Since the times of the data points are given, I'm pretty sure that is necessary information. Needs some kind of insight regarding Kepler's 2nd law (radius sweeps out equal areas in equal times) maybe?

I think that the shape of an ellipse is completely determined by three points (pos vectors), so theoretically you should be able to find the equation of the ellipse and the tangent at a given point - but the math for doing that in 3D is looking really gnarly, so I'm not even thinking about that :-)

Anyway, I'm having fun playing with this one.
Cheers,
Brian

P.S. Your first two Velocity values are negative - typo?
BrianJ is offline   Reply With Quote
Old 06-06-2019, 06:08 PM   #5
ncc1701d
Orbinaut
Default

the checkmarked answer I found here might provide the answer.
https://stackoverflow.com/questions/...elocity-vector

Although not sure its that simple yet. I am still researching.
ncc1701d is offline   Reply With Quote
Old 07-12-2019, 02:35 PM   #6
MontBlanc2012
Orbinaut
Default

Quote:
Originally Posted by BrianJ View Post
 Hi,
thanks for the interesting problem to think about :-)

I can get the semimajor axis (a) and the standard gravitational parameter (mu) from the vis-viva equation;

v^2 = mu*(2/r - 1/a)

But then I get stuck!
I'm a bit suspicious of the orbital speeds given in the opening post, but here is the general procedure for solving this problem for Keplerian orbits where the points are separated by less than half an orbit:

Let's suppose, as above, that you are given two (3-vector) points on an Keplerian orbit \mathbf{R}_i and \mathbf{R}_f and that you know the (scalar) orbital speeds at those points, V_i and V_f respectively.

Now calculate:

\begin{aligned}r_i &= \sqrt{\mathbf{R}_i.\mathbf{R}_i} \\r_f &= \sqrt{\mathbf{R}_f.\mathbf{R}_f} \\\theta &=\cos^{-1}\left(\frac{\mathbf{R}_i.\mathbf{R}_f}{r_i\,r_f}\right) \\A &= r_i + r_f \\B &= 2\,\sqrt{r_i\,r_f}\,\cos\left(\frac{\theta}{2}\right)\end{aligned}

In addition, calculate the semi-major axis at either the first or last point using:

V_i^2 = \mu\,\left(\frac{2}{r_i}-\frac{1}{a}\right)

or

V_f^2 = \mu\,\left(\frac{2}{r_f}-\frac{1}{a}\right)

Then (and here's the magic!), solve for the two roots of the quadratic equation:

A-B\,X = 2\,a\,(1 - X^2)

Call these roots X_1 and X_2. These two numbers encode information about the two Keplerian arcs that join the points \mathbf{R}_i and \mathbf{R}_f with semi-major axis a.

Then, the radial and transverse components of the velocity vector at the \mathbf{R}_i are given by:

\left(v_{i,r},v_{i,\theta}\right) = \sqrt{\frac{2\,\mu}{A-B\,X}}\,\left(\sqrt{\frac{r_f}{r_i}}\,\cos\left(\frac{\theta}{2}\right)- X,\sqrt{\frac{r_f}{r_i}}\,\sin\left(\frac{\theta}{2}\right)\right)

and the radial and transverse components of the velocity vector at the \mathbf{R}_f are given by:

\left(v_{f,r},v_{f,\theta}\right) = \sqrt{\frac{2\,\mu}{A-B\,X}}\,\left(X - \sqrt{\frac{r_i}{r_f}}\,\cos\left(\frac{\theta}{2}\right),\sqrt{\frac{r_i}{r_f}}\,\sin\left(\frac{\theta}{2}\right)\right)

As mentioned, this gives two possible solutions - corresponding to the two possible Keplerian arcs that pass through the initial and final points with set-major-axis a. For each solution, we calculate the eccentricity in the usual way. In practice, for problems involving nearly circular orbits, one value of X will give a high eccentricity; and the other will give one close to zero. Obviously, we choose the solutions corresponding to the low orbital eccentricity - thereby removing the degeneracy.

(N.B., I have a proof of the above but it's a little long-winded - too long for this post. But I can share if anyone is sufficiently interested)

---------- Post added at 02:35 PM ---------- Previous post was at 01:23 PM ----------

Let's take a simple worked example:

Suppose that we choose units such that \mu=1 with:

\begin{aligned}
\mathbf{R}_i &= \left(-0.9803802, 0.241571, 0\right) \\
\mathbf{R}_f &= \left(-0.9972285, -0.159303, 0\right) \\
\\
V_i &= 0.9903428 \\
V_f &= 0.9901760
\end{aligned}

Find the velocity vectors at \mathbf{R}_i and \mathbf{R}_f.

Using the above algorithm we calculate:

\begin{aligned}
r_i &= \sqrt{\mathbf{R}_i.\mathbf{R}_i} = 1.009704\\
r_f &= \sqrt{\mathbf{R}_f.\mathbf{R}_f} = 1.009872 \\
\theta &=\cos^{-1}\left(\frac{\mathbf{R}_i.\mathbf{R}_f}{r_i\,r_f}\right) = 0.4000000\\
A &= r_i + r_f = 2.019576\\
B &= 2\,\sqrt{r_i\,r_f}\,\cos\left(\frac{\theta}{2}\right) = 1.979319
\end{aligned}

We next calculate the semi-major axis using

V_i^2 = \mu\,\left(\frac{2}{r_i}-\frac{1}{a}\right)

or

V_f^2 = \mu\,\left(\frac{2}{r_f}-\frac{1}{a}\right)

and in both cases we calculate the same semi-major axis a = 1.00000. (You may find this value a little suspicious - but it's just a result of me setting up a test problem with simple parameters.)

Next, we solve the quadratic:

A-B\,X = 2\,a\,(1 - X^2)

The two roots of this are X_1 = 0.00999118 and X_2 = 0.9796683.

With X = X_1 = 0.00999118, the radial and transverse components at \mathbf{R}_i are:

\left(v_{i,r},v_{i,\theta}\right) = \sqrt{\frac{2\,\mu}{A-B\,X}}\,\left(\sqrt{\frac{r_f}{r_i}}\,\cos\left(\frac{\theta}{2}\right)- X,\sqrt{\frac{r_f}{r_i}}\,\sin\left(\frac{\theta}{2}\right)\right) = \left(0.9702056, 0.1986958\right)

and the radial and transverse components at \mathbf{R}_f are:


\left(v_{f,r},v_{f,\theta}\right) = \sqrt{\frac{2\,\mu}{A-B\,X}}\,\left(X - \sqrt{\frac{r_i}{r_f}}\,\cos\left(\frac{\theta}{2}\right),\sqrt{\frac{r_i}{r_f}}\,\sin\left(\frac{\theta}{2}\right)\right) = \left(-0.9700421, 0.1986627\right)

And with X = X_2 = 0.9796683, the radial and transverse components at \mathbf{R}_i are:

\left(v_{i,r},v_{i,\theta}\right) = \sqrt{\frac{2\,\mu}{A-B\,X}}\,\left(\sqrt{\frac{r_f}{r_i}}\,\cos\left(\frac{\theta}{2}\right)- X,\sqrt{\frac{r_f}{r_i}}\,\sin\left(\frac{\theta}{2}\right)\right) = \left(0.002392613, 0.9903399\right)

and the radial and transverse components at \mathbf{R}_f are:


\left(v_{f,r},v_{f,\theta}\right) = \sqrt{\frac{2\,\mu}{A-B\,X}}\,\left(X - \sqrt{\frac{r_i}{r_f}}\,\cos\left(\frac{\theta}{2}\right),\sqrt{\frac{r_i}{r_f}}\,\sin\left(\frac{\theta}{2}\right)\right) = \left(-0.001577536, 0.9901747\right)

Based on this, if we go through the straightforward exercise of calculating the eccentricity with X=X_1 we find that the orbital eccentricity is 0.9796683; and with X=X_2 we find that the orbital eccentricity is 0.010000. Since I know that the orbit is nearly circular (i.e., so that the eccentricity is close to zero), the appropriate solution to use is the one corresponding to X=X_2.

So, finally, the radial and transverse components of the velocity vector at \mathbf{R}_i are (0.002392613, 0.9903399); and the radial and transverse components of the velocity vector at \mathbf{R}_f are (-0.001577536, 0.9901747).

Last edited by MontBlanc2012; 07-12-2019 at 02:10 PM.
MontBlanc2012 is offline   Reply With Quote
Thanked by:
Old 07-12-2019, 08:49 PM   #7
BrianJ
Addon Developer
Default

Hi,
wow! Thanks so much, I would never have got those later steps. Magic indeed
BrianJ is offline   Reply With Quote
Thanked by:
Old 07-13-2019, 03:38 AM   #8
MontBlanc2012
Orbinaut
Default

Actually, it occurs to me that the solution to ncc1701d's original problem is actually quite trivial - given that we are given the time of flight on each of the arc segments (120 seconds each). This means that the information about speeds is superfluous and one just needs to use a standard Lambert Solver to solve the problem.

However, we can also use the maths set out above to do the same job of solving the Lambert Problem, - so, in outline, here's the solution:

In SI units, the first and second Jupiter-centric positions of Amalthea are:

\begin{aligned}
R_i &= \left(-178932619.,-28063045.,-17448755.\right)\\
R_f &= \left(-178337419.,-30878143.,-18771071.\right)
\end{aligned}

We also know that for Jupiter, the accepted value (in SI units) of the gravitational parameter is \mu = 1.26687\times 10^{17}. And we know the time of flight is 120 seconds.

As before calculate:

\begin{aligned}
r_i &= \sqrt{\mathbf{R}_i.\mathbf{R}_i} = 181958444.95\\
r_f &= \sqrt{\mathbf{R}_f.\mathbf{R}_f} = 181961665.85\\
\theta &=\cos^{-1}\left(\frac{\mathbf{R}_i.\mathbf{R}_f}{r_i\,r_f}\right) = 0.01740310807\\
A &= r_i + r_f = 363920110.79\\
B &= 2\,\sqrt{r_i\,r_f}\,\cos\left(\frac{\theta}{2}\right) = 363906333.39
\end{aligned}

Now set X and Y as functions of an as yet unknown semi-major axis a:

\begin{aligned}
X_+(a) &= \frac{B+\sqrt{16\,a^2 - 8\,a\,A + B^2}}{4\,a} \\
X_-(a) &= \frac{B-\sqrt{16\,a^2 - 8\,a\,A + B^2}}{4\,a}
\end{aligned}

where X_+ and X_- are the two roots of the equation A - B\,X =2\,a\,(1-X^2).

As it turns out, for elliptical orbits, the transfer time from the initial to final points can be calculated as:

\Delta T =120 = 2\,\sqrt{\frac{a^3}{\mu}}\,\left(\cos^{-1}X_+(a) + X_-(a)\,\sqrt{1-X_+(a)^2}\right)

So, now we have an expression that is a function of a only. And since we know that the orbit of Amalthea is almost circular, we use a standard Newton root-finding algorithm with an initial guess of 180,000 km. From this, we deduce that:

a = 181997386.85

and hence

X_+ = 0.99996215721

and

X_- = -0.0002051286

Then, the radial and transverse components of the velocity vector at the \mathbf{R}_i are given by:

\left(v_{i,r},v_{i,\theta}\right) = \sqrt{\frac{2\,\mu}{A-B\,X_+}}\,\left(\sqrt{\frac{r_f}{r_i}}\,\cos\left(\frac{\theta}{2}\right)- X_+,\sqrt{\frac{r_f}{r_i}}\,\sin\left(\frac{\theta}{2}\right)\right) = (26.79327, 26389.15)

and the radial and transverse components of the velocity vector at the \mathbf{R}_f are given by:

\left(v_{f,r},v_{f,\theta}\right) = \sqrt{\frac{2\,\mu}{A-B\,X_+}}\,\left(X_+ - \sqrt{\frac{r_i}{r_f}}\,\cos\left(\frac{\theta}{2}\right),\sqrt{\frac{r_i}{r_f}}\,\sin\left(\frac{\theta}{2}\right)\right) = (26.88698, 26388.68)

And from this it is straightforward to convert back to the x-y-z coordinates of the original coordinate system.

N.B. In case you were wondering what the orbital eccentricity is, we get this from the expression:

e = \sqrt{X_-^2 + \frac{(r_1-r_2)^2}{A^2 - B^2}\,(1 - X_-^2)} = 0.001037615

i.e., very nearly a circular orbit.

Yay!


---------- Post added at 03:38 AM ---------- Previous post was at 02:15 AM ----------

And finally, just a quick comment on ncc1701d's summation of velocity vector components (V = VX + VY + VZ). Initially, I thought this was a typo and ncc1701d really meant V^2 = VX^2 + VY^2 + VZ^2 - which is, by far and away the more conventional way of 'adding' velocity components.

As it turns, out, I was wrong: ncc1701d actually did mean the straight arithmetic sum of the velocity components. How do I know this? Well, I can calculate the x-y-z representation of the radial and transverse unit vectors at that point. Although I'm not going to bother showing the intermediate calculations, in x-y-z coordinates the radial unit vector is:

\hat{\mathbf{r}} = \left(-0.9833708, -0.1542278, -0.09589418\right)

and the transverse unit vector is:

\hat{\mathbf{t}} = \left(0.1804087, -0.8901994, -0.4183273\right)

So, the velocity vector at the first point in ncc1701d's list can be written as:

\mathbf{V}_i = 26.79327\,\hat{\mathbf{r}} + 26389.15\,\hat{\mathbf{t}} = \left(4734.483, -23495.74, -11041.87\right)

So, at that point, VX = 4734.483; VY = -23495.74; and VZ = -11041.87. And the straight arithmetic sum of those is -29803.12 which is suspiciously close to the value given by ncc1701d of -29804.11.

Last edited by MontBlanc2012; 07-17-2019 at 03:43 AM.
MontBlanc2012 is offline   Reply With Quote
Thanked by:
Old 07-14-2019, 03:26 PM   #9
BrianJ
Addon Developer
Default

Thanks for the further insight.
I was on the wrong track to start with, since I had assumed the OP actually meant V^2 = VX^2 + VY^2 + VZ^2.
Note to self: always read the question carefully!
Cheers,
Brian
BrianJ is offline   Reply With Quote
Thanked by:
Old 07-16-2019, 07:35 AM   #10
ncc1701d
Orbinaut
Default

thank you everyone for all the feedback!
ncc1701d is offline   Reply With Quote
Old 09-25-2019, 06:15 AM   #11
MontBlanc2012
Orbinaut
Default

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, \mathbf{R}_i and \mathbf{R}_f. 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, \mathbf{R}_i = (x_i, y_i, z_i) and \mathbf{R}_f = (x_f, y_f, z_f), and as before we also know the semi-major axis, a of the Keplerian arc; and the gravitational parameter, \mu.

(N.B., in the above posts, the semi-major axis was calculated from a time-of-flight equation.)

Then, calculate the following quantities:

\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}

Here, r_i and r_f are just the distances of the initial and final points from the centre of the gravitating body; A is just the sum of these two distances; and \ell is the length of the chord connecting the initial and final point. The quantity B doesn't have a simple geometric interpretation, but it is just a simple function of A and \ell.

Note, though, the \pmin front of the square root in the expression for B. Here, we have to make a choice about the sign of B 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 \mathbf{R}_i to the final point \mathbf{R}_f, as seen from the gravitating body, the satellite traces out a trajectory across the sky that spans some angle \theta. If that angle is between 0 and \pi (a 'short' arc), we use the positive sign to calculate B. And if the angle traced out is between \pi and 2\,\pi (a 'long' arc), we use the negative sign to calculate B. In other words, by choosing the sign of B 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:

a \ge \frac{1}{4}\,\left(r_i + r_f + \ell\right)

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

(N.B. If the semi-major axis equals the lowest permissible value - i.e., a = (r_i + r_f + \ell)/4 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:

\begin{aligned}
  \phi_i &= \sqrt{(z_i + r_i) / 2} \\
  \phi_f &= \sqrt{(z_f + r_f) / 2} \\
\end{aligned}

\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}

\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}

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:

\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}

(N.B. An arbitrary quaternion \mathcal{Q} can be written in the form: \mathcal{Q} = q_0\, + q_1\,\hat{\mathbf{i}} + q_2\,\hat{\mathbf{j}} + q_3\,\hat {\mathbf{k}} and it's conjugate as \mathcal{Q}^\dagger = q_0\, - q_1\,\hat{\mathbf{i}} - q_2\,\hat{\mathbf{j}} - q_3\,\hat{\mathbf{k}})


Step 5
As in the earlier posts, calculate the two roots, X_\pm, of the quadratic, A - B\,X = 2\,a\,(1 - X^2).
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 a - given our choice of the sign of B.

For parabolic orbits where we would expect a divergence in the semi-major axis a, we can side-step solving for the roots of the quadratic and simply set X=\pm 1. Also note that for elliptical orbits |X|<1; and for hyperbolic orbits |X|>1.


Step 6
For each root, X_\pm, calculate the quaternions \mathcal{V}_i and \mathcal{V}_f:

\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}

The resulting quaternions are 'pure' quaternions. To extract the velocity vector at the initial point, one takes the coefficients of \hat{\mathbf{i}}, \hat{\mathbf{j}} and \hat{\mathbf{k}} of \mathcal{V}_i 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.
MontBlanc2012 is offline   Reply With Quote
Reply

  Orbiter-Forum > Far Side of the Moon > Math & Physics


Thread Tools

Posting Rules
BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts
Forum Jump


All times are GMT. The time now is 05:23 AM.

Quick Links Need Help?


About Us | Rules & Guidelines | TOS Policy | Privacy Policy

Orbiter-Forum is hosted at Orbithangar.com
Powered by vBulletin® Version 3.8.11
Copyright ©2000 - 2019, vBulletin Solutions Inc.
Copyright 2007 - 2017, Orbiter-Forum.com. All rights reserved.