search

hpparvi / PyTransit

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

Stellar oblateness computation #135

Closed delinea closed 3 years ago

delinea commented 3 years ago

The computation of the stellar oblateness f in pytransit/models/numba/gdmodel.py assumes a density computed from a spherical star: rho = M_star * (4/3 * pi * R_star**3) to simplify the original expression of f: f = w**2 * R_star**3 / (2 * G * M_star)

One can however include the stellar oblateness in the expression of the density: rho = M_star * (4/3 * pi * R_star**3 * (1-f)) and then solve the resulting quadratic equation in f that gives a single possible value: f = 0.5 * (1 - sqrt(1 - 3.*w*w / (2.*G*pi*rho))) that should be returned by the stellar_oblateness(w, rho) function.

Of course, for slow rotators, w << 1 => f ~ 3.*w*w/(8.*G*pi*rho), which is the current implementation.

hpparvi commented 3 years ago

Changed as suggested!

delinea commented 3 years ago

Erratum: I recomputed the stellar oblateness from scratch instead of the reference formula f = w**2 * R_star**3 / (2 * G * M_star) and I found: f = 1 / (1 + (2 * G * M_star) / (w**2 * R_star**3))

Injecting the density of an oblate star rho = M_star * (4/3 * pi * R_star**3 * (1-f)) into the correct expression of f leads to f = 3*w**2 / (8 * G * pi * rho), which was the expression implemented initially...

Apologies for the erroneous suggestion I made earlier.

hpparvi commented 3 years ago

Thanks for double-checking this, I've reverted the function back to the original.