Project Math problems with RPOS and RVEL

[Albinon]

I am developing an add-on that creates a Galaxy of randomly created stars and planets.

Part of this add-on will take an existing scenario and modify it for one of the systems I have randomly created.

It works up to the point of calculating the RPOS and RVEL parameters.

I do not have the math skills to do this myself so I have "borrowed" the math from other add-ons.

This add-on is written in Visual Studio 2010 Visual Basic (which handles forms a lot easier than VC does imho). It is not meant to interface directly with orbiter.

Here is the code that is causing me a problem:

Dim r AsDouble = Val(BodyRadius.Text) + Val(OrbitSMA.Text)
Dim lng AsDouble = Val(Long_In.Text)
Dim clong AsDouble = System.Math.Cos(lng)
Dim slng AsDouble = System.Math.Sin(lng)
Dim lat AsDouble = Val(Lat_In.Text)
Dim clat AsDouble = System.Math.Cos(lat)
Dim slat AsDouble = System.Math.Sin(lat)
Dim M AsDouble = Val(BodyMass.Text)
' position in cartesian coordinates
Dim x AsDouble = r * System.Math.Cos(lat) * System.Math.Cos(lng)
Dim z AsDouble = r * System.Math.Cos(lat) * System.Math.Sin(lng)
Dim y AsDouble = r * System.Math.Sin(lat)
' Base velocity of orbit
Dim G AsDouble = 0.0000000000667 ' Gravitational Constant
Dim v AsDouble = System.Math.Sqrt(G * M / r)
Dim inc AsDouble = Val(OrbitInc.Text) / (180 / System.Math.PI) ' to radians
' assuming a standard (non-retrograde) orbit
Dim vx AsDouble = -v * System.Math.Cos(inc)
Dim vy AsDouble = v * System.Math.Sin(inc)
Dim vz AsDouble = 0.0
' derivatives in cartesian coordinates
RPOS_Out.Text = x & " " & y & " " & z
RVEL_Out.Text = vx & " " & vy & " " & vz

Can anyone give me some help in getting the correct values for RVEL and RPOS?

(I set vz to 0.0 as I could not find the math to properly calulate it)

 

[ADSWNJ]

Did you want your users to input long and lat in radians rather than degrees? If not, you need to convert by PI/180 like you did later.

What do you want RPOS and RVEL to have in them?

 

[SolarLiner]

Yeah, it depends on what unit you want the user to type the elements.

But on a side note, you just solved my problems for my AOSP project ! And all that syntax coloring, it must have took so much time !

 

[Albinon]

Did you want your users to input long and lat in radians rather than degrees?

Oops [face-slap] Thank you. User enters degrees.

What do you want RPOS and RVEL to have in them?

The actual parameters used by the scenario parser to position the ship in orbit over the planet when it reads the scenario. (for RPOS, the X, Y and Z positions. For RVEL, the X, Y and Z velocities.)

---------- Post added at 08:06 AM ---------- Previous post was at 08:02 AM ----------

And all that syntax coloring, it must have took so much time !

Naw... just cut and paste from VB 2010

 

[ADSWNJ]

So ... just to see if I have this right ...

you have an input inclination val(OrbitInc.Text), and altitude val(OrbitSMA.Text) (??? SMA would have the planet radius in it already), presumably assuming spherical planet and a 0.000 eccentricity ...

and you want to calculate the XYZ position and XYZ velocity in planet coords for a user-input LAT/LONG.

You are missing the Longitude of Ascending Node (LAN), to determine the orientation of the orbit - i.e. INC by itself does not lock down a circular orbit, without a LAN.

 

[SolarLiner]

I mean, if you copy paste the code from VS to the forum, you have plain text. Here it is beautifully colored !

 

[Albinon]

(??? SMA would have the planet radius in it already)

This is the variable I named. It has no connection to orbiter. It is the height in meters above the planet's surface in this case.

and you want to calculate the XYZ position and XYZ velocity in planet coords for a user-input LAT/LONG.

Basically Correct.

You are missing the Longitude of Ascending Node (LAN), to determine the orientation of the orbit - i.e. INC by itself does not lock down a circular orbit, without a LAN.

As I said, I don't do math very well. I assimulate other people's code into my projects when I do math problems. I am familiar with basic trig (distance formulae and 3D angles) but when it comes to orbital mechanics, I'm pretty much out of my depth. I know the terminology and have the basic idea, I just don't know how to apply it.

Can you help me?

---------- Post added at 09:54 AM ---------- Previous post was at 09:52 AM ----------

I mean, if you copy paste the code from VS to the forum, you have plain text. Here it is beautifully colored !

You are correct when using VC++. I am using VB. Different text ediing.

I just copied and pasted straingt from the code window.

 

[SolarLiner]

I just copied and pasted straingt from the code window.

I use VB with VS2010 (as well as VC#) but if you say so...

 

[Albinon]

but if you say so...

I just don't want credit for doing something that I did not do.

I do appreciate you comments, though.

I'm curious, what was your problem with your project that this thread helped you solve?

 

[SolarLiner]

I'm curious, what was your problem with your project that this thread helped you solve?

My AOSP Project, a way to handle Orbiter scenarios into a .NET app: http://www.orbiter-forum.com/showthread.php?t=31101

 

[Albinon]

My AOSP Project, a way to handle Orbiter scenarios into a .NET app: http://www.orbiter-forum.com/showthread.php?t=31101

Here is an updated code packet with some derivatives that I did not include in my original post.

Dim r AsDouble = Val(BodyRadius.Text) + Val(OrbitSMA.Text) ' Planet radius plus height of orbit above planet surface in meters
Dim lng AsDouble = Val(Long_In.Text) / (180 / System.Math.PI) ' to radians
Dim clong AsDouble = System.Math.Cos(lng)
Dim slng AsDouble = System.Math.Sin(lng)
Dim lat AsDouble = Val(Lat_In.Text) / (180 / System.Math.PI) ' to radians
Dim clat AsDouble = System.Math.Cos(lat)
Dim slat AsDouble = System.Math.Sin(lat)
Dim M AsDouble = Val(BodyMass.Text)
' position in cartesian coordinates
Dim x AsDouble = r * System.Math.Cos(lat) * System.Math.Cos(lng)
Dim z AsDouble = r * System.Math.Cos(lat) * System.Math.Sin(lng)
Dim y AsDouble = r * System.Math.Sin(lat)
' Base velocity of orbit
Dim G AsDouble = 0.0000000000667 ' Gravitational Constant
Dim v AsDouble = System.Math.Sqrt(G * M / r)
Dim inc AsDouble = Val(OrbitInc.Text) / (180 / System.Math.PI) ' to radians
' assuming a standard (non-retrograde) orbit
Dim vx AsDouble = v * System.Math.Cos(inc)
Dim vy AsDouble = v * System.Math.Sin(inc)

I still do not know how to calculate vz properly so I set it to 0.

Dim vz AsDouble = 0.0
' derivatives in cartesian coordinates or Delta-V
Dim dxdt AsDouble = vx * clat * clong - vy * r * clat * slng - vz * r * slat * clong
Dim dzdt AsDouble = vx * clat * slng + vy * r * clat * clong - vz * r * slat * slng
Dim dydt AsDouble = vx * slat + vz * r * clat
' RPOS and RVEL values
RPOS_Out.Text = x & " " & y & " " & z
RVEL_Out.Text = vx & " " & vy & " " & vz

 

[ADSWNJ]

As I said, I don't do math very well. I assimulate other people's code into my projects when I do math problems. I am familiar with basic trig (distance formulae and 3D angles) but when it comes to orbital mechanics, I'm pretty much out of my depth. I know the terminology and have the basic idea, I just don't know how to apply it.

Sure - happy to help, and when I get stuck, I'm sure some of the math gods on this forum can step in! I'll number my thinking for ease of reference...

1. - let's define our Cartesian coordinate system as "ECEF" (see Wiki). It's a right-hand system, with the X-axis going from the planet origin to 0 lat / 0 long (Prime Meridian intersecting Equator). Y-axis is 0 lat / 90E. Z-axis is the North Pole.

2. - consider our ship, mass "m" at planet radius "r", at tangential velocity "v", in circular orbit, with a spherical planet of mass "M". To be in this orbit, the gravitational force (defined as G . M. m / r^2, where G is the gravitational constant) has to provide exactly the centripetal force ( m . v^2 / r) to pull it round the orbit.

3. - from 2 ... velocity = square root of (G.M / r)

4. - let's put this ship exactly in a zero inc equatorial orbit, exactly over the prime meridian. At this instant, the position in ECEF-space will be (r, 0, 0) and the velocity vector will be (0,v,0).

5. - we can now describe the position at longitude A (A = 0 to but not including 360) as ... (r.cos(A),r.sin(A),0), and the velocity vector as(-v.sin(A), -v.cos(A),0). OK so far I hope?!

6. - next let's bring in the inclination. Consider inclination "B" at LAN of 0. I.e. the orbit crosses the equatorial plane from south to north exactly at the prime meridian. Therefore the latitude relative to longitude A is B.sin(A) ... i.e. you reach max latitude at A=90, back to the equator at A=180, then min latitude at A=270 and back round to the LAN at A=0.

7. - we need to rotate our equatorial orbit by angle B around the x-axis. Look up rotation matrices in Wiki, and you'll find that the X-axis rotation "Mx" is:

{ { 1 0 0 }
{ 0 cosB -sinB}
{0 sinB cosB} }.

8. Matrix multiplying Mx . Pos gives (x, y.cosB-z.sinB, y.sinB+z.cosB). You can apply this rotation to the velocity vector as well. (Think of v simply as a point in space 1 second from your current position ... although it's not quite right as the vector will be dragged round a bit, but in a typical 5000 second orbit, it's close enough!).

9. We need one more rotation ... this time around the z-axis (polar axis). This moves the LAN to the desired longitude. Let's say the LAN is "C" degrees. The z-axis rotation matrix "Mz" is

{ { cosC -sinC 0}
{ sinC cosC 0}
{ 0 0 1} }.

10. So we want to put this all together ... we started with a unit point at (r,0,0). We applied the A rotation around the Z-axis to find the position around the equatorial orbit. We next need to rotate by MINUS C around Z to bring the LAN back to 0 degrees, then rotate by B around the X axis, then by PLUS C to put back the LAN where it should B.


The interesting thing is that your user LAT angle can be used as a check. I originally thought that you could determine the LAN from the LAT angle, but it's only true where the LAT is exactly the +Inc or the -Inc. At all other values, it could be one of two points on the orbit - ascending and descending through your user LAT.

I hope this is useful for starters. See how you get on with this.

 

[Hielor]

I use VB with VS2010 (as well as VC#) but if you say so...

What browser do you use? Some browsers may not support pasting formatted text.

 

[SolarLiner]

I am using Chrome and Firefox.

 

[Hielor]

There you go--IE supports it, I guess those two don't.

 

[ADSWNJ]

one other thing to factor in ... the planet is rotating! I was thinking about geostationary orbits. in ECEF coords, they are static, but obviously rotating in phase with the Earth. You need to get the planet's rotation period and turn that into a base velocity offset.

 

[Albinon]

I hope this is useful for starters. See how you get on with this.

Ungh.... I barely passed 1st semester calculus after 3 attempts.

Matrix math in wiki is like learning greek with out a rosetta stone. I'll trip over it a while and see what comes out other than Nan or -1.#AnF

 

[ADSWNJ]

Yup - I know how you feel. Initially, these rotation matrices with both direct multiplication (mul()) and transpose multiplication (tmul()) definitely seems like gobbledegook, but unfortunately there's no way round it if you want to convert between coordinate systems and do the math that you asked for in your OP.

Start with more simple 2D coordinate systems and see if you can work through a 2x2 rotation matrix times a 2 x 1 vector to give a 2x1 vector output. Try it with 90 degree, then 45 degree, then 30 degree inputs to make sure the results match what you are visualizing. Then step into 3D, and do it just in the X-Z plane, then just in the Y-Z plane, then try all three.

Let me know if you need help.

 

[Albinon]

one other thing to factor in ... the planet is rotating! I was thinking about geostationary orbits. in ECEF coords, they are static, but obviously rotating in phase with the Earth. You need to get the planet's rotation period and turn that into a base velocity offset.

Eeps!... I for got that one. It was in my initial research, just didn't make it into the final code.

 

[Albinon]

Let me know if you need help.

I got a bit further. The orbit is eccentric so the RVEL VECTOR is still off, but I got the RPOS by barrowing from Wiki.

Dim r AsDouble = Val(BodyRadius.Text) + Val(OrbitSMA.Text)
Dim lng AsDouble = Val(Long_In.Text) / (180 / System.Math.PI) ' to radians
Dim clong AsDouble = System.Math.Cos(lng)
Dim slng AsDouble = System.Math.Sin(lng)
Dim lat AsDouble = Val(Lat_In.Text) / (180 / System.Math.PI) ' to radians
Dim clat AsDouble = System.Math.Cos(lat)
Dim slat AsDouble = System.Math.Sin(lat)
Dim M AsDouble = Val(BodyMass.Text)
Dim inc AsDouble = Val(OrbitInc.Text) / (180 / System.Math.PI) ' to radians
Dim LAN AsDouble = Val(OrbitLAN.Text) / (180 / System.Math.PI) ' to radians
Dim Peri AsDouble = Val(OrbitPeri.Text) / (180 / System.Math.PI) ' to radians

' position in cartesian coordinates
'Dim x As Double = r * System.Math.Cos(lat) * System.Math.Cos(lng)
'Dim z As Double = r * System.Math.Cos(lat) * System.Math.Sin(lng)
'Dim y As Double = r * System.Math.Sin(lat)
Dim x1 AsDouble = r * (System.Math.Cos(LAN) * System.Math.Cos(Peri) - System.Math.Sin(LAN) * System.Math.Cos(inc) * System.Math.Sin(Peri))
Dim x2 AsDouble = r * (System.Math.Sin(LAN) * System.Math.Cos(Peri) + System.Math.Cos(LAN) * System.Math.Cos(inc) * System.Math.Sin(Peri))
Dim x3 AsDouble = r * (System.Math.Sin(inc) * System.Math.Sin(Peri))
Dim y1 AsDouble = r * (-System.Math.Cos(LAN) * System.Math.Sin(Peri) - System.Math.Sin(LAN) * System.Math.Cos(inc) * System.Math.Cos(Peri))
Dim y2 AsDouble = r * (-System.Math.Sin(LAN) * System.Math.Sin(Peri) + System.Math.Cos(LAN) * System.Math.Cos(inc) * System.Math.Cos(Peri))
Dim y3 AsDouble = r * (System.Math.Sin(inc) * System.Math.Cos(Peri))
Dim z1 AsDouble = r * (System.Math.Sin(inc) * System.Math.Sin(LAN))
Dim z2 AsDouble = r * (-System.Math.Sin(inc) * System.Math.Cos(LAN))
Dim z3 AsDouble = r * (System.Math.Cos(inc))

Dim x AsDouble = x1 + x2 + x3
Dim y AsDouble = y1 + y2 + y3
Dim z AsDouble = z1 + z2 + z3

' Base velocity of orbit
Dim G AsDouble = 0.0000000000667 ' Gravitational Constant
Dim v AsDouble = System.Math.Sqrt(G * M / r)

' assuming a standard (non-retrograde) orbit
Dim vx AsDouble = v * (System.Math.Sin * System.Math.Cos(x))
Dim vy AsDouble = v * System.Math.Sin(x)

'Dim vz As Double = 0.0
Dim vz AsDouble = v * (-System.Math.Cos * System.Math.Cos(x))
' derivatives in cartesian coordinates
Dim dxdt AsDouble = vx * clat * clong - vy * r * clat * slng - vz * r * slat * clong
Dim dzdt AsDouble = vx * clat * slng + vy * r * clat * clong - vz * r * slat * slng
Dim dydt AsDouble = vx * slat + vz * r * clat

' RPOS and RVEL values
RPOS_Out.Text = x & " " & y & " " & z
RVEL_Out.Text = vx & " " & vy & " " & vz
　
Can you give me a hand with the RVEL Vector?

---------- Post added at 08:25 PM ---------- Previous post was at 08:23 PM ----------

Still haven't done Sidereal Rotation calc yet, I wanted to get one thing at a time so I know where the code errors are.