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 09-19-2011, 03:29 PM   #1
HarvesteR
Orbinaut
 
HarvesteR's Avatar
Default Calculate (not derive) Orbital Velocity Vector

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 by HarvesteR; 09-19-2011 at 03:34 PM.
HarvesteR is offline   Reply With Quote
Old 09-19-2011, 06:07 PM   #2
Moach
Crazy dude with a rocket
 
Moach's Avatar
Default

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

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 by Moach; 09-19-2011 at 06:12 PM.
Moach is offline   Reply With Quote
Old 09-19-2011, 09:33 PM   #3
HarvesteR
Orbinaut
 
HarvesteR's Avatar
Default

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
HarvesteR is offline   Reply With Quote
Old 09-20-2011, 05:16 PM   #4
Calsir
Orbinaut
Default

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 by Calsir; 09-20-2011 at 07:17 PM. Reason: added mathworld link
Calsir is offline   Reply With Quote
Old 09-22-2011, 03:51 PM   #5
Miner1
Orbinaut
Default

The radial and along track velocity components are

v_r = \frac{\mu e sin\nu}{h}
v_\theta = \frac{h}{r}

where h = \sqrt[]{\mu a(1-e^2)}

You can then use \Omega, i, and \theta = \omega+\nu to transform these components into global coordinates. Hope this helps.
Miner1 is offline   Reply With Quote
Old 10-18-2011, 03:54 AM   #6
HarvesteR
Orbinaut
 
HarvesteR's Avatar
Default

Quote:
Originally Posted by Calsir View Post
 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?

Quote:
Originally Posted by Miner1 View Post
 The radial and along track velocity components are

v_r = \frac{\mu e sin\nu}{h}
v_\theta = \frac{h}{r}

where h = \sqrt[]{\mu a(1-e^2)}

You can then use \Omega, i, and \theta = \omega+\nu to transform these components into global coordinates. Hope this helps.
I'm going to try this too! thanks!

Cheers
HarvesteR is offline   Reply With Quote
Old 10-18-2011, 05:08 AM   #7
n72.75
Donator
 
n72.75's Avatar


Default

Read this: http://spiff.rit.edu/richmond/nbody/...ungeKutta4.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.
n72.75 is offline   Reply With Quote
Old 10-25-2011, 07:24 PM   #8
HarvesteR
Orbinaut
 
HarvesteR's Avatar
Default

Quote:
Originally Posted by n72.75 View Post
 Read this: http://spiff.rit.edu/richmond/nbody/...ungeKutta4.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
HarvesteR is offline   Reply With Quote
Old 10-26-2011, 04:50 PM   #9
HopDavid
Hop David
 
HopDavid's Avatar
Default

Quote:
Originally Posted by HarvesteR View Post
 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



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

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



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



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.



(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.
HopDavid is offline   Reply With Quote
Old 10-30-2011, 03:50 PM   #10
HarvesteR
Orbinaut
 
HarvesteR's Avatar
Default

Thanks man! This look like it will work perfectly.

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 by HarvesteR; 10-30-2011 at 03:57 PM.
HarvesteR is offline   Reply With Quote
Old 10-31-2011, 03:29 PM   #11
HopDavid
Hop David
 
HopDavid's Avatar
Default

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.
HopDavid is offline   Reply With Quote
Old 11-01-2011, 06:22 PM   #12
HarvesteR
Orbinaut
 
HarvesteR's Avatar
Default

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!

Cheers
HarvesteR is offline   Reply With Quote
Old 11-03-2011, 07:23 PM   #13
Miner1
Orbinaut
Default

Or you could do this. Start with the equation for a conic section, i.e.

r = \frac{p}{1+e\cos{\nu}}

then take the derivative with respect to time to get the radial component of velocity as

v_r = \frac{pe\omega\sin{\nu}}{(1+e\cos{\nu})^2}

Now multiply and divide by p and use the first equation above to get

v_r = \frac{e\omega\sin{\nu}}{p} \frac{p^2}{(1+e\cos{\nu})^2} = \frac{e\sin{\nu}}{p}\omega r^2

Recall that the magnitude of an orbit's specific angular momentum is h = \omega r^2 = \sqrt{\mu p} and you can simplify the radial velocity component to

v_r = \frac{e\sin{\nu}}{h^2/\mu}h = \frac{\mu e\sin{\nu}}{h}

For the along track component, it is simply

v_\theta = \omega r = \frac{\omega r^2}{r} = \frac{h}{r}

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

v_x = v_r \cos{\nu}-v_\theta \sin{\nu}
v_y = v_r \sin{\nu}+v_\theta \cos{\nu}

To find the angle \alpha note from the previously posted diagrams that

\gamma = \frac{\pi}{2}-\frac{\pi-\alpha}{2} = \frac{\alpha}{2}

Therefore \alpha = 2\gamma where \gamma is the flight path angle and it is related to the radial and along track velocity components by

\tan{\gamma} = \frac{v_r}{v_\theta}

Using the above expressions for v_r and v_\theta one can readily show

\tan{\gamma} = \frac{\mu e\sin{\nu}}{h} \frac{r}{h} = \frac{er\sin{\nu}}{h^2/\mu} = \frac{er\sin{\nu}}{p}

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 \nu. This means that the velocity can be calculated quite efficiently in computer code.
Miner1 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 04:39 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.6
Copyright ©2000 - 2017, Jelsoft Enterprises Ltd.
Copyright 2007 - 2017, Orbiter-Forum.com. All rights reserved.