Possible problem with rotation equations in forces.c

[Larry1019]

I think I see a problem with the rotation equations in the J2,J3,J4 subroutine -- lines 492-495 in forces.c.

For starters, assume the rotations were exactly zero in both ra and dec (ie no rotation at all). I would think that you would expect that the rotation equations would yield: dxp = dx dyp = dy dzp = dz

But if you set a=0 and d=0 with cosa=1, sina=0, cosd=1, sind=0 the equations give:

dxp = dy dyp = dz dzp = dx

which CAN'T be right and implies that the rotation equations are wrong.

Second, you're taking the value of a and d to be the actual RA and Dec of the Earth's pole (at J2000), ie 0 and 90 respectively and not the rotation angle of the pole from the J2000 position, both of which should both be very small angles and are actually both zero if you want to fix things at J2000. But, you're effectively saying you want to rotate the pole through 90 degrees in Dec and 0 degrees in RA.

And a very minor point, you calculate costheta2 = dz*dz/r2. But it's actually the sine of an angle, not a cosine. And I think it's actually the sine of the "dec" of the object (but calculated from the non-light time corrected positions).

Larry Wasserman Lowell Observatory lhw@lowell.edu

[hannorein]

I'll let Matt comment on this, just one small thing: note that an inverse rotation is performed a few lines further down.

[Larry1019]

I was just looking at this. You apparently are rotating x,y,z geocentric equatorial coordinates to some other coordinate frame. Davide told me that it should be the coordinates defined by the moving position of the polar axis of the earth with time so that the rotation should be through very small angles in ra and dec unless you;re integrating over thousands of years. The accelerations are computed using the rotated x,y,z in that frame and then the accelerations are rotated back to equatorial coordinates. So I would expect that the acceleration rotations should be through -a, -d to get back to an equatorial coordinate frame if the first rotations for x,y,z were through a, d. But the equations for the accelerations clearly have a different structure than the rotations for the coordinates. If this is not the case, I'd really like an explanation of what is happening here.

Larry Wasserman

On Sun, Jul 16, 2023 at 4:29 PM Hanno Rein @.***> wrote:

I'll let Matt comment on this, just one small thing: note that an inverse rotation is performed a few lines further down https://github.com/matthewholman/assist/blob/main/src/forces.c#L534C9-L534C14 .

— Reply to this email directly, view it on GitHub https://github.com/matthewholman/assist/issues/86#issuecomment-1637233280, or unsubscribe https://github.com/notifications/unsubscribe-auth/BBKAAVMK4FPPSEOHHG4IDW3XQSBPLANCNFSM6AAAAAA2MH4RKA . You are receiving this because you authored the thread.Message ID: @.***>

[Larry1019]

Upon thinking about it, I'm not sure if my reply earlier this evening is correct. So I will retract that message until we sort out my original post which asks (1) Why if the two rotations ( a and d) are zero is the output x,y,z not identical to the input x,y,z? and (2) Why if you are correcting for the change in the obliquity of the axis of the Earth over time (a small amount) are you rotating through 90 degrees in Dec?

Larry

On Sun, Jul 16, 2023 at 5:29 PM Hanno Rein @.***> wrote:

I'll let Matt comment on this, just one small thing: note that an inverse rotation is performed a few lines further down https://github.com/matthewholman/assist/blob/main/src/forces.c#L534C9-L534C14 .

— Reply to this email directly, view it on GitHub https://github.com/matthewholman/assist/issues/86#issuecomment-1637233280, or unsubscribe https://github.com/notifications/unsubscribe-auth/BBKAAVMK4FPPSEOHHG4IDW3XQSBPLANCNFSM6AAAAAA2MH4RKA . You are receiving this because you authored the thread.Message ID: @.***>