Instability for eccentric orbit

[odemangeon]

I noticed a non reproducibility of the code for eccentric orbit. I made a small script to show what is the problem:

import pytransit as pt
import numpy as np
import matplotlib.pyplot as pl
from math import acos, pi

m = pt.MandelAgol()

t = np.linspace(-0.3, 0.3, num=500)

k = 0.093390
u = [0.489582, 0.127812]
t0 = 0
p = 3.337757
a = 10.779853
i = 1.52
e = 0.13
# e = 0.0
w = pi * 0.5

l_res = []
for ii in range(10):
    pl.figure()
    l_res.append(m.evaluate(t, k, u, t0, p, a, i, e, w))
    pl.plot(t, l_res[ii])

l_res = np.asarray(l_res)

print(np.abs((l_res[:-1, :] - l_res[1:, :])).max())

The for loop execute the same computation 10 times and then I compute and print the maximum absolute difference between two consecutive executions. If the code is stable/reproducible it should always give me 0.0

When I run this script with 'e = 0.0', I get 0.0, but when I run it with e = 0.13, I get 0.00480200621537.

I made several test varying the other parameters and it looks that eccentricity is the parameter responsible for this instability. I also tried different values of e and it looks like the higher e is the higher is the maximum difference.

In the script there is a plot of the light-curve a every execution and it looks likes the resulting light-curve is made of two different transit shapes (which seems stable), but each time sample can be taken randomly from one of those two transit shapes.

Do you know where this comes from ?

PS: I tried reinstalling pytransit, but it didn't solve the problem.

[hpparvi]

Thanks for reporting this! I haven't seen this behaviour before, and will have a look as soon as I can.

[hpparvi]

The default solver for eccentric anomaly is misbehaving on some platforms. I haven't really found the reason behind this yet, since the Fortran code is rather straightforward. I've now changed the default solver to another while I'll try to find out what's causing the issue.