Calculate (not derive) Orbital Velocity Vector

HarvesteR

Member
Joined
Apr 22, 2008
Messages
386
Reaction score
15
Points
18
Hi,

So, I was wondering here, is there a way to calculate an orbital velocity vector without taking the orbital position's derivative?

I searched wikipedia a lot, but haven't found anything other than the usual derivative approach.

I know the orbital velocity magnitude can be calculated using the vis-viva equation, but I haven't been able to find anything as elegant for the velocity direction.

The reason I'd rather not use a derivative is noise and floating point inaccuracies... If there were a way to get the velocity vector like that, my orbital code would be a lot more reliable.

All other orbital parameters are already calculated here. I'm propagating the orbits simply by stepping through time, so the velocity here is pretty much just a result, not an input. I need it to be as precise as possible, though, because later this calculated velocity will be used to throw the propagated satellite back into a simulated orbit.

I'd greatly appreciate any help :)

Thanks in advance!

Cheers
 
Last edited:

Moach

Crazy dude with a rocket
Addon Developer
Joined
Aug 6, 2008
Messages
1,581
Reaction score
62
Points
63
Location
Vancouver, BC
well, if you can get the parameters to propagate your position over the orbit, then it should be simple enough to calculate two points close enough to each other so you can use the normalized delta between them as the velocity direction, no?

yeah, it's dirty... but should work :rolleyes:

and it's not really the same as a derivative, i suppose... well, we're not really propagating over time, but you're still using a time-lapse in order to get something akin to a derivative, yes... but instead of propagating it, you simply normalize and use it as a direction


better yet, re-reading your post....

calculate the precise delta between the "launch back to real orbit" position and the position at current + 1 times the current simulation delta-time... THAT's your velocity - no need to do anything else
 
Last edited:

HarvesteR

Member
Joined
Apr 22, 2008
Messages
386
Reaction score
15
Points
18
Yeah, that's what I'm doing now, in fact. My question was whether there wasn't a way to calculate the orbital velocity direction without having to take a delta.

This was mostly because of problems with numerical precision and things like that. I was able to kind of get around the problem, by replacing unity's float-based vectors and quaternions with my own, double-precision versions.

They're not as good performance-wise (since Unity's data types are embedded in the engine code), but they're much, much more accurate.

All the orbital math now is being done in double precision, so these problems are not likely to reappear soon.

Cheers
 

Calsir

New member
Joined
Jan 8, 2009
Messages
71
Reaction score
0
Points
0
If you know the orbital elements, then you know the shape of your orbit and where you are upon it.

The direction of your velocity is unit vector which is tangent to your orbit in your position. Such a tangent vector can be calculated analitically for Keplerian orbits.

Edit: look at the last paragraph here: http://mathworld.wolfram.com/Ellipse.html .
 
Last edited:

Miner1

New member
Joined
Jun 9, 2008
Messages
25
Reaction score
0
Points
0
The radial and along track velocity components are

[math]v_r = \frac{\mu e sin\nu}{h}[/math][math]v_\theta = \frac{h}{r}[/math]
where [math]h = \sqrt[]{\mu a(1-e^2)}[/math]
You can then use [math]\Omega[/math], [math]i[/math], and [math]\theta = \omega+\nu[/math] to transform these components into global coordinates. Hope this helps.
 

HarvesteR

Member
Joined
Apr 22, 2008
Messages
386
Reaction score
15
Points
18
If you know the orbital elements, then you know the shape of your orbit and where you are upon it.

The direction of your velocity is unit vector which is tangent to your orbit in your position. Such a tangent vector can be calculated analitically for Keplerian orbits.

Edit: look at the last paragraph here: http://mathworld.wolfram.com/Ellipse.html .

Those equations do look useful. I'm just not sure what that 't' term is supposed to be... Is it true anomaly?

The radial and along track velocity components are

[math]v_r = \frac{\mu e sin\nu}{h}[/math][math]v_\theta = \frac{h}{r}[/math]
where [math]h = \sqrt[]{\mu a(1-e^2)}[/math]
You can then use [math]\Omega[/math], [math]i[/math], and [math]\theta = \omega+\nu[/math] to transform these components into global coordinates. Hope this helps.

I'm going to try this too! thanks!

Cheers
 

n72.75

Move slow and try not to break too much.
Orbiter Contributor
Addon Developer
Tutorial Publisher
Donator
Joined
Mar 21, 2008
Messages
2,696
Reaction score
1,353
Points
128
Location
Saco, ME
Website
mwhume.space
Preferred Pronouns
he/him
Read this: http://spiff.rit.edu/richmond/nbody/OrbitRungeKutta4.pdf

It's the way orbiter does it and shouldn't take you long to program.

If you're worried about numerical accuracy, use a shorter time step and a higher order numerical approximation.

One way to test for accuracy is to run multiple cases of a two-body point-pass system. Record data for the semi-major axis of the orbit and take the standard deviation of it for the whole case. In cases with longer time steps you will eventually see the standard-deviation spike rapidly, if you stay under this limit you're all set.
 

HarvesteR

Member
Joined
Apr 22, 2008
Messages
386
Reaction score
15
Points
18
Read this: http://spiff.rit.edu/richmond/nbody/OrbitRungeKutta4.pdf

It's the way orbiter does it and shouldn't take you long to program.

If you're worried about numerical accuracy, use a shorter time step and a higher order numerical approximation.

One way to test for accuracy is to run multiple cases of a two-body point-pass system. Record data for the semi-major axis of the orbit and take the standard deviation of it for the whole case. In cases with longer time steps you will eventually see the standard-deviation spike rapidly, if you stay under this limit you're all set.


Thanks, but this setup needs to work as a precalculated 2-body system. It's used to move planets and moons around, and a stable, and most importantly, deterministic system is more important than physical accuracy.

I'm aware that I could run the physical integration at startup to put the solar system in the required time frame, but it would be much more computationally expensive (and possibly unstable) than just setting a universal time variable and letting the orbit propagator calculate a position.

Cheers
 

HopDavid

Hop David
Joined
Feb 1, 2009
Messages
63
Reaction score
0
Points
0
Location
Ajo
Website
clowder.net
Hi,

So, I was wondering here, is there a way to calculate an orbital velocity vector without taking the orbital position's derivative?

I searched wikipedia a lot, but haven't found anything other than the usual derivative approach.

I know the orbital velocity magnitude can be calculated using the vis-viva equation, but I haven't been able to find anything as elegant for the velocity direction.

The reason I'd rather not use a derivative is noise and floating point inaccuracies... If there were a way to get the velocity vector like that, my orbital code would be a lot more reliable.

All other orbital parameters are already calculated here. I'm propagating the orbits simply by stepping through time, so the velocity here is pretty much just a result, not an input. I need it to be as precise as possible, though, because later this calculated velocity will be used to throw the propagated satellite back into a simulated orbit.

I'd greatly appreciate any help :)

Thanks in advance!

Cheers

Start with a coordinate system that puts perihelion on x axis and the sun at the origin

VelVect1.jpg


Position vector is |r| * ( cos(f), sin(f) )

Now draw a line segment to 2nd focus. Call angle F2 P F1 alpha

VelVect2.jpg


The velocity vector lies on the tangent line at point P. Using a well known reflective property of the ellipse and the fact a straight line subtends pi radians, it can be deduced angle between r and v is
(pi - alpha) / 2.

The angle of v is f + ((pi - alpha) / 2)

But what is alpha?

A well known property of the ellipse is
F2P + PF1 = 2a
since PF1 = r
F2P = 2a - r
And of course, the distance between foci is 2ea

VelVect3.jpg


Knowing the 3 sides of the above triangle allows us find alpha using the law of cosines. However plugging the above into the law of cosines gives an unruly mess that gives me a headache just looking at it.

Setting k = r/a, we can get a similar triangle that's much nicer.

VelVect4.jpg


(2 - k)2 + k2 - 2k(2-k)cos(alpha) = (2e)2

alpha = acos(((2-2e2)/(k(2-k)))-1)

Hopefully that gives you enough to find direction of velocity vector. You mentioned you know the Vis-Viva equation which gives you magnitude.

I seem to remember seeing a rotation matrix using the elements OM w and i where you could put this back in the standard coordinate system. I could have sworn it was in Fundamentals of Astrodynamics by Bate Mueller and White. Will look for it when I have time.
 

HarvesteR

Member
Joined
Apr 22, 2008
Messages
386
Reaction score
15
Points
18
Thanks man! This look like it will work perfectly. :thumbup:

The angle is really all I need. With that I can use a quaternion to rotate the position vector, and the magnitude I can get from the Vis-Viva equation:
Code:
Vector3d vel = Quaternion.AngleAxis(alpha + ((PI - alpha) / 2.0), h) * -pos.normalized * velMag; // h is my orbit plane normal vector
Thanks again!

Cheers
 
Last edited:

HopDavid

Hop David
Joined
Feb 1, 2009
Messages
63
Reaction score
0
Points
0
Location
Ajo
Website
clowder.net
My pleasure! I've always enjoyed playing with conic sections.

It's a good idea to check my work. There were 4 or 5 steps between the law of cosines and the last equation and I'm prone to making stupid copying errors.
 

HarvesteR

Member
Joined
Apr 22, 2008
Messages
386
Reaction score
15
Points
18
Just to follow up:

I used your cosine law scheme here, and got a much, much more precise velocity vector here. I had to make a few adjustments (if trueAnomaly > PI I had to invert the angle, because Acos returns a value from -PI to PI), but it works really well now.

To finally solve the inaccuracy when launching the propagated object into a physically integrated trajectory, I also had to get the position at the next frame. I found that there was a significant shift in velocity because the orbit propagation was giving a one-frame-late velocity vector.

well, I don't think it can get any more precise than this now, short of improving the orbit tracking precision, which is kinda hard since the physics simulation works at a lower precision level, and there's really not much I can do there.

Again, many thanks for all that! :thumbup:

Cheers
 

Miner1

New member
Joined
Jun 9, 2008
Messages
25
Reaction score
0
Points
0
Or you could do this. Start with the equation for a conic section, i.e.

[math]r = \frac{p}{1+e\cos{\nu}}[/math]
then take the derivative with respect to time to get the radial component of velocity as

[math]v_r = \frac{pe\omega\sin{\nu}}{(1+e\cos{\nu})^2}[/math]
Now multiply and divide by [math]p[/math] and use the first equation above to get

[math]v_r = \frac{e\omega\sin{\nu}}{p} \frac{p^2}{(1+e\cos{\nu})^2} = \frac{e\sin{\nu}}{p}\omega r^2[/math]
Recall that the magnitude of an orbit's specific angular momentum is [math]h = \omega r^2 = \sqrt{\mu p}[/math] and you can simplify the radial velocity component to

[math]v_r = \frac{e\sin{\nu}}{h^2/\mu}h = \frac{\mu e\sin{\nu}}{h}[/math]
For the along track component, it is simply

[math]v_\theta = \omega r = \frac{\omega r^2}{r} = \frac{h}{r}[/math]
These are the equations from my first post. Now to rotate these components into the frame suggested by HopDavid, the following relations are readily found

[math]v_x = v_r \cos{\nu}-v_\theta \sin{\nu}[/math][math]v_y = v_r \sin{\nu}+v_\theta \cos{\nu}[/math]
To find the angle [math]\alpha[/math] note from the previously posted diagrams that

[math]\gamma = \frac{\pi}{2}-\frac{\pi-\alpha}{2} = \frac{\alpha}{2}[/math]
Therefore [math]\alpha = 2\gamma[/math] where [math]\gamma[/math] is the flight path angle and it is related to the radial and along track velocity components by

[math]\tan{\gamma} = \frac{v_r}{v_\theta}[/math]
Using the above expressions for [math]v_r[/math] and [math]v_\theta[/math] one can readily show

[math]\tan{\gamma} = \frac{\mu e\sin{\nu}}{h} \frac{r}{h} = \frac{er\sin{\nu}}{h^2/\mu} = \frac{er\sin{\nu}}{p}[/math]
The benefit to using these expressions is that they involve several constants, utilize multiplication with a single division, and the only angle you need to take the sine and cosine of is [math]\nu[/math]. This means that the velocity can be calculated quite efficiently in computer code.
 

se5a

New member
Joined
Jul 1, 2019
Messages
1
Reaction score
0
Points
0
So I've been trying to do this exact same thing, and been going around in circles on this for months.
HopDavid's post was very helpful, however I'm stuck.
I have a 2d velocity vector that is on the orbital plane, which is fine if my orbits are 2d and not retrograde. But I need to turn that into a 3d vector ralitive to the parent body/barycenter, and take inclination into account.
I think this is what HopDavid was referring to at the end of the post when he mentions a matrix rotation, and what HavesterR is doing with the quaternion, but I don't quite understand what he's doing there. Can anyone enlighten me?
 

MontBlanc2012

Active member
Joined
Aug 13, 2017
Messages
138
Reaction score
26
Points
28
I have a 2d velocity vector that is on the orbital plane, which is fine if my orbits are 2d and not retrograde. But I need to turn that into a 3d vector relative to the parent body/barycenter, and take inclination into account.
I think this is what HopDavid was referring to at the end of the post when he mentions a matrix rotation, and what HavesterR is doing with the quaternion, but I don't quite understand what he's doing there. Can anyone enlighten me?

To convert the 2d vector to a 3d vector (relative to the coordinate system of a parent body), one needs 'rotate' to the x-y-z reference frame of the parent body from the natural, perifocal reference frame of the elliptical/hyperbolic orbit. This rotation is defined in terms of three angles: the longitude of the ascending node ([MATH]\Omega[/MATH]); the orbital inclination ([MATH]\iota[/MATH]); and the argument of periapsis ([MATH]\omega[/MATH]). These are just the three standard orbital elements that define the orientation of the orbit in space. They are also the Euler angles that define a rotation from one reference frame (the perifocal) to another.

The perifocal reference frame is defined by three unit vectors (which we can call [MATH]\hat{p}[/MATH], [MATH]\hat{q}[/MATH] and [MATH]\hat{w}[/MATH]) each at right angles to each other - in exactly the same way as the three unit vectors [MATH]\hat{x}[/MATH], [MATH]\hat{y}[/MATH] and [MATH]\hat{z}[/MATH] are right angles to each other in the reference frame of the parent body. Here, though, the [MATH]\hat{p}[/MATH] unit vector points from the origin (i.e. the position of the gravitating body) towards the orbital periapsis; the [MATH]\hat{q}[/MATH] unit vector points at right-angles to that vector but in the direction of the orbital motion and also lying in the orbital plane; and the[MATH]\hat{w}[/MATH] unit vector is perpendicular to the orbital plane (and to the other two vectors). For example, if one knows the orbital radius ([MATH]r[/MATH]) and true anomaly ([MATH]\nu[/MATH]) of a satellite, the perifocal coordinates of the satellite is just [MATH](r\,\cos(\nu),\,r\,\sin(\nu),\,0)[/MATH]). This is why the perifocal reference frame is useful: it leads to simple expressions for the position and velocity vectors of a satellite.

In terms of the x-y-z coordinate system, the x-y-z components of the perifocal unit vectors, [MATH]\hat{p}[/MATH], [MATH]\hat{q}[/MATH] and [MATH]\hat{w}[/MATH], can be written as:


[MATH] \begin{aligned} p_{x}&=\cos(\omega)\,\cos(\Omega)-\cos(\iota)\,\sin(\omega)\,\sin(\Omega)\\p_{y}&=\cos(\iota)\,\sin(\omega)\,\cos(\Omega)+\cos(\omega)\,\sin(\Omega)\\p_{z}&=\sin(\iota)\,\sin(\omega) \end{aligned} [/MATH]

[MATH] \begin{aligned} q_{x}&=-\sin(\omega)\,\cos(\Omega)-\cos(\iota)\,\cos(\omega)\,\sin(\Omega)\\q_{y}&=\cos(\iota)\,\cos(\omega)\,\cos(\Omega)-\sin(\omega)\,\sin(\Omega)\\q_{z}&=\sin(\iota)\,\cos(\omega) \end{aligned} [/MATH]

[MATH] \begin{aligned} w_{x}&=\sin(\iota)\,\sin(\Omega)\\w_{y}&=-\sin(\iota)\,\cos(\Omega)\\w_{z}&=\cos(\iota) \end{aligned} [/MATH]

From this, it's easy to calculate the position (or velocity) vector of the satellite in the 3d x-y-z coordinates of the body-centric coordinate system. Taking the position example above, the x-y-z coordinates are given by:

[MATH]r\,\cos(\nu)\,(p_x, p_y, p_z) + r\,\sin(\nu)\,(q_x, q_y, q_z) + 0\,(w_x, w_y, w_z)[/MATH]

The net effect of these transformations is to perform a rotation of the perifocal reference frame using the three Euler angles, [MATH]\Omega[/MATH], [MATH]\iota[/MATH] and [MATH]\omega[/MATH].

As an example, let's suppose that in the perifocal reference frame, the vector you wish to transform is [MATH](1.0, 2.0, 0.0)[/MATH]. Let's also suppose that the orbital longitude of the ascending node is 180 degrees; the orbital inclination is 25 degrees; and the argument of periapsis is 10 degrees. We calculate the x-y-z coordinates of the unit vectors as:

[MATH]\hat{p} = (-0.984808,-0.157379,0.0733869)[/MATH]
[MATH]\hat{q} = (0.173648,-0.892539,0.416198)[/MATH]
[MATH]\hat{w} = (0,0.422618,0.906308)[/MATH]
so that the x-y-z coordinates of the perifocal vector [MATH](1.0, 2.0, 0.0)[/MATH] is:

[MATH]1.0\times(-0.984808,-0.157379,0.0733869) + 2.0\times(0.173648,-0.892539,0.416198) + 0.0\times (0,0.422618,0.906308)[/MATH]
or just:

[MATH](-0.637511,-1.94246,0.905782)[/MATH]

Does that help?
 
Top