Project Math problems with RPOS and RVEL

[ADSWNJ]

You've totally lost me! No idea what Peri is (as an angle). No comments on what X1,X2,X3 are all about. You seem to be trying to jump directly to XYZ coordinates taking in all the variables, but it would be much easier to read if you did it piece by piece, explaining what each intermediate XYZ coordinate represents.

Also - can you run it and print out intermediate results too?

 

[Albinon]

No idea what Peri is (as an angle).

Argument of Periapsis

No comments on what X1,X2,X3 are all about.

Not entirely sure myself, but it seems to have to work with the 3 Axis.

Also - can you run it and print out intermediate results too?

Tell me which data you want to see and how to format it so it means something to you.

Refer to:
https://en.wikipedia.org/wiki/Orbital_elements

---------- Post added at 11:03 AM ---------- Previous post was at 07:42 AM ----------

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.

:compbash:
OK... This code is a red-herring. After futzing with it, it turns out to be a relationship between two coordinate systems.

Apologies to all.


I still would like to know how to get the correct values for RPOS and RVEL

Anyone have some code that they can toss at me? Thanks!

 

[ADSWNJ]

OK... This code is a red-herring. After futzing with it, it turns out to be a relationship between two coordinate systems.

Apologies to all.


I still would like to know how to get the correct values for RPOS and RVEL

Anyone have some code that they can toss at me? Thanks!

Dammit!! After scratching my head looking at those Euler Angles and trying to determine how it related to your code, I came to the same conclusion - that it's not what you were looking for.

Why not go back to my numbered post and work that through piece by piece?

 

[Albinon]

Why not go back to my numbered post and work that through piece by piece?

I've been mangling the matrix math. I hit a brick wall when I discovered that your solution works only at a MJD of March 21.... (Vernal Equinox)

You need to factor in the changing angle of the earth's axis through out the year.

My "red-Herring" code may be of some use after all because it calculates the adjustment between the current body angle and the required euler reference plane. - ugh.... Which means I can't do the coding in VB 2010 as all the required angle info is in the Oriter simulation. And if I do the coding in C++ using the Orbiter libraries, the oapi I need is already written. Along with the tools i need for Matrix math.

I was using VB so I can do the forms better. I don't like setting up forms in C++. Oh well.

 

[ADSWNJ]

I've been mangling the matrix math. I hit a brick wall when I discovered that your solution works only at a MJD of March 21.... (Vernal Equinox)

Huh? Why's that? I thought the default planetary XYZ coordinates were ECEF, not Vernal Equinox centric.

 

[Albinon]

Huh? Why's that? I thought the default planetary XYZ coordinates were ECEF, not Vernal Equinox centric.

From Orbiter's reference guide on scenario editing:

2.3 Orbital elements​

Use this page to place a vessel into orbit around a celestial body.

First select the object (planet, moon or sun) around which you want your vessel to orbit in the ​

Orbit reference box.​

Next, in the ​

Frame box, select the coordinate frame in which the elements will be expressed. The choices are ecliptic and ref. equator.

The ecliptic frame is defined by the plane of Earth’s orbit, and the direction of the vernal equinox (at epoch J2000.0). This frame is useful e.g. if you want to place a vessel into an orbit for an interplanetary transfer.

The equatorial frame is defined by the equatorial plane of the reference body. It is useful if you want to place the vessel into a specific orbit with respect to the planet surface, e.g. geostationary or polar.

The euler axis in orbiter is based on being centered on the Vernal point. for x, aproxmatley Polaris for y and somewhere between the Taurus and Gemini constellation for z. The axial tilt of earth comes into play as it travels in its orbit around the sun since Earth's x only points at the vernal point on March 21.

(part of my project is to place the vehicle in orbit above a particular point of the referenced celestial body - I had not mentioned it earlier as I was trying to get the basic math down first)

Another issue is that my check results blow up trying to doing the matrix math. I did not do well on it in school.

I am thanking you for your assistence and perseverance with my problem here. I can't do the math itself. And some Math Notation (like matrix) is greek to me so I don't know how to encode it. What I actually need is the finished math in code form or the notation in simple form so I can encode it.. ​

---------- Post added 07-19-13 at 07:46 AM ---------- Previous post was 07-18-13 at 03:00 PM ----------

' position in cartesian coordinates
' Adjustments to location on the equatorial plane of the planet
Dim xlng AsDouble = r * System.Math.Cos(lng)
Dim ylng AsDouble = r * System.Math.Sin(lng)
Dim zlng AsDouble = 0

' Adjustments to location on the Euler plane (-LAN)
Dim xmLan AsDouble = ((xlng * System.Math.Cos(-LAN)) - (ylng * System.Math.Sin(-LAN)))
Dim ymLan AsDouble = ((xlng * System.Math.Cos(-LAN)) - (ylng * System.Math.Sin(-LAN)))
Dim zmLan AsDouble = zlng

' Adjustments to location on the Euler plane (Inclination)
Dim xinc AsDouble = xmLan
Dim yinc AsDouble = ((ymLan * System.Math.Cos(inc)) - (zmLan * System.Math.Sin(inc)))
Dim zinc AsDouble = ((ymLan * System.Math.Sin(inc)) + (zmLan * System.Math.Cos(inc)))

' Adjustments to location on the Euler plane (+LAN)
Dim xLan AsDouble = ((xinc * System.Math.Cos(LAN)) - (yinc * System.Math.Sin(LAN)))
Dim yLan AsDouble = ((xinc * System.Math.Cos(LAN)) - (yinc * System.Math.Sin(LAN)))
Dim zLan AsDouble = zinc

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

' Adjustments to velocity on the equatorial plane of the planet
Dim xvlng AsDouble = -v * System.Math.Sin(lng)
Dim yvlng AsDouble = -v * System.Math.Cos(lng)
Dim zvlng AsDouble = 0

' Adjustments to velocity on the Euler plane (-LAN)
Dim xmvLan AsDouble = ((xvlng * System.Math.Cos(-LAN)) - (yvlng * System.Math.Sin(-LAN)))
Dim ymvLan AsDouble = ((xvlng * System.Math.Cos(-LAN)) - (yvlng * System.Math.Sin(-LAN)))
Dim zmvLan AsDouble = zvlng

' Adjustments to velocity on the Euler plane (Inclination)
Dim xvinc AsDouble = xmvLan
Dim yvinc AsDouble = ((ymvLan * System.Math.Cos(inc)) - (zmvLan * System.Math.Sin(inc)))
Dim zvinc AsDouble = ((ymvLan * System.Math.Sin(inc)) + (zmvLan * System.Math.Cos(inc)))

' Adjustments to velocity on the Euler plane (+LAN)
Dim xvLan AsDouble = ((xvinc * System.Math.Cos(LAN)) - (yvinc * System.Math.Sin(LAN)))
Dim yvLan AsDouble = ((xvinc * System.Math.Cos(LAN)) - (yvinc * System.Math.Sin(LAN)))
Dim zvLan AsDouble = zvinc

' assuming a standard (non-retrograde) orbit
Dim x AsDouble = xLan
Dim z AsDouble = yLan
Dim y AsDouble = zLan
Dim vx AsDouble = xvLan
Dim vz AsDouble = yvLan
Dim vy AsDouble = zvLan

' RPOS and RVEL values
'RPOS_Out.Text = FormatNumber(r, 2, , , False) & " 0.00 0.00"
RPOS_Out.Text = x & " " & y & " " & z
'RVEL_Out.Text = "0.00 " & FormatNumber(v, 2, , , False) & " 0.00"
RVEL_Out.Text = vx & " " & vy & " " & vz
　

Puts me in a parabolic orbit. Can you show me my math/logic error?

 

[Albinon]

Working code

This code gives an orbit. (finally). Not the one I want, but it is working a lot better.

' Euler X is the Vernal Equinox or Vernal Point
' Euler y is approximately Polaris
' Euler z is between Taurus and Gemini
' this matches the Orbiter system of coordinates.

Dim ToRadian AsDouble = Math.PI / 180

Dim lng AsDouble = Val(Long_In.Text) * ToRadian ' to radians
Dim clong AsDouble = System.Math.Cos(lng)
Dim slng AsDouble = System.Math.Sin(lng)

Dim lat AsDouble = Val(Lat_In.Text) * ToRadian ' to radians
Dim clat AsDouble = System.Math.Cos(lat)
Dim slat AsDouble = System.Math.Sin(lat)

Dim inc AsDouble = Val(OrbitInc.Text) * ToRadian ' to radians
Dim cinc AsDouble = System.Math.Cos(inc)
Dim sinc AsDouble = System.Math.Sin(inc)

Dim LAN AsDouble = Val(OrbitLAN.Text) * ToRadian ' to radians
Dim cmlan AsDouble = System.Math.Cos(-LAN)
Dim smlan AsDouble = System.Math.Sin(-LAN)

Dim clan AsDouble = System.Math.Cos(LAN)
Dim slan AsDouble = System.Math.Sin(LAN)

r = Val(BodyRadius.Text) + Val(OrbitSMA.Text)

' position in cartesian coordinates
Dim x AsDouble = r
Dim y AsDouble = 0
Dim z AsDouble = 0

' Adjustments to location on the equatorial plane of the planet
Dim xlng AsDouble = ((x * clong) + (y * slng)) ' x is 0
Dim ylng AsDouble = ((x * -slng) + (y * clong))
Dim zlng AsDouble = z

' Adjustments to location on the Euler plane (-LAN)
Dim xmLAN AsDouble = ((xlng * cmlan) + (ylng * smlan))
Dim ymLAN AsDouble = ((xlng * -smlan) + (ylng * cmlan))
Dim zmLAN AsDouble = zlng

' Adjustments to location on the Euler plane (Inclination)
Dim xinc AsDouble = xmLAN
Dim yinc AsDouble = ((ymLAN * cinc) + (zmLAN * sinc))
Dim zinc AsDouble = ((ymLAN * -sinc) + (zmLAN * cinc))

' Adjustments to location on the Euler plane (+LAN)
Dim xLAN AsDouble = ((xinc * clan) + (yinc * slan))
Dim yLAN AsDouble = ((xinc * -slan) + (yinc * clan))
Dim zLAN AsDouble = zinc

' Base velocity of orbit
M = Val(BodyMass.Text)
v = System.Math.Sqrt((G * M) / r)

' velocity in cartesian coordinates
Dim vx AsDouble = 0
Dim vy AsDouble = v
Dim vz AsDouble = 0

' Adjustments to location on the equatorial plane of the planet
Dim vxlng AsDouble = ((vx * clong) + (vy * slng)) ' x is 0
Dim vylng AsDouble = ((vx * -slng) + (vy * clong))
Dim vzlng AsDouble = vz

' Adjustments to location on the Euler plane (-LAN)
Dim vxmLAN AsDouble = ((vxlng * cmlan) + (vylng * smlan))
Dim vymLAN AsDouble = ((vxlng * -smlan) + (vylng * cmlan))
Dim vzmLAN AsDouble = vzlng

' Adjustments to location on the Euler plane (Inclination)
Dim vxinc AsDouble = vxmLAN
Dim vyinc AsDouble = ((vymLAN * cinc) + (vzmLAN * sinc))
Dim vzinc AsDouble = ((vymLAN * -sinc) + (vzmLAN * cinc))

' Adjustments to location on the Euler plane (+LAN)
Dim vxLAN AsDouble = ((vxinc * clan) + (vyinc * slan))
Dim vyLAN AsDouble = ((vxinc * -slan) + (vyinc * clan))
Dim vzLAN AsDouble = vzinc
　
' RPOS and RVEL values
RPOS_Out.Text = FormatNumber(xLAN, 2, , , False) & " " & FormatNumber(yLAN, 2, , , False) & " " & FormatNumber(zLAN, 2, , , False)
RVEL_Out.Text = FormatNumber(vxLAN, 2, , , False) & " " & FormatNumber(vyLAN, 2, , , False) & " " & FormatNumber(vzLAN, 2, , , False)


---------- Post added at 08:25 AM ---------- Previous post was at 08:16 AM ----------

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

Forced me to take remedial analytical algebra (something I needed to do anyway) but I finally got matrix math! Yee-Haw!

All I got to do now is figure out why everything is off-kilter. The orbit doesn't match the one in the original scenario.

Thanks to ADSWNJ for getting me this far.:tiphat:

---------- Post added at 10:37 AM ---------- Previous post was at 08:25 AM ----------

All I got to do now is figure out why everything is off-kilter.

Hah!:facepalm:

This code just moves the orbiting body around the "sphere" and applys the proper relative velovity vectors so it stays in a circular orbit. It doesn't actually place it into a specific Keplerian orbit, which is what I want to do.

At least I got the math working so I can fiddle around with euler angles and matrix math

 

[Albinon]

' Velocity in cartesian coordinates
Dim vx AsDouble = 0
Dim vy AsDouble = v
Dim vz AsDouble = 0

Is wrong. Should be:
' pVelocity in cartesian coordinates
Dim vx AsDouble = 0
Dim vy AsDouble = 0
Dim vz AsDouble = v

Still working out how to get a proper Keplerian orbit. Anyone got an idea?

 

[indy91]

Let me break this down again. You want to calculate RPOS and RVEL of an orbit around a planet with the inputs: Longitude, Latitude, Inclination (Ecliptic or equatorial frame?), Longitude of the Ascending Node (same question), Altitude. You calculate the velocity for a circular orbit with that radius. Basically the vehicle is supposed to be directly above a specific point (base or whatever) on the planet in a specific circular orbit.

If you have to convert equatorial coordinates in ecliptic coordinates (I think you have to do this) then you'll need even more calculations. I would suggest reading the precession.pdf in Doc/Technotes. It is quite complicated though. When you create stars and planets, you do have to set the parameters of Table 1 in the document, right? Also do you have access to API functions? Sorry, I am good with the math, the matrices and the orbital mechanics, but not the programming :lol:

Because the problem is RPOS and RVEL are the state vector in the ecliptic frame (I don't know if you can use another format and put the equatorial state vector in the scenario file?) When you have the state vector in the equatorial, rotating frame you have to calculate it first in the non-rotating equatorial frame and then in the ecliptic frame. This is probably the state vector that you have to use in the scenario file. It sounds difficult, but it could be worse

 

[Albinon]

Explan

Also do you have access to API functions?

Not in Visual Studios Visual Basic 2010

Sorry, I am good with the math, the matrices and the orbital mechanics, but not the programming

I am the opposite, good with the programming, not good with the math. I do understand the theory on orbital mechanics, just don't ask me to do the math part.

My problem now is setting up the orientation of the 2 planes. I think I got the X-Y (vernal point - Polaris) because the inclination is coming out correct as long as LAN and AgP plus TrA are both 0. However, the X-Z (vernal point - Gemini/Taurus) is proving to be difficult. The orbit is circular, but weird numbers show up in the Orbit MFD. I've tried both X-Z and Y-Z orientations but both will not set up correctly.

 

[indy91]

I have still these questions: Were all assumptions in my first paragraph right? Do you have to set SidRotPeriod and Obliquity etc. when you create the solar systems? Because these parameters are needed to calculate the Rotation matrix for the coordinate transformation. If the answers to the questions are yes, then the next paragraphs will make sense (kind of):

Let me start at the beginning. I think you first have to get the spherical polar coordinates. That's easy, because it is just the three components lat, lng and radius.

Then you have to transform these in cartesian coordinates. This is similar to what you already do to calculate e.g. ylng (and xlng and zlng).

Now you have to do the same calculations as the function oapiGetPlanetObliquityMatrix does and transform the coordinates from equatorial in ecliptic coordinates. Here you need all the information from the config file of the planet as specified under "Rotation and precession parameters". I don't know how your program will work, but I guess it creates these information? I can explain you in detail what you would have to do.

The velocity is another story. I think you don't even need the LAN. With a given position and the velocity there are 0,1 or 2 possibilites for the LAN. For example, an inclination of 30° with the latitude of 50° would never work, because the spacecraft will never go so high north/south.

I don't know if I am overcomplicating things, but I had to do all of this for a lua script, where the nice API functions aren't avaiable either. I hope I could/can help.

 

[Albinon]

Yes, the planets are set up using normal Orbiter parameters. I place them into a cfg file according to orbiter specifications. The parameters are stored in a database that my program has access to. The system, star, planet and moons work normally.

As a programmer, I can work out the processes visually if I see it. Is it possible to see your LUA script? I don't know that programming language, but I pick up new programming languages fairly quickly. (I have been a programmer since 1979 and know a "few" languages, even some that are now considered obsolete)

 

[indy91]

Ok, let's say I felt challenged to solve your problem. Here is a lua script (I had to change it to txt so I could upload it) that should do exactly what you want. I hope you can use it/rewrite it for your code.

 

[Albinon]

Ok, let's say I felt challenged to solve your problem. Here is a lua script (I had to change it to txt so I could upload it) that should do exactly what you want. I hope you can use it/rewrite it for your code.

OK, I'll give it a shot and let you know.

Also, I did a bit of research on my own.
Since most orbits are considered ellipses, I got velocity as

G = 0.0000000000667 ' Gravitational Constant
a = Val(OrbitSMA.Text) 'semi-major axis (furthest from center of orbit)
M = Val(BodyMass.Text)
U = G * M
v = System.Math.Sqrt(U * ((2 / math.sqrt(x^2 + y^2 + z^2)) - (1 / a))) ' for a elliptical orbit


With more research I got the following 2D formula for a point on an ellipse
'x = acos(θ)
'y = rsin(θ)

I am now working on setting up a 3D matrix to use. Once that is done and I have translated your script (Thankyou!), I'll see what comes out the other side and post my findings back here.

 

[Albinon]

Ran into a circular argument...

a=(b√(1-e^2 ))/(1-e^2 )
e=√(1+(2Eh^2)/u^2 )
E=v^2/2-u/√(x^2+y^2+z^2 )
h=√(a*u*(1-e^2))

It seems that a semi-major axis can't be calculated with out already knowing what it is...

Does anyone have a different method?

 

[indy91]

Well, there are different ways to calculate the SMA. Here are some ways, maybe one is helpful for you:

usual definition: a = (r_p+r_a)/2, r_a and r_p are apoapsis and periapsis radius
orbital elements: a = (h*h/mu)*1/(1-e*e) with mu being the gravitational parameter
through orbital period: a=(T*sqrt(mu)/(2*PI))^(2/3)

 

[Albinon]

Well, there are different ways to calculate the SMA. Here are some ways, maybe one is helpful for you:

usual definition: a = (r_p+r_a)/2, r_a and r_p are apoapsis and periapsis radius
orbital elements: a = (h*h/mu)*1/(1-e*e) with mu being the gravitational parameter
through orbital period: a=(T*sqrt(mu)/(2*PI))^(2/3)

I thank you for this, however research shows me that all three methods require the knowledge of the semi-major axis (I could be wrong, I'm pretty weak in math). I am trying to calculate it without measurements or orbiter's API suite.

r_a (Apoapsis) = a(1+e)
h =√(a*u*(1-e^2))
T (orbital period) = 2*pi*√(a^3/u)

 

[indy91]

Ok, then I don't quite understand what parameters you have to calculate the SMA. Usually (and in Orbiter) it is one of the basic orbital parameters and other parameters are calculated with the SMA. What parameters do you have to calculate it from?

 

[Albinon]

Ok, then I don't quite understand what parameters you have to calculate the SMA. Usually (and in Orbiter) it is one of the basic orbital parameters and other parameters are calculated with the SMA. What parameters do you have to calculate it from?

Not many, just the basic calculation that puts it at RPOS r, 0, 0 and RVEL 0, 0, v with r giving me semi-minor, LAN and AgP/TrA

I was trying to over-simplify and calculate it from those, but I now understand that I cannot cut corners.

Remember, math is not my forte'. I need a "pretty picture" to understand what is going on.

I'm about half done translating your script. It has given me some interesting programming ideas (VB 2010 is not well adapted to matrix math, your LUA script showed me alternatives)

---------- Post added at 02:27 PM ---------- Previous post was at 01:16 PM ----------

Ok, let's say I felt challenged to solve your problem. Here is a lua script (I had to change it to txt so I could upload it) that should do exactly what you want. I hope you can use it/rewrite it for your code.

OK, VB 2010 does not have an equivalent to the following statements and I definitely don't have the math to work them out. What are they? and how do they work? I believe they are some form of matrix math.

local Q_Xp = mat.mmul(R2,R3_W)
local V_equ = mat.tmul(Q_Xp,vp)
local RPOS = mat.mul(Rot,R_equ)
local VPOS = mat.mul(Rot,V_equ)

As for your Descending "beta" I think the answer is you have to transpose the math as the angle is set from a different direction. You'd have to test for a positive or negative Inc to do it. (Tell me if I'm wrong)
I haven't had the chance to work it out yet, but I think its simple trig and within what I can do.

 

[indy91]

Not many, just the basic calculation that puts it at RPOS r, 0, 0 and RVEL 0, 0, v with r giving me semi-minor, LAN and AgP/TrA

I'm not sure, but I think you are confusing semi-minor axis with periapsis and semi-major axis with apoapsis. These are not the same. Remember, we are talking about ellipses, however the planet is in one of the foci(or focus' ?) of the ellipse and not in the center.

Mathwise you could use my second formula a = (h*h/mu)*1/(1-e*e).

H=RxV and h=norm(H), so h is the length of the crossproduct of the R and V vectors.

E=(VxH)/mu-R/r. This is two vector subtracted and consists of the crossproduct of V and H and R/r would be the unit vector (length 1) of the radius vector. r is again just the length of R. e=norm(E) so again the length of the vector. Now you can calculate a.

In this notation a vector is always in upper case and the length of a vector (or another parameter) is usually in lower case.

OK, VB 2010 does not have an equivalent to the following statements and I definitely don't have the math to work them out. What are they? and how do they work? I believe they are some form of matrix math.

mat.mmul is two matrices multiplicated, mat.mul is a matrix multiplicated with a vector and mat.tmul is the transposed matrix multiplicated with a vector.