hpparvi / PyTransit

Fast and easy exoplanet transit light curve modelling.
GNU General Public License v2.0
99 stars 23 forks source link

Issue with oblate star model #86

Closed shashankdholakia closed 3 years ago

shashankdholakia commented 3 years ago

Getting different lightcurves with the latest version of PyTransit (compared to version 2.4.0). Lightcurves in general have the same duration but the extent of the gravity darkening is far less. Plots (line in black) showing same parameter LC before and after update. lc_full_omega08 pytransit_update

hpparvi commented 3 years ago

Thanks for pointing this out Shashank! Can you send me the code that you used to produce the plots and I'll have a look at what's happening.

shashankdholakia commented 3 years ago

The code I used is not too different from the simple working example in the oblate star notebook:

import pytransit
from pytransit import OblateStarModel, QuadraticModel
print(pytransit.__file__)
tmo = OblateStarModel(sres=1000, pres=100, tres=100, rstar = 1.561)
tmc = QuadraticModel(interpolate=False)

inc = 0
i_p = 88.01
long_p = 92
r_p = 0.1087
G_mks = 6.67e-11
Msun = 1.989e+30
Rsun = 6.95700e8
M_star = 1.59
R_star = 1.561

time = np.linspace(1.2198696-0.15, 1.2198696+0.15, 500)
tmo.set_data(time)
tmc.set_data(time)

omega_s = np.abs(0.8)
k = np.array([])
p_rot = (2*np.pi*(R_star*Rsun)**(3/2))/(omega_s*(G_mks*M_star*Msun)**(1/2)*(60*60*24))

t0, p, a, i, az, e, w = 0.0, 1.2198681, 3.614, i_p*(np.pi/180), long_p*(np.pi/180.), 0.0, 0.0 
rho, rperiod, tpole, phi, beta = 0.778, p_rot,7600, (inc)*(np.pi/180.), 0.22

flux_om = tmo.evaluate_ps(k, rho, rperiod, tpole, phi, beta, ldc, t0, p, a, i, az, e, w)
shashankdholakia commented 3 years ago

After a look-through of the source and comparison to Barnes et al. 2009 (which I am also using), I found a potential bug. On line 223 of pytransit/models/numba/osmodel: y = y * rstar * (1 - feff)

I think when plugging the values x and y into the Von Zeipel Law there shouldn't be this extra factor of (1-f_eff). When I remove the (1-f_eff), the lightcurves match PyTransit 2.4.0 and also my generated lightcurves.

Let me know if you think this is indeed a bug and I can submit a pull request.

hpparvi commented 3 years ago

You're right, this is a bug I introduced when I changed the way the model is evaluated.

Thanks again for noticing (and fixing) this! I've merged the fix and will work on the model more this and next week.