bad bug in hotspot offset

[rodluger]

If I add a latitudinal offset Phi = 30 degrees to a planet in an orbit with inclination = 80 degrees, I get something that looks good:

But when the inclination is close to zero, I get the following:

I think my math in Appendix B.3 may be wrong: specifying an offset in latitude at low inclination leads to an offset in longitude, with some other weird results.

To reproduce this bug, run

python -c "import planetplanet as pp; import matplotlib.pyplot as pl; pp.DrawOrbit(size = 2, inc = 5, Phi = 30, nz = 11); pl.show()"

[rodluger]

Equation (B11) is missing a minus sign in the expression for z_star (if we're in a left-handed coordinate system), but I don't think that's the full issue. I need to re-derive (B12). Here's the correct version of (B11) and the current version of (B12), which is where I think I've made a mistake:

[rodluger]

I think I know what the problem is. The Cartesian coordinates of the planet are not enough to apply these offsets in the general case, since I'm defining "longitude" along the orbital plane and "latitude" perpendicular to it. The problem is therefore underdetermined, as I need to know the velocity vector (3 additional variables) of the planet to define these directions.

Alternatively, I can go back to my original routine for doing this, which involved transformations by the orbital angles. I would start in the prime frame, as above, then rotate along the y axis by the true anomaly, along the x axis by the inclination and along the z axis by the longitude of ascending node.

@ericagol, what do you think?

[ericagol]

I think that sounds good.

Eric Agol Astronomy Professor University of Washington

On Aug 26, 2017, at 7:07 PM, Rodrigo Luger notifications@github.com wrote:

I think I know what the problem is. The Cartesian coordinates of the planet are not enough to apply these offsets in the general case, since I'm defining "longitude" along the orbital plane and "latitude" perpendicular to it. The problem is therefore underdetermined, as I need to know the velocity vector (3 additional variables) of the planet to define these directions.

Alternatively, I can go back to my original routine for doing this, which involved transformations by the orbital angles. I would start in the prime frame, as above, then rotate along the z axis by the true anomaly, along the x axis by the inclination and along the z axis by the longitude of ascending node.

@ericagol, what do you think?

― You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or mute the thread.

[rodluger]

OK, I figured it out! Here's the correct solution for an inclination of 5 degrees and a latitudinal offset of 30 degrees:

The approach using the orbital angles would have worked, but REBOUND doesn't automatically compute the true anomaly, so that would have been slow. I figured out how to use the position and velocity cartesian components to define the latitudinal and longitudinal directions. This is something @ericagol suggested a while back but I didn't quite know how to do.

The idea is to compute the position of the substellar point,

which defines a vector rstar. We then use the Rodrigues formula to rotate by Lambda along the orbital plane (defined by the normal vector r x v) or by Phi in a plane defined by the normal vector (r x v) x r.

Compound transformations (in both latitude and longitude) are a bit trickier. To keep things simple, I'm currently applying the rotations sequentially in longitude and then in latitude. I rotate by Lambda about the vector r x v to get rhotspot' then rotate by Phi about the vector (r x v) x rhotspot' to get the final rhotspot vector. Is this sensible?

I had to adapt the C code to return the array of velocity vectors, but it's very straightforward. I still need to implement this for the Kepler solver, since I don't know the equations to get the velocity from the orbital elements. I'll figure it out from Murray & Dermott this week.

These changes are still in the dev branch but I'll push them to master once I'm done. The function that computes the eyeball angles is here.

[rodluger]

The sky-projected velocity is this:

For reference, I derived this from expression (8) in this link and the linear transformation P3P2P1 on page 51 of Murray and Dermott.

[rodluger]

Done!