hannorein / rebound

💫 An open-source multi-purpose N-body code.
https://rebound.readthedocs.io/
GNU General Public License v3.0
846 stars 219 forks source link

Rotations warning? #671

Closed tigerchenlu98 closed 1 year ago

tigerchenlu98 commented 1 year ago

Hi Hanno,

This is very minor, but I thought I'd bring it to your attention. Creating a reb_rotation object to the same axes leads to NaNs in the rotation object. For instance,

struct reb_vec3d newz = reb_vec3d_add(reb_tools_angular_momentum(sim), rebx_tools_spin_angular_momentum(rebx)); struct reb_vec3d newx = reb_vec3d_cross((struct reb_vec3d){.z =1}, newz); struct reb_rotation rot = reb_rotation_init_to_new_axes(newz, newx);

returns NaNs if newz & newx are the same as the original z & x axes. Clearly this rotation does nothing - but maybe some sort of warning should be put in here?

Thanks! Tiger

hannorein commented 1 year ago

I'm not sure I understand. If the new z is the same as the old z, then the cross product

reb_vec3d_cross((struct reb_vec3d){.z =1}, newz);

is not well defined. So I think every function here is doing everything right.

tigerchenlu98 commented 1 year ago

Oh yes everything is absolutely working correctly. I'm just wondering if it would be good to add warnings in cases like this. But the user should know when it's appropriate to use the rotation so maybe not - up to you!

hannorein commented 1 year ago

Where does this code appear?

tigerchenlu98 commented 1 year ago

This is in one of the example folders in my REBOUNDx fork:

https://github.com/tigerchenlu98/reboundx/blob/chains/examples/hatp13/problem.c

the issue occurs when you try to shift a completely coplanar system to the invariant frame (which it's already in)

hannorein commented 1 year ago

I don't think it makes sense to add a warning message to any of the reb_ functions. Again, they do exactly what they are supposed to do. Returning NaN when the cross product is not defined is good. It's the same as if you divide by zero. I think you should add an exception to your side of the code. Maybe as simple as

if isnan(rot.r) {
   rot = reb_rotation_identity();
}

or

if (reb_vec3d_length_squared(newx)<some_value) {
   rot = reb_rotation_identity();
} 
tigerchenlu98 commented 1 year ago

That's fair! Closing the issue.