Question Equatorial inc in lua

pclaurent

Daydreamer
Joined
Dec 21, 2014
Messages
49
Reaction score
0
Points
0
I try to get the equatorial inclination of my orbit within an LUA script. The LUA implementation of function GetElements does not seem to support the frame parameter, and the only elements I get are the ecliptical ones. So no option for getting this value directly from OAPI library... Any idea?
 

pclaurent

Daydreamer
Joined
Dec 21, 2014
Messages
49
Reaction score
0
Points
0
Thanks very much for the link. The algorithm seems to perfectly fit my needs.

However the example code does'nt work as is. Matrix.lua requires another module, complex.lua. And, even if I copy these 2 modules (using the full matrix package, see attached), the execution fails at first line :
Rref = matrix{{1,0,0},{0,1,0},{0,0,1}}

As my knowledge in lua is somehow limited, perhaps I missed somehting obvious to more experienced people.

If someone has a working example of code (if possible in a stand-alone form), it would be more than cool...
 

Attachments

  • matrix-0.2.8.zip
    20.4 KB · Views: 8
Last edited:

dgatsoulis

ele2png user
Donator
Joined
Dec 2, 2009
Messages
1,924
Reaction score
340
Points
98
Location
Sparta
Ah yes, I had forgotten about the complex.lua dependency.

It works fine for me.
I placed matrix.lua and complex.lua in the Orbiter root directory, and used this script in a scenario:
(same as before, I just added a "note" in the end to show the eq Inc and LAN on screen.

Code:
matrix = require "matrix"
v = vessel.get_focusinterface()
function Get_Target_Elements_eq()
Rref = matrix{{1,0,0},{0,1,0},{0,0,1}}
Erel = 0.4090928023 --obliquity [rad]
T_p = -9413040.4 --precession period [days]
T_0 = 51544.5 --LAN_MJD
T_s = 0.99726968 --Siderial rotation period [days]
f_0 = 4.88948754 --SidRotOffset [rad]
L_0 = 0 --LAN
T = oapi.get_simmjd()
Lrel = L_0 + (2*PI*(T-T_0)/T_p)
Rrel_1 = matrix{{math.cos(Lrel), 0 , -math.sin(Lrel)},{0,1,0},{math.sin(Lrel),0,math.cos(Lrel)}}
Rrel_2 = matrix{{1,0,0},{0,math.cos(Erel),-math.sin(Erel)},{0,math.sin(Erel),math.cos(Erel)}}
Rrel = matrix.mul(Rrel_1,Rrel_2)
RrefRrel = matrix.mul(Rref,Rrel) 
f = f_0 + (2*PI*(T-T_0)/T_s)+(L_0 - Lrel)*math.cos(Erel)
Rrot = matrix{{math.cos(f),0,-math.sin(f)},{0,1,0},{math.sin(f),0,math.cos(f)}}
R = matrix.mul(RrefRrel,Rrot)
V1 = {{0},{1},{0}}
OS_st = matrix.mul(R,V1)
E_ecl = math.acos(OS_st[2][1])
L_ecl = math.atan(-OS_st[1][1]/OS_st[3][1])
E = E_ecl + L_ecl
el,op = v:get_elementsex()
Inc = el['i']
LAN = el['theta']
omegab = el['omegab']
MnA = op['MnA']
MnL = el['L']
AgP = MnL - LAN - MnA
if LAN < 0 then LAN = LAN + 2*PI end
if AgP < 0 then AgP = AgP + 2*PI end
if MnA < 0 then MnA = MnA + 2*PI end
if MnL < 0 then MnL = MnL + 2*PI end
I_e = math.acos(-math.sin(E_ecl)*math.sin(Inc)*math.cos(LAN-L_ecl)+math.cos(E_ecl)*math.cos(Inc))
A = math.sqrt(math.sin(E_ecl)*math.sin(E_ecl)*(math.cos(Inc)*math.cos(Inc)+math.sin(Inc)*math.sin(Inc)*math.sin(LAN-L_ecl)*math.sin(LAN-L_ecl))+0.5*math.sin(2*E_ecl)*math.sin(2*Inc)*math.cos(LAN-L_ecl)+math.cos(E_ecl)*math.cos(E_ecl)*math.sin(Inc)*math.sin(Inc))
LAN_e = math.acos((math.sin(E_ecl)*math.cos(Inc)+math.cos(E_ecl)*math.sin(Inc)*math.cos(LAN-L_ecl))/A)
BETA = math.acos((math.sin(E_ecl)*math.cos(Inc)*math.cos(LAN-L_ecl)+math.cos(E_ecl)*math.sin(Inc))/A)
B = math.sin(E_ecl)*math.sin(Inc)*math.sin(LAN-L_ecl)
if B < 0 then B = -1 end
if B == 0 then B = 0 end
if B > 0 then B = 1 end
BETA = BETA*B
AgP_e = AgP+BETA 
MnL_e = MnA + LAN_e + AgP_e
TrL_e = op['TrL'] + MnL_e - MnL
check = math.sin(Inc)*math.sin(LAN-L_ecl)/A
if check < 0 then LAN_e = 2*PI - LAN_e end
if AgP_e < 0 then AgP_e = AgP_e + 2*PI end
if AgP_e > 2*PI then AgP_e = AgP_e - 2*PI end
if MnL_e < 0 then MnL_e = MnL_e + 2*PI end
if LAN_e > PI then MnL_e = 2*PI - MnL_e + PI end
if LAN_e > PI then TrL_e = 2*PI - TrL_e + PI end
if MnL_e > 2*PI then MnL_e = MnL_e - 2*PI end
if TrL_e > 2*PI then TrL_e = TrL_e - 2*PI end
if TrL_e < 0 then TrL_e = TrL_e + 2*PI end
SMa = el['a']
Ecc = el['e']
return I_e, LAN_e, AgP_e, MnL_e, TrL_e, SMa, Ecc
end

note = oapi.create_annotation()
note:set_pos (0.7,0.07,1,1)
note:set_colour ({r=76/255,g=207/255,b=1})
note:set_size (1)

goal = 0
while goal < 1 do
I_e, LAN_e, AgP_e, MnL_e, TrL_e, SMa, Ecc = Get_Target_Elements_eq()
msg = string.format('Inc:%.2f°\nLAN:%.2f°',I_e*DEG,LAN_e*DEG)
note:set_text(msg)
proc.skip()
if goal > 0 then end
end

Link to a screenshot.

https://ibb.co/gVX485
 

pclaurent

Daydreamer
Joined
Dec 21, 2014
Messages
49
Reaction score
0
Points
0
OK, works well now. Just added LPe_e = LAN_e+AgP_e to complete the returned elements.
Many thanks..........
 

pclaurent

Daydreamer
Joined
Dec 21, 2014
Messages
49
Reaction score
0
Points
0
FYI, the previously posted code had a bug in managing the position of the LAN. After a few tries, I made a fix. Now the function returns proper elements in all cases - as far as I can see -.

matrix = require "matrix"
v = vessel.get_focusinterface()

function Get_Target_Elements_eq()
Rref = matrix{{1,0,0},{0,1,0},{0,0,1}}
Erel = 0.4090928023 --obliquity [rad]
T_p = -9413040.4 --precession period [days]
T_0 = 51544.5 --LAN_MJD
T_s = 0.99726968 --Sidereal rotation period [days]
f_0 = 4.88948754 --SidRotOffset [rad]
L_0 = 0 --LAN
T = oapi.get_simmjd()
Lrel = L_0 + (2*PI*(T-T_0)/T_p)
Rrel_1 = matrix{{math.cos(Lrel), 0 , -math.sin(Lrel)},{0,1,0},{math.sin(Lrel),0,math.cos(Lrel)}}
Rrel_2 = matrix{{1,0,0},{0,math.cos(Erel),-math.sin(Erel)},{0,math.sin(Erel),math.cos(Erel)}}
Rrel = matrix.mul(Rrel_1,Rrel_2)
RrefRrel = matrix.mul(Rref,Rrel)
f = f_0 + (2*PI*(T-T_0)/T_s)+(L_0 - Lrel)*math.cos(Erel)
Rrot = matrix{{math.cos(f),0,-math.sin(f)},{0,1,0},{math.sin(f),0,math.cos(f)}}
R = matrix.mul(RrefRrel,Rrot)
V1 = {{0},{1},{0}}
OS_st = matrix.mul(R,V1)
E_ecl = math.acos(OS_st[2][1])
L_ecl = math.atan(-OS_st[1][1]/OS_st[3][1])
E = E_ecl + L_ecl
el,op = v:get_elementsex()
Inc = el['i']
LAN = el['theta']
omegab = el['omegab']
MnA = op['MnA']
MnL = el['L']
AgP = MnL - LAN - MnA
if LAN < 0 then LAN = LAN + 2*PI end
if AgP < 0 then AgP = AgP + 2*PI end
if MnA < 0 then MnA = MnA + 2*PI end
if MnL < 0 then MnL = MnL + 2*PI end
I_e = math.acos(-math.sin(E_ecl)*math.sin(Inc)*math.cos(LAN-L_ecl)+math.cos(E_ecl)*math.cos(Inc))
A = math.sqrt(math.sin(E_ecl)*math.sin(E_ecl)*(math.cos(Inc)*math.cos(Inc)+math.sin(Inc)*math.sin(Inc)*math.sin(LAN-L_ecl)*math.sin(LAN-L_ecl))+0.5*math.sin(2*E_ecl)*math.sin(2*Inc)*math.cos(LAN-L_ecl)+math.cos(E_ecl)*math.cos(E_ecl)*math.sin(Inc)*math.sin(Inc))
LAN_e = math.acos((math.sin(E_ecl)*math.cos(Inc)+math.cos(E_ecl)*math.sin(Inc)*math.cos(LAN-L_ecl))/A)
BETA = math.acos((math.sin(E_ecl)*math.cos(Inc)*math.cos(LAN-L_ecl)+math.cos(E_ecl)*math.sin(Inc))/A)
B = math.sin(E_ecl)*math.sin(Inc)*math.sin(LAN-L_ecl)
if B < 0 then B = -1 end
if B == 0 then B = 0 end
if B > 0 then B = 1 end
BETA = BETA*B
AgP_e = AgP+BETA
check = math.sin(Inc)*math.sin(LAN-L_ecl)/A
if check < 0 then LAN_e = 2*PI - LAN_e end
LPe_e = LAN_e+AgP_e
MnL_e = MnA + LAN_e + AgP_e
TrL_e = op['TrL'] + MnL_e - MnL
AgP_e = AgP_e - 2*PI*math.floor(AgP_e/2/PI)
LPe_e = LPe_e - 2*PI*math.floor(LPe_e/2/PI)
MnL_e = MnL_e - 2*PI*math.floor(MnL_e/2/PI)
TrL_e = TrL_e - 2*PI*math.floor(TrL_e/2/PI)

return I_e, LAN_e, AgP_e, LPe_e, MnL_e, TrL_e
end
 
Last edited:
Top