Problems with distrot near coplanar retrograde obliquity about vplanet HOT 9 CLOSED

Thanks for identifying this issue, Billy, and for the detailed error report. We'll look into it asap.

from vplanet.

Thanks for this additional example. It looks to me like it's the same problem: when i reaches 179 degrees, the NaNs appear. To be clear, this issue isn't just about initial conditions, should the natural behavior (or buggy behavior!) drive i to 179, then the NaNs will appear. It's not clear to me yet if the large eccentricity growth is realistic, but I suspect not, too. I should have been more careful in my previous comment. Perhaps you should keep i_0 < 30 degrees for now. The model is only approximate at large e and i, so you are probably pushing it beyond its limits. We need to implement some more checks and halts to make sure the code is stable for this type of application. I'm sorry it's given you trouble.

from vplanet.

SpiNBody+DistRot+POISE coupling has not been validated yet, so I wouldn't recommend it at this time. You can, however, run a different N-body model, and then use DistRot's capability to upload those data and compute the rotational axis. An example of this functionality is here: https://github.com/VirtualPlanetaryLaboratory/vplanet/tree/master/examples/SSDistOrbDistRot/marshnb. You can then couple in POISE as usual to compute the climate.

from vplanet.

I've taken a closer look here and there is indeed an issue, but I think the origin is in large inclinations. It seems if i > 179, then the denominator of a derivative can go to 0, and trigger a cascade of NaNs. As the derivative of i factors into the derivative of obliquity, it can look like the obliquity is the source, but it's inclination. Until we can identify a proper fix, extremely high inclinations should be avoided.

from vplanet.

Thanks for looking into this. I had performed simulations at other inclinations and ran into what seemed like a similar problem with the NaNs for other retrograde inclinations. I assume it is part of the same problem, but here's an example just in case it is something slightly different.
Inc[150].zip

from vplanet.

Since the problem is with the orbital integration using distorb, could I use the spinbody module instead? I expect that it would only be slower integration, right?

from vplanet.

PR #74 adds better handling of large mutual inclinations. A new option, dHaltMaxMutualInc, allows the user to set a mutual inclination to halt the execution. If it is set to 0 (the default), then this halt is not checked, which avoid the costly trig function calculations. Additionally, a warning is now issued if the mutual inclination is above the recommended limit for the Laplace-Lagrange (LL2) model (10 deg), or the 4th order (RD4) model (35 degrees). While it's difficult to identify exactly what conditions break the models, this patch at least warns the user of a likely breakdown, and allows users to set their own limits. It's been tested for DistOrb, but not for SpiNBody. Please give it a try and let us know if these additions help and/or are working as you expect.

from vplanet.

I used a separate integrator (Mercury) for the orbital evolution and was able to load those results into VPlanet as shown in the marshnb example with HNBody. That appears to work, where Distorb returns the expected obliquity variation when comparing to other software. It still returns 'nan' if I start with inc = 180 deg. (exact retrograde coplanar), so there might still be a problem for that special case.

from vplanet.

It's great to hear you find the NBody + DistRot functionality working. For the case of exactly coplanar orbits, we decided to leave the CassiniOne and CassiniTwo parameters as nans since they are undefined. The code should report that when executed, and it is also described in the long help for those two parameters.

from vplanet.