Moon.cfg in orbiter beta

[tblaxland]

Just one more thought - is it the case that for planets, the precessional period is in the order of 10,000's to 100,000's of years? If so, is it really worth bothering about? Is the Moon a special case? (in respect of it's short precessional period and importance as a target for historical simulations)

I don't think it would be worth worrying about if it was just for the planets. Yes, the Moon is a special case because of its historical value but not because of its precession period. I would be surprised if its short precession period was unique, although I don't have any data to hand. I expect that any smallish body formed by an energetic collision would display significant precession. EDIT: If you are going to model the precession of the Moon, it being a special case, I don't why such a model should not be made use of for other bodies (asteroids was what I had in mind though it would be purely hypothetical since I don't think there is any real data available).

 

[jarmonik]

I have been trying to search more information from about precession rates of other bodies without having much success. But I did encounter some interesting documents. We could have a good chances to compute approximate precession rates for many other bodies.

This http://www.es.ucsc.edu/~fnimmo/website/Bills_Titan.pdf could be interesting in a hands of someone who understand these things.

Here is an other one http://www.gtti.it/GTTI08/papers/persi.pdf

 

[BrianJ]

I don't think it would be worth worrying about if it was just for the planets. Yes, the Moon is a special case because of its historical value but not because of its precession period. I would be surprised if its short precession period was unique, although I don't have any data to hand. I expect that any smallish body formed by an energetic collision would display significant precession. EDIT: If you are going to model the precession of the Moon, it being a special case, I don't why such a model should not be made use of for other bodies (asteroids was what I had in mind though it would be purely hypothetical since I don't think there is any real data available).

I agree

I have been trying to search more information from about precession rates of other bodies without having much success.

This paper "Report of the IAU/IAGWorking Group on cartographic coordinates and rotational elements: 2006" is quite comprehensive.
http://www.dtic.mil/cgi-bin/GetTRDoc?AD=ADA471037&Location=U2&doc=GetTRDoc.pdf


Poor Martins! All he wanted to know was whether it would be a problem for Apollo add-on developers if he changed the moon.cfg parameters a bit!!

Since mission-specific moon.cfg's are easy to include, I don't see why it should be.

 

[jarmonik]

This paper "Report of the IAU/IAGWorking Group on cartographic coordinates and rotational elements: 2006" is quite comprehensive.
http://www.dtic.mil/cgi-bin/GetTRDoc?AD=ADA471037&Location=U2&doc=GetTRDoc.pdf

That's a great finding. :speakcool: I suppose this would require implementing the rotational elements into a planet modules.

 

[tblaxland]

This http://www.es.ucsc.edu/~fnimmo/website/Bills_Titan.pdf could be interesting in a hands of someone who understand these things.

Here is an other one http://www.gtti.it/GTTI08/papers/persi.pdf

This paper "Report of the IAU/IAGWorking Group on cartographic coordinates and rotational elements: 2006" is quite comprehensive.
http://www.dtic.mil/cgi-bin/GetTRDoc?AD=ADA471037&Location=U2&doc=GetTRDoc.pdf

Thanks guys, it looks like I've found my lunch time reading for today

That's a great finding. :speakcool: I suppose this would require implementing the rotational elements into a planet modules.

Interesting suggestion. I had suggested earlier incorporating the precession components into a CELBODY2 class. I guess if you are going to do that, why not go the whole hog and include rotation too. Of course, it would be nice to maintain the ability to create simple planets from just cfg files.

---------- Post added at 13:28 ---------- Previous post was at 09:04 ----------

This http://www.es.ucsc.edu/~fnimmo/website/Bills_Titan.pdf could be interesting in a hands of someone who understand these things.

That document implies that Cassini states are common amongst natural satellites - common enough that they feel justified to assume that Titan has settled into such a state around Saturn. There is some discussion of that and they mention timescales of 1e6 years for such damping to occur.

Here is an other one http://www.gtti.it/GTTI08/papers/persi.pdf

This paper confirms the assumption made by the one above.

This paper "Report of the IAU/IAGWorking Group on cartographic coordinates and rotational elements: 2006" is quite comprehensive.
http://www.dtic.mil/cgi-bin/GetTRDoc?AD=ADA471037&Location=U2&doc=GetTRDoc.pdf

Awesome find! Just glancing through the data I calculate the Moons axis to precess with a period of 18.59949 Julian years, so that agrees with the data presented elsewhere (such as at JPL: http://ssd.jpl.nasa.gov/?sat_elem). Using that data also indicates that Phobos and Deimos are also in Cassini states adding further weight to the argument that such states are wide spread for natural satellites. Looking further out shows Europa not in a Cassini state with an orbital node precession period of 1.394 years and an axis precession period of 30.219 years. Other Jupiter satellites are not in Cassini states but they tend to have quite small amounts of precession anyway (Europa has the most of the Galliean satellites with +/- 1.086 deg, Thebe the most of the non-Galliean satellites with +/- 2.11 deg).

On that basis, Martin's proposal of defining a natural satellite's spin axis precession in the rotating frame of its orbit may be acceptable (see my notes on the IAU Working Group report below). Just to note that in such a scheme, the LAN of obliquity can only ever either 0 or 180 when measured in the orbit frame. Other values would be unstable, not Cassini states and would have the body rotating in a physically unlikely way. That said, since a Cassini state can also be modelled by a more generic precession model (by aligning or anti-aligning the LANs and setting the periods the same) I think the more generic precession model would be better, provided said model is not too difficult to implement. Your call I guess Martin, since you are writing the code

 

[martins]

Working on this, although I'm a bit short on time at the moment. I've decided not to define the precession in the rotating frame of the orbital plane after all, but with respect to the ecliptic, as suggested above. I think I've sorted out the necessary transformation. Just need to code it now.

For now, I'll provide cfg file options to define the precession parameters. A more generic API interface could follow later. The most general version would simply return a rotation matrix to transform the planet from local to global coordinates. Then everything (rotation, precession, etc.) would be up to the module.

 

[jarmonik]

For now, I'll provide cfg file options to define the precession parameters. A more generic API interface could follow later. The most general version would simply return a rotation matrix to transform the planet from local to global coordinates. Then everything (rotation, precession, etc.) would be up to the module.

That's a great news.:speakcool: If I understood correctly, precession parameters would also include LAN and Inclination of Laplace plane respect to ecliptic. (i.e. precession pole). But how would the Orbiter define rotation period ? I suppose a period from LAN to LAN would make sense.

 

[tblaxland]

Working on this, although I'm a bit short on time at the moment. I've decided not to define the precession in the rotating frame of the orbital plane after all, but with respect to the ecliptic, as suggested above. I think I've sorted out the necessary transformation. Just need to code it now.

For now, I'll provide cfg file options to define the precession parameters. A more generic API interface could follow later. The most general version would simply return a rotation matrix to transform the planet from local to global coordinates. Then everything (rotation, precession, etc.) would be up to the module.

Sounds great :speakcool:. I haven't reviewed the graphic you posted because it keeps crashing my browser. I'll try again when I'm on a different machine.

But how would the Orbiter define rotation period ? I suppose a period from LAN to LAN would make sense.

I agree.

---------- Post added at 19:27 ---------- Previous post was at 19:15 ----------

Sounds great :speakcool:. I haven't reviewed the graphic you posted because it keeps crashing my browser. I'll try again when I'm on a different machine.

OK, got it working now. That looks pretty good and the same transformation can be used for the Laplace plane and the equatorial plane. EDIT: I'm curious what program you have done that animation in?

 

[martins]

EDIT: I'm curious what program you have done that animation in?

The frames were generated by Matlab, and then baked together into an animated GIF with a freeware utility called Ulead GIF animator LE. If you are interested, here is the matlab script (ugly as hell, sorry, wasn't intended for public consumption ):

function fig1

clear all
close all

RAD=pi/180.0;
DEG=180.0/pi;

obliq_ref = 30*RAD;
lan_ref = 40*RAD;
obliq_axis = 19.5*RAD;
lan_axis = 320*RAD;

rad = 1.5*tan(obliq_axis);

Oref = [[1 0 0];[0 cos(obliq_ref) -sin(obliq_ref)];[0 sin(obliq_ref) cos(obliq_ref)]];
Lref = [[cos(lan_ref) -sin(lan_ref) 0];[sin(lan_ref) cos(lan_ref) 0];[0 0 1]];
Rref = Lref*Oref;

Oaxis = [[cos(obliq_axis) 0 sin(obliq_axis)];[0 1 0];[-sin(obliq_axis) 0 cos(obliq_axis)]];
Laxis = [[cos(lan_axis) -sin(lan_axis) 0];[sin(lan_axis) cos(lan_axis) 0];[0 0 1]];
Raxis = Laxis*Oaxis;
ref_axis = Rref * [0 0 1.5]';

figure;
set(gcf,'Position',[0 0 900 450]);
subplot(1,2,1);
line([-1 1],[0 0],[0 0]);
line([0 0],[-1 1],[0 0]);
line([0 0],[0 0],[0 1.5]);
patch([-1 -1 1 1],[-1 1 1 -1],[0 0 0 0],[.7 .7 .7],'FaceAlpha',0.2);
text(1.05,0,0,'vernal equinox');
text(0,0,1.6,'ecliptic north pole','HorizontalAlignment','right');
text(-0.9,1.4,0,'ecliptic');
line([0 ref_axis(1)],[0 ref_axis(2)],[0 ref_axis(3)],'LineWidth',2);
line([ref_axis(1) ref_axis(1)],[ref_axis(2) ref_axis(2)],[0 ref_axis(3)]);
line([0 ref_axis(1)],[0 ref_axis(2)],[0 0]);
c=Lref*[1.3 0 0]';
line([c(1) -c(1)],[c(2) -c(2)],[0 0],'LineStyle','--');
text(c(1)*1.05,c(2)*1.05,0,'AN_{ref}');

c = circle(rad);
c(:,3) = c(:,3)+1.5;
for i=1:size(c,1), c(i,:) = Rref*c(i,:)'; end
line(c(:,1),c(:,2),c(:,3),'Color',[.3 .3 .3]);

c = arc(0.7,pi,pi+obliq_ref);
for i=1:size(c,1), c(i,:) = Lref * [[0 0 1];[0 1 0];[-1 0 0]] * c(i,:)'; end
line(c(:,1),c(:,2),c(:,3));
c = arc(0.6,0,lan_ref);
line(c(:,1),c(:,2),c(:,3));
text(0.05,0.0,0.5,'\omega_{ref}');
text(0.7,0.3,0,'L_{ref}');

text(-1,0,2.3, ['\omega_{ref} = ' num2str(obliq_ref*DEG) '\circ'], 'Color', 'blue');
text(-1,0,2.1, ['L_{ref} = ' num2str(lan_ref*DEG) '\circ'], 'Color', 'blue');

axis equal tight
axis off
view(-12,16)
set(gca,'Projection','perspective');
%set(gca,'Position',[0.0 0.05 0.5 1])

lanrel = [0:2*pi/200:2*pi];
for i=1:length(lanrel)
    [obl(i) lan(i)] = axis_ecl(obliq_ref,lan_ref,obliq_axis,lanrel(i));
end
subplot(2,2,2);
plot (lanrel,obl);
hold on
xlabel('LAN_{rel}');
ylabel('obliquity_{ecl}');
axis([0 2*pi 0.1 0.9]);
%set(gca,'Position',[0.65 0.6010 0.33 0.324]);
subplot(2,2,4);
plot (lanrel,lan);
hold on
xlabel('LAN_{rel}');
ylabel('LAN_{ecl}');
axis([0 2*pi -0.1 1.5]);
%set(gca,'Position',[0.65 0.1100 0.33 0.324]);


h = [];
pause
for dphi = 0:2:360
    lan_axis = RAD*dphi;
    Laxis = [[cos(lan_axis) -sin(lan_axis) 0];[sin(lan_axis) cos(lan_axis) 0];[0 0 1]];

    [obliq_axis_tot lan_axis_tot] = axis_ecl(obliq_ref, lan_ref, obliq_axis, lan_axis);
    Oaxistot = [[1 0 0];[0 cos(obliq_axis_tot) -sin(obliq_axis_tot)];[0 sin(obliq_axis_tot) cos(obliq_axis_tot)]];
    Laxistot = [[cos(lan_axis_tot) -sin(lan_axis_tot) 0];[sin(lan_axis_tot) cos(lan_axis_tot) 0];[0 0 1]];
    Raxistot = Laxistot*Oaxistot;
    rot_axis = Raxistot * [0 0 norm([rad 0 1.5])]';

    subplot(1,2,1);
    delete(h);
    h(1) = line([0 rot_axis(1)],[0 rot_axis(2)],[0 rot_axis(3)],'LineWidth',2,'Color','red');
    h(2) = line([rot_axis(1) rot_axis(1)],[rot_axis(2) rot_axis(2)],[0 rot_axis(3)],'Color','red');
    h(3) = line([0 rot_axis(1)],[0 rot_axis(2)],[0 0],'Color','red');
    h(4) = line([ref_axis(1) rot_axis(1)],[ref_axis(2) rot_axis(2)],[ref_axis(3) rot_axis(3)],'LineWidth',2,'Color',[0 0.5 0]);

    c = [rad 0 0]';
    c(3) = 1.5;
    c = [Rref*c,Rref*[-c(1);-c(2);c(3)]];
    h(5) = line([c(1,1) c(1,2)],[c(2,1) c(2,2)],[c(3,1) c(3,2)]);
    
    c = Laxis*[rad 0 0]';
    c(3) = 1.5;
    c = [Rref*c,Rref*[-c(1);-c(2);c(3)]];
    h(6) = line([c(1,1) c(1,2)],[c(2,1) c(2,2)],[c(3,1) c(3,2)],'LineStyle','--','Color',[0 0.5 0]);
    h(7) = text(c(1,1)*1.1, c(2,1)*1.1, c(3,1),'AN_{rel}','Color',[0 0.5 0]);
    
    c=Laxistot*[1.3 0 0]';
    h(8) = line([c(1) -c(1)],[c(2) -c(2)],[0 0],'Color','red','LineStyle','--');
    h(9) = text(c(1)*1.05,c(2)*1.05,0,'AN_{axis}','Color','red');

    c = arc(0.7,pi,pi+obliq_axis_tot);
    for i=1:size(c,1), c(i,:) = Laxistot * [[0 0 1];[0 1 0];[-1 0 0]] * c(i,:)'; end
    h(10) = line(c(:,1),c(:,2),c(:,3),'Color','red');
    h(11) = text(0.15,0.0,0.75,'\omega_{ref}','Color','red');

    c = arc(0.7,0,lan_axis_tot);
    nh = 0; % size(c,1);
    h(12) = line(c(:,1),c(:,2),c(:,3),'Color','red');

    c = arc(0.3,0,lan_axis);
    c(:,3) = c(:,3)+1.5;
    for i=1:size(c,1) c(i,:) = Rref * c(i,:)'; end
    h(13) = line(c(:,1),c(:,2),c(:,3),'Color',[0 0.5 0]);

    h(14) = text(-0.2,0,2.25, ['\omega_{rel} = ' num2str(obliq_axis*DEG) '\circ'], 'Color', [0 0.5 0]);
    h(15) = text(-0.2,0,2.05, ['L_{rel} = ' num2str(lan_axis*DEG) '\circ'],'Color',[0 0.5 0]);
    h(16) = text( 0.6,0,2.2,  ['\omega_{tot} = ' num2str(obliq_axis_tot*DEG) '\circ'], 'Color','red');
    h(17) = text( 0.6,0,2.0,  ['L_{tot} = ' num2str(lan_axis_tot*DEG) '\circ'], 'Color','red');
    
    subplot(2,2,2);
    h(18+nh) = plot(lan_axis,obliq_axis_tot,'ro','MarkerSize',10);
    subplot(2,2,4);
    h(19+nh) = plot(lan_axis,lan_axis_tot,'ro','MarkerSize',10);

    drawnow
    frame = getframe(gcf);
    img = frame2im(frame);
    name = sprintf ('img_%03d.bmp',dphi);
    imwrite(img,name,'bmp');
    pause(0.01);
end

%movie2avi(frame,'precession.avi','compression','MSVC');
end



function vtx=circle(rad)

nseg=64;
for i=0:nseg
    phi = i/nseg*pi*2;
    vtx(i+1,:) = [rad*cos(phi) rad*sin(phi) 0];
end

end


function vtx=arc(rad,phi0,phi1)

nseg=64;
vtx = [0 0 0];
nstep = ceil((phi1-phi0)/(2*pi)*nseg);
if nstep > 0
    dphi=(phi1-phi0)/nstep;
    i=1;
    for phi=phi0:dphi:phi1
        vtx(i,:) = [rad*cos(phi) rad*sin(phi) 0];
        i = i+1;
    end
end

end


% Calculate obliquity and LAN of axis w.r.t. ecliptic
% from obliquity and LAN of precession reference axis
% and relative obliquity and LAN
function [obliq lan]=axis_ecl(obliq_ref,lan_ref,obliq_rel,lan_rel)

Oref = [[1 0 0];[0 cos(obliq_ref) -sin(obliq_ref)];[0 sin(obliq_ref) cos(obliq_ref)]];
Lref = [[cos(lan_ref) -sin(lan_ref) 0];[sin(lan_ref) cos(lan_ref) 0];[0 0 1]];
Rref = Lref*Oref;
Orel = [[1 0 0];[0 cos(obliq_rel) -sin(obliq_rel)];[0 sin(obliq_rel) cos(obliq_rel)]];
Lrel = [[cos(lan_rel) -sin(lan_rel) 0];[sin(lan_rel) cos(lan_rel) 0];[0 0 1]];
Rrel = Lrel*Orel;
v=[0 0 1];

rot_axis = Rref * Rrel * v';
obliq = acos(rot_axis(3));
lan = atan2(rot_axis(1),-rot_axis(2));
end

Btw, what is the browser that crashed on the animation? I thought any browser would be capable of showing animated gifs.

 

[martins]

That's a great news.:speakcool: If I understood correctly, precession parameters would also include LAN and Inclination of Laplace plane respect to ecliptic. (i.e. precession pole). But how would the Orbiter define rotation period ? I suppose a period from LAN to LAN would make sense.

Do you mean the period of rotation of the axis, or the period of rotation of the planet itself? The latter may need a bit of thought, because at the moment the transformations are cumulative. So this would add one full planet rotation to each precession cycle. If I keep the siderial rotation periods to calculate angular velocity, the moon would show its back side after 9.3 years. Is the solution as simple as subtracting the precession angular velocity from the rotation angular velocity?

 

[jarmonik]

Do you mean the period of rotation of the axis, or the period of rotation of the planet itself?

Yes, I was talking about the planet itself.

The latter may need a bit of thought, because at the moment the transformations are cumulative. So this would add one full planet rotation to each precession cycle. If I keep the siderial rotation periods to calculate angular velocity, the moon would show its back side after 9.3 years. Is the solution as simple as subtracting the precession angular velocity from the rotation angular velocity?

It could be that simple.

I was just thinkin that defining a sidereal rotation period respect to fixed star in a case when the rotation axis will change isn't actually so self-evident. Defination of angular velocity relative to LAN would make more sense.

But, since, the Orbiter is allready using the sidereal rotation period it's probably good idea to maintain that and the compatibility with current configuration files.

Also I had in mind that should one multibly the precession velocity with Cos(Obliquity) before subtracting it. For the Moon that would make a difference of 4 seconds.

 

[martins]

So this would add one full planet rotation to each precession cycle. If I keep the siderial rotation periods to calculate angular velocity, the moon would show its back side after 9.3 years. Is the solution as simple as subtracting the precession angular velocity from the rotation angular velocity?

Actually, scratch that. I realised that this is only true if the precession is around the ecliptic, otherwise (like in the animation), it just means that absolute angular velocity is wobbling around the mean value. Probably the better strategy is to subtract the current LAN from the SidRotOffset value. What do you think?

Another point: I had been thinking about applying the precession calculation only once at the start of the simulation, but I am beginning to think that this would be a bad idea, because it could introduce inconsistencies (lauching a saved simulation would not return to the same physical state).

If I update the precession axis in-game, then I can either do it continuously at each frame, or (to save a few processor cycles), periodically at longer intervals. The update frequency would have to be chosen individually for each body so that no significant jumps are introduced (e.g. as seen from a vessel hovering close to the surface).

 

[martins]

Defination of angular velocity relative to LAN would make more sense.

Re-reading your post, I realised that this is a point I don't quite understand. The motion of the LAN is not necessarily uniform (only if the precession reference axis coincides with the ecliptic normal), so it seems to me that defining the angular velocity with respect to LAN would not be a good idea. Am I missing something?

 

[tblaxland]

The frames were generated by Matlab, and then baked together into an animated GIF with a freeware utility called Ulead GIF animator LE. If you are interested, here is the matlab script (ugly as hell, sorry, wasn't intended for public consumption ):

Thanks for the info. One of these days, I really should learn Matlab or one of its alternatives.

Btw, what is the browser that crashed on the animation? I thought any browser would be capable of showing animated gifs.

Firefox 3.0.10. I think it was some sort of memory issue since the crash could be prevented if I stopped the browser before it loaded all the frames. I saved the link and viewed it in Irfanview no problem.

Probably the better strategy is to subtract the current LAN from the SidRotOffset value. What do you think?

That sounds right to me.

If I update the precession axis in-game, then I can either do it continuously at each frame, or (to save a few processor cycles), periodically at longer intervals. The update frequency would have to be chosen individually for each body so that no significant jumps are introduced (e.g. as seen from a vessel hovering close to the surface).

Another option would be to dynamically update the precession calculation frequency based on the camera distance to the object and the angular velocity of precession. I'm not sure if that would really save many cycles overall though.... Otherwise, a manually specified update frequency in the config file would be fine.

 

[jarmonik]

Originally Posted by martins
Re-reading your post, I realised that this is a point I don't quite understand. The motion of the LAN is not necessarily uniform (only if the precession reference axis coincides with the ecliptic normal), so it seems to me that defining the angular velocity with respect to LAN would not be a good idea. Am I missing something?

I was meaning the LAN of the laplace plane and the equatorial plane of a planet. That's the (AN_rel) in the animation.

Originally Posted by martins
Probably the better strategy is to subtract the current LAN from the SidRotOffset value. What do you think?

That's a good idea as long as the SidRotOffset is the angle between AN_rel and meridian in the epoch. If not then I am not so sure about it.

Also, in that case if the paramaters L_ref, omega_ref and LAN-rate are not defined in configuration files, the Orbiter could assume them to be zero. Then, the Laplace plane would be equal to the ecliptic and old configuration files would work with the new model same way as they are working right now. Right ?

 

[martins]

Note: Public beta 090520 has been released with support for axis precession. Please have a look and let me know if this looks reasonable. Currently I've inserted rather preliminary precession data for Earth and Moon. Any improvements on these, and data for any other bodies, are most welcome. I would hope that we can come up with definitive rotation/precession data, so that major changes in the future will no longer be necessary.

 

[tblaxland]

I would hope that we can come up with definitive rotation/precession data, so that major changes in the future will no longer be necessary.

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.

---------- Post added at 21:25 ---------- Previous post was at 17:58 ----------

Further to my last post, since we are working on a mean value for precession rate and a static direction for the precession axis, I have some questions:

1. What direction should I take for the axis of precession? The value at J2000 or some other value?

2. The rate of precession is not constant, so I am thinking we should take a mean value. Over what period should I take the mean? J2000 +/- 100 years? 4 Oct 1957 to 4 Oct 2057? Longer? shorter?

Anyone have an opinion?

 

[martins]

1. What direction should I take for the axis of precession? The value at J2000 or some other value?

Yes, J2000 should be the invariant reference frame.

2. The rate of precession is not constant, so I am thinking we should take a mean value. Over what period should I take the mean? J2000 +/- 100 years? 4 Oct 1957 to 4 Oct 2057? Longer? shorter?

I guess J2000 +/- 100 years is a good choice, since probably most scenarios fall within this time frame. How fast does the rate diverge?

 

[jarmonik]

I have compared current rotation model to AGC model and there are some problems.

1. Since, the L_{rel} is substracted from SidRotOffset (in eq. 7), the L_0 must be added in SidRotOffset values in configuration files to compensate that. The problem is that L_{rel} includes L_0. I suppose this could be sovled by replacing L_{rel} with 2PI(t-t0)/T_p in eq. 7

2. Current rotation model seems to drift a little form the AGC model, but the problem is that witch one is more correct ?

AGC computes the rotation form r(t) = L + dL(t-t0) + F + dF(t-t0),

where dF = 2PI / 2351139.374 (angular velocity relative to LAN)


Orbit period form node to node can be computed from

T_nn = (1/T_s - cos(Obliquity)/T_p)^{-1}

Using SidRotPeriod T_s = 2360591.672 this will give 2351139.359 for T_nn witch is pretty close to the one used by AGC. Where as (1/T_s - 1/T_p)^{-1} would give 2351135.98

This will lead into a suggestion...

Any thoughts ?

There is an other explanation for cos(Obliquity) but I'll need to sketch a drawing.

---------- Post added at 11:02 PM ---------- Previous post was at 10:56 PM ----------

Of course, it would be possible to get an other reference system by implementing lunar rotation model from http://www.dtic.mil/cgi-bin/GetTRDoc?AD=ADA471037&Location=U2&doc=GetTRDoc.pdf into a MFD.

---------- Post added at 01:01 AM ---------- Previous post was Yesterday at 11:02 PM ----------

Here's the sketch.

'L' is the longitude of ascending node and 'R' is the rotation offset. Now, if the 'L' is modified by nodal precession 'dL' then the rotation offset must be increased to keep the meridian in the same place. However, 'dR' is not equal to 'dL'. From the law of right angled spherical triangles: sin(dR) = sin(dL)*cos(Obliquity). This could be approximated as dR = dL*cos(Obliquity), at least my pocket calculator doesn't make a difference. Taking the cos(Obliquity) in account that would improve stability by a factor of 1000 (compared to AGC model). Current deviation is about 0.2 degrees within 32 years. It would be reduced to 0.002 degrees.

AGC model can be viewed with LTMFD by using the NASSP program. It does require this patch http://koti.mbnet.fi/jarmonik/SSE2Test.zip Debug programs are disabled from official build.

 

[martins]

Ok, I think I see where you are coming from. So just to clarify: you think that the last term in Eq. 7 of precession.pdf should be

L_{ecl} * cos(\eps_{rel})

instead of just L_{ecl} ? I can well belive that this is true. I'll give it some more thought tonight ...

I guess your expression for \phi in your post above is only valid for the case where the precession reference coincides with the ecliptic normal, since the last term is just a linear function in t, so in the general case I couldn't replace L_{ecl} with 2\pi(t-t_0)/T_p.