dtamayo / reboundx

A library for adding additional forces to the REBOUND N-body integration package
GNU General Public License v3.0
80 stars 60 forks source link

Test particle spin EoM bug #108

Closed tigerchenlu98 closed 1 year ago

tigerchenlu98 commented 1 year ago

Hey @dtamayo, found a bug: there's no check for test particles in the spin vector equation of motion, so a single test particle would introduce NaNs in all spin vector EoMs. This should fix that!

hannorein commented 1 year ago

Good catch!

dtamayo commented 1 year ago

Thanks @tigerchenlu98! Hope you had a good memorial day weekend. Does it matter that mu_ij will be nan when both masses are zero (perhaps too edge an edge case?)?

I haven't fully wrapped my head around all possibilities. If you set m=0 but moment of inertia I > 0, will the particle's spin vector evolve and the simulation doesn't conserve angular momentum? I can also look at this more carefully tomorrow morning.

hannorein commented 1 year ago

I think the right way would be to use the specific moment of inertia so that mu/I is well defined in the test particle case.

tigerchenlu98 commented 1 year ago

It is a (very) edge case that doesn't make much physical sense - the spin axis of a test particle being affected by tides. But I guess in principle I could see some cases where people may have use for it... I'll think about it more deeply tonight!

dtamayo commented 1 year ago

I think it's a reasonable limit that people doing analytical work might want to compare against N-body. You might for example want to model the spins of two earth mass planets going around a star, and not want the planetary orbits to interact gravitationally with one another, but you still want their spins to evolve. Isn't it analogous to using test particles in N-body simulations and wanting them to feel accelerations even if the force acting on them is zero?

I agree with Hanno that the right way to treat this would be to set the specific moment of inertia, so that it works whether or not the user then sets the mass to zero. We could add it as an extra parameter option, but I think it's much better if we just replace 'I' with 'Is' or something like that. Even though the paper says I, it points to all the documentation, and as long as the docs and all the examples show how to do that, I think it's totally fine. What do you think @tigerchenlu98 ?

hannorein commented 1 year ago

How about we just assume I is the specific moment of inertia if m=0?

tigerchenlu98 commented 1 year ago

Yeah exactly @dtamayo, it would be useful to have the effect in principle, I see reason that it would be useful.

I think I agree with @hannorein's suggestion - I think the fully dimensional moment of inertia will be more intuitive in most cases. Perhaps a warning when setting I for a test particle, telling the user that we will assume it's the specific MoI? If we're all in agreement with that I can add it in.

dtamayo commented 1 year ago

I'm curious what @hannorein thinks but I would implement it without a warning, and just document it explicitly in the docs and examples. If you're doing a bunch of test particle integrations it's a pain to keep seeing them, and it's better than getting nans back anyway

hannorein commented 1 year ago

Yeah. I think I would not bother with a warning message. But yes to documenting it!

tigerchenlu98 commented 1 year ago

Fixed in the code and noted in the documentation! Let me know if you think it's not clear enough/the implementation could be improved upon.

dtamayo commented 1 year ago

Thanks Tiger! The docs look perfect, and I added a couple notes to the ipynbs too