Transit time computation errors - outlier in a single transit

[ericagol]

For some initial orbital elements, the transit times show outliers.

Here is a minimum working example (this error occurs in Julia v1.7.2):

using NbodyGradient
t0 = 0.0; h=0.15; nplanet=2
elements = [1.0 0.0 0.0 0.0 0.0 1.5707963267948966 3.141592653589793; 0.0001206292283762164 3.3451735966291682 1.033198235742513 0.0894955105620105 -0.3857638886160483 1.5707963267948966 3.141592653589793; 2.6676800298570147e-5 6.742505903045666 5.086090365093386 0.09469923166430225 -0.38473827738282235 1.5707963267948966 3.141592653589793]
ic = ElementsIC(t0, nplanet+1, elements);
intr = Integrator(h, 0.0, 1100.0);
s = State(ic);
tt = TransitTiming(intr.tmax,ic);
intr(s,tt,grad=false)
println(tt.tt[2,241:244] .- tt.tt[2,240:243])  # Transit for the inner planet is late by ~14 days.
println(tt.tt[2,240:244]

If I rerun with a slightly different integration time, the error disappears:

intr = Integrator(h*0.99999, 0.0, 1100.0);
s = State(ic);
tt = TransitTiming(intr.tmax,ic);
intr(s,tt,grad=false)
println(tt.tt[2,241:244] .- tt.tt[2,240:243])
println(tt.tt[2,240:244])

But, decreasing changing h by a smaller factor, the outlier remains:

intr = Integrator(h*0.999999, 0.0, 1100.0);
s = State(ic);
tt = TransitTiming(intr.tmax,ic);
intr(s,tt,grad=false)
println(tt.tt[2,241:244] .- tt.tt[2,240:243])
println(tt.tt[2,240:244])

[ericagol]

It appears that the initial guess is within the time step, but then Newton's method causes the times to go out of the timestep. The cause appears to be that the g function becomes larger than the initial time step boundaries, rather than getting smaller. So, perhaps a hybrid finder, such as a bracketing followed by Newton's, or Brent's method, would fix this.

[ericagol]

A couple of options I am thinking of:

• [ ] If Newton's method takes us outside of the bounds of the time step, then carry out a bracketing and repeat until Newton's method is within bracketing bounds.
• [ ] If Newton's method takes us outside of the bounds of the time step, carry out a quadratic fit and find the zero within the time step.

[ericagol]

@langfzac I found the problem. This initial guess should be found by interpolation from the current time (where the "prior" state is saved), and then integrated backwards. But, the current formula is assuming that the integration is forwards in time:

https://github.com/ericagol/NbodyGradient.jl/blob/ca63dedbd4b4588c021fc8ba860b9fee1bbbedfa/src/transits/timing.jl#L20

Working through the interpolation, this formula should look like:

                dt0 = -gi*intr.h/(gi-tt.gsave[i]) # Starting estimate

i.e. tt.gsave[i] in the numerator should be gi instead. When I fix this, it fixes the outlier issue, as well as the other problems we were finding before!

[ericagol]

So, all to say, I think Newton's method works fine, and it was only because the starting guess was wrong that caused a problem with convergence and the occasional outlier. The only thing I'm puzzled about is why we didn't catch this before when we were checking the code... I thought we put in a temporary check that the interpolated initial guess was negative, and was greater than one time step before (i.e. between -h and 0). If we fix this, the code should run a bit faster since it will start much closer to the transit-time!