Rotations

[dtamayo]

Hi Hanno, @tigerchenlu98, @shadden

Playing around with tides_spin (and debugging my own issues!) I think we do need to add some functionality for doing rotations in order for people to be able to use it easily. Obviously none of this is specific to spins, so I think it makes sense for these general functions to be in REBOUND. Not looking to merge now, but wanted to give you working code so you could try it out.

Probably the best way to get a feel for what might be a good API is to try out different things one might to track/plot in the examples. I love all the great examples Tiger's put together--I thought I'd add an even simpler example to more slowly go through some of the issues someone new would need to think through carefully. I think the issues with what might be the best API might be best resolved by looking through that example and how I'm calculating things. You can find it on the newpr91 branch of REBOUNDx (I've also sent Tiger a pull request to add to the main pull request).

https://github.com/dtamayo/reboundx/blob/newpr91/ipython_examples/SpinsIntro.ipynb

Here are my main questions:

• I like the idea of adding this functionality to the reb_vec3d class. Then it's really easy to do things we want to do a ton, like calculate the magnitude and spherical angles, e.g. vec.theta, vec.mag, or get it as a unit vector (vec.hat). We can then also have methods to rotate the vector into different reference planes (see example).

The main question I have is whether we want properties and methods to return vec3ds or lists. For example I've added a hvec and evec attribute to reb_orbit (both as reb_vec3ds) so we can easily get the orbit normal and eccentricity vectors. On the one hand it's nice if these are returned as reb_vec3ds if we want to subsequently want to rotate them using the built-in methods. But we're typically also trying to store outputs from calls to sim.integrate in lists. So we want to say something like h[i] = ps[1].hvec, but this doesn't work if hvec returns a reb_vec3d.

• For now I have the Python properties for hvec and evec returning a list rather than reb_vec3d (since that's also what we get with sim.calculate_angular_momentum()). If we want to rotate etc., then we just have to turn that list into a reb_vec3d. Anyway, it probably helps to see the example to play around with different syntaxes.

• Another option would be to not add this functionality to reb_vec3d and just have separate rotation functions (that take a list and return a list). I have some possible syntax for that option commented out in the example notebook.

• I've generalized align_simulation to take an arbitrary normal vector (so if you want the invariable plane great, but can also now do others). Once we decide on a syntax, I might add a rotations.ipynb notebook to REBOUND to show the basic functionality, e.g. rotating between the invariable plane and ecliptic. I also split it up between the simulation and rebx parameters. I think that's nice since I think it's a useful function to have in REBOUND by itself. One pitfall I ran into though is that if in the notebook you by mistake did, e.g.,

sim.rotate_simulation(normalvec=rebx.calculate_total_angular_momentum())
rebx.rotate_params(normalvec=rebx.calculate_total_angular_momentum())

you would get the wrong behavior, since you'd have rotated the orbits in the first call, and when you called rebx.calculate_total_angular_momentum() the second time, you would get a different answer (with rotated orbits but unrotated spins). Not sure what the best approach is there.

• I'm also not sure what the most useful way to plot these things is at the bottom. I like the idea of doing 3d arrow plots. The most effective visualizations of this I've seen in talks in movies of how these arrows move...perhaps there's a better static representation than what I made!

[dtamayo]

Definitely! I think the clearest thing to say in the documentation is just what you said: rot * [1,0,0] maps [1,0,0] into the orbital plane, i.e., it gives a vector pointing toward pericenter in the reference xyz coordinates.

[hannorein]

Ok. Making that change and updating the tests...

Amazingly long discussion for what can literally be written in 4 lines of code 😉:

struct reb_rotation P1 = reb_rotation_init_angle_axis(omega, (struct reb_vec3d){.z=1});
struct reb_rotation P2 = reb_rotation_init_angle_axis(inc, (struct reb_vec3d){.x=1});
struct reb_rotation P3 = reb_rotation_init_angle_axis(Omega, (struct reb_vec3d){.z=1});
return reb_rotation_mul(P3, reb_rotation_mul(P2, P1));

[dtamayo]

Haha, just means it's well written code, and that this is confusing! I personally don't want to figure this out again :)

[dtamayo]

Oops, I realized that I didn't update to_new_axes to make newx perpendicular. Let me know when you're done and I'll add that!

[hannorein]

Yes, please add that!

[dtamayo]

(oops, fixing now)

[hannorein]

Amazing, I think this means we're ready to merge!

[dtamayo]

Maybe you already fixed it, but I pushed the syntax fix and updated the unit test for orbital to match.

Thank you!!

[hannorein]

Yes. Thank you!