Moon.cfg in orbiter beta

[jarmonik]

Looks like I may have fundamentally misunderstood the point here. So, do you mean that orbiter is using L_{ref}, \eps_{ref}, L_{rel} and \eps_{rel} just to compute L_{ecl} and \eps_{ecl}.

And then uses L_{ecl} and \eps_{ecl} to define a planet orientation and computes the planet rotation relative to L_{ecl}

If that's the case, then cos(\eps_{rel}) may need to be applied during computation of L_{ecl} but not necessarily in the final value of L_{ecl}. I'll need to think about this.

 

[martins]

This will lead into a suggestion...

Any thoughts ?

The thing that bugs me about this equation is that it doesn't seem to be cyclic in longitude. To see what I mean, for simplicity assume T_s = infinite and Delta \phi_{t_0} = 0. Then

\phi(t) = - 2\pi (t-t_0)/T_p cos \eps_rel

and thus

\phi(t_0) = 0

Now I would expect \phi to revert to its original value (0, or a multiple of 2\pi) after one precession cycle, i.e.

\phi(t_0 + T_p) = 0

but in your equation that doesn't seem to be the case. Instead, you get an offset of 2\pi cos \eps_rel, and I just can't figure out where this would come from.

 

[martins]

Looks like I may have fundamentally misunderstood the point here. So, do you mean that orbiter is using L_{ref}, \eps_{ref}, L_{rel} and \eps_{rel} just to compute L_{ecl} and \eps_{ecl}.

And then uses L_{ecl} and \eps_{ecl} to define a planet orientation and computes the planet rotation relative to L_{ecl}

Yes, precisely.

 

[jarmonik]

Yes, precisely.

In that case the equation 7 would become:

I also tried the lunar model from GetTRDoc.pdf but there is a reference frame problem. Any ideas what are the main axies of ICRF reference frame in the Orbiter ? I suppose it shoud be pretty close to equator of the Earth.

This is the AGC model I have used. (It's right-handed)

[SIZE=2][COLOR=#0000ff][SIZE=2][COLOR=#0000ff]#define[/COLOR][/SIZE][/COLOR][/SIZE][SIZE=2] _I_ECL _V(1, 0, 0) [/SIZE]
[SIZE=2][COLOR=#0000ff][SIZE=2][COLOR=#0000ff]#define[/COLOR][/SIZE][/COLOR][/SIZE][SIZE=2] _K_ECL _V(0, 0, 1) [/SIZE]
[SIZE=2][COLOR=#0000ff][SIZE=2][COLOR=#0000ff]#define[/COLOR][/SIZE][/COLOR][/SIZE][SIZE=2] _J_ECL _V(0, 1, 0) [/SIZE]
 
[SIZE=2][COLOR=#0000ff][SIZE=2][COLOR=#0000ff]#define[/COLOR][/SIZE][/COLOR][/SIZE][SIZE=2] lnEpoch 40038.000000[/SIZE]
[SIZE=2][COLOR=#0000ff][SIZE=2][COLOR=#0000ff]#define[/COLOR][/SIZE][/COLOR][/SIZE][SIZE=2] lnRad 1738090.0[/SIZE]
[SIZE=2][COLOR=#0000ff][SIZE=2][COLOR=#0000ff]#define[/COLOR][/SIZE][/COLOR][/SIZE][SIZE=2] lnRotOff 5.755387533590 [/SIZE]
[SIZE=2][COLOR=#0000ff][SIZE=2][COLOR=#0000ff]#define[/COLOR][/SIZE][/COLOR][/SIZE][SIZE=2] lnRotPer 2351139.374       // node to node rotation period[/SIZE]
[SIZE=2][COLOR=#0000ff][SIZE=2][COLOR=#0000ff]#define[/COLOR][/SIZE][/COLOR][/SIZE][SIZE=2] lnLanOff 0.250691117180 [/SIZE]
[SIZE=2][COLOR=#0000ff][SIZE=2][COLOR=#0000ff]#define[/COLOR][/SIZE][/COLOR][/SIZE][SIZE=2] lnLanPer -586955667.0[/SIZE]
[SIZE=2][COLOR=#0000ff][SIZE=2][COLOR=#0000ff]#define[/COLOR][/SIZE][/COLOR][/SIZE][SIZE=2] lnObliquity 0.0267907[/SIZE]
 
 
[SIZE=2][COLOR=#0000ff][SIZE=2][COLOR=#0000ff]double[/COLOR][/SIZE][/COLOR][/SIZE][SIZE=2] t = (oapiGetSimMJD() - lnEpoch) * 86400.0;[/SIZE]
[SIZE=2][COLOR=#0000ff][SIZE=2][COLOR=#0000ff]double[/COLOR][/SIZE][/COLOR][/SIZE][SIZE=2] rt = limit(lnRotOff + PI2 * t / lnRotPer);[/SIZE]
[SIZE=2][COLOR=#0000ff][SIZE=2][COLOR=#0000ff]double[/COLOR][/SIZE][/COLOR][/SIZE][SIZE=2] ob = lnObliquity;[/SIZE]
[SIZE=2][COLOR=#0000ff][SIZE=2][COLOR=#0000ff]double[/COLOR][/SIZE][/COLOR][/SIZE][SIZE=2] ln = limit(lnLanOff + PI2 * t / lnLanPer); [/SIZE]
 
[SIZE=2]VECTOR3 _R; [/SIZE]
[SIZE=2]hVessel->GetRelativePos(hMoon, _R); [/SIZE]
 
[SIZE=2]_R = _V(_R.x, _R.z, _R.y);[/SIZE]
 
[SIZE=2]VECTOR3 _i = _I_ECL*cos(ln) + _J_ECL*sin(ln);[/SIZE]
[SIZE=2]VECTOR3 _j = unit(crossp(_K_ECL, _i));[/SIZE]
 
[SIZE=2]VECTOR3 _K = _K_ECL*cos(ob) + _j*sin(ob);[/SIZE]
[SIZE=2]VECTOR3 _I = _i*cos(rt) + unit(crossp(_K, _i))*sin(rt);[/SIZE]
[SIZE=2]VECTOR3 _J = unit(crossp(_K, _I));[/SIZE]
 
[SIZE=2]_R = _V( dotp(_R,_I), dotp(_R,_J), dotp(_R,_K));[/SIZE]
 
[SIZE=2][COLOR=#0000ff][SIZE=2][COLOR=#0000ff]double[/COLOR][/SIZE][/COLOR][/SIZE][SIZE=2] lat = atan( _R.z / sqrt(_R.x*_R.x + _R.y*_R.y) );[/SIZE]
[SIZE=2][COLOR=#0000ff][SIZE=2][COLOR=#0000ff]double[/COLOR][/SIZE][/COLOR][/SIZE][SIZE=2] lng = atan2( _R.y, _R.x );[/SIZE]
 
[SIZE=2]TextA(5,pos,[/SIZE][SIZE=2][COLOR=#a31515][SIZE=2][COLOR=#a31515]"Lat "[/COLOR][/SIZE][/COLOR][/SIZE][SIZE=2],DEG*lat,6); pos+=ld;[/SIZE]
[SIZE=2]TextA(5,pos,[/SIZE][SIZE=2][COLOR=#a31515][SIZE=2][COLOR=#a31515]"Lng "[/COLOR][/SIZE][/COLOR][/SIZE][SIZE=2],DEG*lng,6); pos+=ld;[/SIZE]

Moon configuration file parameters:

SidRotPeriod = 2360591.672      
;SidRotOffset = 4.765658435       
SidRotOffset = 0.657331652
 
LAN = 2.174858298    
Obliquity = 0.0267907
PrecessionPeriod = -6793.4683    ; precession period (days)
LAN_MJD = 51544.5           ; reference date for LAN

 

[martins]

In that case the equation 7 would become:

Ok, but if we look at the case L_{ref} = 0, \eps_{ref} = 0, then this equation reverts to the previous one, and my question still stands: If we assume T_s = \infty and \Delta \phi_{t_0}=0, then why is \phi(t_0 + T_p) != 0 ?

 

[jarmonik]

Ok, but if we look at the case L_{ref} = 0, \eps_{ref} = 0, then this equation reverts to the previous one, and my question still stands: If we assume T_s = \infty and \Delta \phi_{t_0}=0, then why is \phi(t_0 + T_p) != 0 ?

You are absolutely right, if a planet doesn't rotate it should be right back where it started from, after passing N-number of full precession periods. That's a very good point so there must be a flaw in that equation.

---------- Post added at 04:11 AM ---------- Previous post was at 03:10 AM ----------

It looks like cos(\eps_{rel}) is not a factor in the equation. It just compensated an other error and provided pretty good results.

The node to node rotation period used in AGC is T_nn = 2351139.374, Sidereal rotation period can be computed from:

SidRotPeriod = (1/T_nn + 1/T_p)^{-1} = 2360595.093

So, this is the value that need to be used in Moon.cfg if syncronization with AGC is desired. This is 3.496 seconds higher than initial IAU J2000 value. Effect of cos(\eps_{rel}) term was 3.551 seconds and that mostly compensated the previous error.

 

[martins]

It looks like cos(\eps_{rel}) is not a factor in the equation. It just compensated an other error and provided pretty good results.

Ah, good. I was scratching my head if I had missed something.

So now I guess the question is which value for SidRotPeriod to use. For what it's worth, as a sanity check I tried both values in orbiter and projected them 10 000 years into the future. The IAU value seems pretty stable, while the AGC value introduces a significant rotation offset over this period.

 

[jarmonik]

So now I guess the question is which value for SidRotPeriod to use. For what it's worth, as a sanity check I tried both values in orbiter and projected them 10 000 years into the future. The IAU value seems pretty stable, while the AGC value introduces a significant rotation offset over this period.

That's true. Rotation period used by AGC was about 4 seconds longer than the true one. Obviously they didn't know it more accurately at that time.

In theory AGC should compensate the longer rotation period with different rotation offset in order to make it accurate near it's own epoch of 01-Jul-1968. So, that's the point where the true model and AGC model should actually meet each-other. But that's currently not the case with the values I have.

NASSP will distribute it's own version of Moon.cfg anyway mostly due to surface bases and observers. Therefore, parameters those are matching the AGC requirements isn't top priority.

...I tried both values in orbiter and projected them 10 000 years into the future. The IAU value seems pretty stable...

What values did you use ?

 

[martins]

What values did you use ?

SidRotPeriod = 2360591.597 from my default Moon.cfg (actually I am not even sure if this really is an IAU value or where else I picked that up), and

SidRotPeriod = 2360595.093 as per your previous post.

For both, SidRotOffset = 0.657362170. (This is different from the pre-precession support value of 4.78084319 because of the subtraction of the L_{ecl} value), but in any case I am not at all confident about this one, so unless there are documented values for it, I would be happy to adjust it so that it agrees with the AGC values in 1968.

I haven't studied the IAU document that BrianJ mentioned yet, but would that give an answer to the correct rotation offset?

 

[jarmonik]

I haven't studied the IAU document that BrianJ mentioned yet, but would that give an answer to the correct rotation offset?

Yes, I got some preliminary results from that model but I'll still need to implement the PA-ME rotations and find out what is the exact orientation of the Earth's mean polar axis respect to ecliptic.

These are the current values from IAU document but they could still change a pit:

SidRotPeriod 2360591.5721 (2PI/dW) = (2PI / 2.66169945763e-6)
SidRotOffset 0.6688486548 (4.7808379472 + 2.1711960153)
Obliquity 0.0274042311
LAN_MJD 51544.5
LAN 2.1711960153


These are the values those should match with the AGC in 01-Jul-1968, but I'll still need to double check them.

PrecessionPeriod -6793.4683
SidRotPeriod 2360591.597
SidRotOffset 0.66124264
Obliquity 0.0267907
LAN 2.174858394
LAN_MJD 51544.5

Here are the values for lunar model but the rest of the implementation is still somewhat unclear.

[SIZE=2][COLOR=#0000ff][SIZE=2][COLOR=#0000ff]double[/COLOR][/SIZE][/COLOR][/SIZE][SIZE=2] d = oapiGetSimMJD() - J2000;[/SIZE]
[SIZE=2][COLOR=#0000ff][SIZE=2][COLOR=#0000ff]double[/COLOR][/SIZE][/COLOR][/SIZE][SIZE=2] T = d / 36525.0;[/SIZE]
 
[SIZE=2][COLOR=#0000ff][SIZE=2][COLOR=#0000ff]double[/COLOR][/SIZE][/COLOR][/SIZE][SIZE=2] E1 = (125.045 - 0.0529921*d)*RAD; [/SIZE]
[SIZE=2][COLOR=#0000ff][SIZE=2][COLOR=#0000ff]double[/COLOR][/SIZE][/COLOR][/SIZE][SIZE=2] E2 = (250.089 - 0.1059842*d)*RAD;[/SIZE]
[SIZE=2][COLOR=#0000ff][SIZE=2][COLOR=#0000ff]double[/COLOR][/SIZE][/COLOR][/SIZE][SIZE=2] E3 = (260.008 + 13.0120009*d)*RAD; [/SIZE]
[SIZE=2][COLOR=#0000ff][SIZE=2][COLOR=#0000ff]double[/COLOR][/SIZE][/COLOR][/SIZE][SIZE=2] E4 = (176.625 + 13.3407154*d)*RAD;[/SIZE]
[SIZE=2][COLOR=#0000ff][SIZE=2][COLOR=#0000ff]double[/COLOR][/SIZE][/COLOR][/SIZE][SIZE=2] E5 = (357.529 + 0.9856003*d)*RAD;[/SIZE]
[SIZE=2][COLOR=#0000ff][SIZE=2][COLOR=#0000ff]double[/COLOR][/SIZE][/COLOR][/SIZE][SIZE=2] E6 = (311.589 + 26.4057084*d)*RAD;[/SIZE]
[SIZE=2][COLOR=#0000ff][SIZE=2][COLOR=#0000ff]double[/COLOR][/SIZE][/COLOR][/SIZE][SIZE=2] E7 = (134.963 + 13.0649930*d)*RAD; [/SIZE]
[SIZE=2][COLOR=#0000ff][SIZE=2][COLOR=#0000ff]double[/COLOR][/SIZE][/COLOR][/SIZE][SIZE=2] E8 = (276.617 + 0.3287146*d)*RAD;[/SIZE]
[SIZE=2][COLOR=#0000ff][SIZE=2][COLOR=#0000ff]double[/COLOR][/SIZE][/COLOR][/SIZE][SIZE=2] E9 = (34.266 + 1.7484877*d)*RAD; [/SIZE]
[SIZE=2][COLOR=#0000ff][SIZE=2][COLOR=#0000ff]double[/COLOR][/SIZE][/COLOR][/SIZE][SIZE=2] E10 = (15.134 - 0.1589763*d)*RAD;[/SIZE]
[SIZE=2][COLOR=#0000ff][SIZE=2][COLOR=#0000ff]double[/COLOR][/SIZE][/COLOR][/SIZE][SIZE=2] E11 = (119.743 + 0.0036096*d)*RAD; [/SIZE]
[SIZE=2][COLOR=#0000ff][SIZE=2][COLOR=#0000ff]double[/COLOR][/SIZE][/COLOR][/SIZE][SIZE=2] E12 = (239.961 + 0.1643573*d)*RAD;[/SIZE]
[SIZE=2][COLOR=#0000ff][SIZE=2][COLOR=#0000ff]double[/COLOR][/SIZE][/COLOR][/SIZE][SIZE=2] E13 = (25.053 + 12.9590088*d)*RAD;[/SIZE]
 
[SIZE=2][COLOR=#0000ff][SIZE=2][COLOR=#0000ff]double[/COLOR][/SIZE][/COLOR][/SIZE][SIZE=2] a0 = (269.9949 + 0.0031*T - 3.8787*sin(E1) - 0.1204*sin(E2)[/SIZE]
[SIZE=2]+ 0.0700*sin(E3) - 0.0172*sin(E4) + 0.0072*sin(E6)[/SIZE]
[SIZE=2]- 0.0052*sin(E10) + 0.0043*sin(E13))*RAD;[/SIZE]
 
[SIZE=2][COLOR=#0000ff][SIZE=2][COLOR=#0000ff]double[/COLOR][/SIZE][/COLOR][/SIZE][SIZE=2] r0 = (66.5392 + 0.0130*T + 1.5419*cos(E1) + 0.0239*cos(E2) [/SIZE]
[SIZE=2]- 0.0278*cos(E3) + 0.0068*cos(E4) - 0.0029*cos(E6)[/SIZE]
[SIZE=2]+ 0.0009*cos(E7) + 0.0008*cos(E10) - 0.0009*cos(E13))*RAD;[/SIZE]
 
[SIZE=2][COLOR=#0000ff][SIZE=2][COLOR=#0000ff]double[/COLOR][/SIZE][/COLOR][/SIZE][SIZE=2] W0 = (38.3213 + 13.17635815*d - 1.4e-12*d*d + 3.5610*sin(E1)[/SIZE]
[SIZE=2]+ 0.1208*sin(E2) - 0.0642*sin(E3) + 0.0158*sin(E4) [/SIZE]
[SIZE=2]+ 0.0252*sin(E5) - 0.0066*sin(E6) - 0.0047*sin(E7)[/SIZE]
[SIZE=2]- 0.0046*sin(E8) + 0.0028*sin(E9) + 0.0052*sin(E10)[/SIZE]
[SIZE=2]+ 0.0040*sin(E11) + 0.0019*sin(E12) - 0.0044*sin(E13))*RAD;[/SIZE]

 

[Chode]

Hi Martin,

I've started some testing of the new precession calculations, and I've found what looks to be a bug. If you load up this scenario, you can see what I found. It looks like the body orientation is not being calculated correctly for bodies that do not supply precession data. I've tried this scenario in some previous betas and release, and they all show the correct orientation, which is that Phobos should have it's plus x axis pointed at Mars.

Regards

 

[Chode]

OK, I think I see what needs to be fixed. With the new precession calculation, SidRotOffset needs to have LAN added to it to compensate for the -L_rel term in the calculation of the rotation angle. If I edit the config, adding LAN to SidRotOffset, everything seems to be OK.

Regards

 

[martins]

Yes, precisely. I think it will be better if I add the reference LAN directly in the Orbiter code, so that the config files don't need to be modified (as Jarmonik suggested above).

So instead of just subtracting the current LAN from the rotation offset, the next beta will also add the LAN at the reference date:

rot_offset = ... +LAN_{ecl}(t_0) - LAN_{ecl}(t)

where t_0 is given by the LAN_MJD config entry, and t is the current time. For bodies which don't define a precession, LAN is time-invariant, and the two terms cancel.


Edit: Also, I am beginning to realise that L_{ecl} (longitude with respect to ecliptic frame) in equation 7 is a mistake. The rotation should be independent of the inertial frame of the observer, so one can always select a frame in which \eps_{ref} = 0. Therefore, L_{ecl} should be replaced by L_{rel} in Eq. 7.

 

[jarmonik]

Based on the lunar model in IAU documentation there seems to be periodic changes in Obliquity and LAN. Obliquity will change between 0.0273 and 0.0265 in periods of 16 days (approximately) and the LAN will drift about +-1.0 degrees from the mean value.

There seems be an error of 0.42 degrees between the AGC model and the IAU model. Due to precession, the node of the Earth's equator and ecliptic has changed about the same amount since 1968.

Originally Posted by martins
Edit: Also, I am beginning to realise that L_{ecl} (longitude with respect to ecliptic frame) in equation 7 is a mistake. The rotation should be independent of the inertial frame of the observer, so one can always select a frame in which \eps_{ref} = 0. Therefore, L_{ecl} should be replaced by L_{rel} in Eq. 7.

Sounds like a good idea.

 

[martins]

There seems be an error of 0.42 degrees between the AGC model and the IAU model. Due to precession, the node of the Earth's equator and ecliptic has changed about the same amount since 1968.

Does that mean that the discrepancy could arise from the fact that the AGC values may not be referenced to the J2000 ecliptic, but to some other epoch? (The previous standard reference date used to be J1950, as far as I know).

Based on the lunar model in IAU documentation there seems to be periodic changes in Obliquity and LAN. Obliquity will change between 0.0273 and 0.0265 in periods of 16 days (approximately) and the LAN will drift about +-1.0 degrees from the mean value.

To be honest, I am inclined not to put support for higher-order axis perturbations into the core code. For that, I think it is better to provide an API interface so that it can be coded into the planet dll.

 

[jarmonik]

Does that mean that the discrepancy could arise from the fact that the AGC values may not be referenced to the J2000 ecliptic, but to some other epoch? (The previous standard reference date used to be J1950, as far as I know).

Yes, that's exactly what I was thinking.

This is how the reference frame is defined in AGC:

The orientation of this coordinate system is defined by the line of intersection of the mean earth equatorial plane and the mean orbit of the earth (the ecliptic ) at the beginning of the Besselian year which starts January 1.2516251, 1972, et. The X-axis is along this intersection with the positive sense in the direction of the ascending node of the ecliptic on the equator (the equinox), the Z-axis is along the mean earth north pole, and the Y-axis completes the right-handed triad.

But I don't know how AGC is implemented in NASSP or the Orbiter and I don't know how this problem is solved or is it solved at all. But there's probably no need to worry about it now.

PS. 01-Jul-1986 is an other reference date that's used in AGC.

 

[Chode]

I've been trying to model the precession of Triton with the Orbiter precession model, comparing it to data generated by the IAU document and have had good success, if I use L_rel insted of L_ecl in equation 7.

I can get the Orbiter model to agree with the IAU model to within 0.09 degrees for at least the next 1000 years (yes!). This is a good test, since the the three planes (ecliptic, reference and precession) are at large angles to one another for Triton. I've done some work towards figuring out how to easily translate the IAU data into the Orbiter model data, the method I used for Triton was somewhat cumbersome.

Regards

 

[jarmonik]

Here are my recommendations for the Moon.

SidRotPeriod = 2360591.597
SidRotOffset = 4.770748
Obliquity = 0.02694055
LAN = 2.181328
PrecessionPeriod = -6792.9076
LAN_MJD = 51544.5

This should fit in the IAU Lunar model.

 

[Chode]

I made a spreadsheet that carries out both the IAU calculation and the
Orbiter calculation of precession so I can directly compare the two.
I've made some plots of the angle between the rotation axes for the two
models.

The first plot shows the difference between the two models using the Moon.cfg values from jarmonik. (degrees vs. years)

For comparison, the second plot is the result if I set the precession period to infinity (no precession in Orbiter).

I've tried some fitting to the IAU data myself, and I concur with jarmoniks numbers as being about as good as your going to get

 

[jarmonik]

The IAU document linked to by BrianJ above seems pretty definitive. I'd be happy to do some work deriving the obliquity precession numbers from that. It does not have data on the Laplace planes though - if anyone can point me to a good source, I'm happy to do some number crunching to convert between coordinate systems, if required.

I am not familiar with concept of Laplace plane but could we just compute mean value of axis of rotation over a precession cycle, shouldn't that be at least pretty close to nodal precession pole.