# PyTransit: fast and easy exoplanet transit modelling in Python.
# Copyright (C) 2010-2020 Hannu Parviainen
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <https://www.gnu.org/licenses/>.
from pathlib import Path
from typing import Union
import pandas as pd
import astropy.units as u
from numba import njit
from numpy import ndarray, pi, atleast_1d, zeros, exp, diff, log, sqrt, nan, vstack, tile, linspace, cos, sin, sum, \
isfinite
from scipy.constants import c,h,k,G
from scipy.optimize import brentq
from ..stars import create_husser2013_interpolator, create_bt_settl_interpolator
#from ..contamination.filter import Filter
NPType = Union[float,ndarray]
mj2kg = u.M_jup.to(u.kg)
ms2kg = u.M_sun.to(u.kg)
d2s = u.day.to(u.s)
[docs]def equilibrium_temperature(tstar: NPType, a: NPType, f: NPType, ab: NPType) -> NPType:
"""Planetary equilibrium temperature [K].
Parameters
----------
tstar
Effective stellar temperature [K]
a
Scaled semi-major axis [Rsun]
f
Redistribution factor
ab
Bond albedo
Returns
-------
Teq : float or ndarray
Equilibrium temperature [K]
"""
return tstar * sqrt(1 / a) * (f * (1 - ab)) ** 0.25
# Thermal emission
# ================
[docs]def planck(t: NPType, l: NPType) -> NPType:
"""Radiance of a black body as a function of wavelength.
Parameters
----------
t
Black body temperature [K]
l
Wavelength [m]
Returns
-------
L : float or ndarray
Back body radiance [W m^-2 sr^-1]
"""
return 2 * h * c ** 2 / l ** 5 / (exp(h * c / (l * k * t)) - 1)
@njit
def summed_planck(teff, wl, tm):
teff = atleast_1d(teff)
flux = zeros(teff.shape[0])
for i in range(flux.size):
flux[i] = sum(tm*(2*h*c**2 / wl**5 / (exp(h*c / (wl*k*teff[i])) - 1.)))
return flux
[docs]def emission(tp: NPType, tstar: NPType, k: NPType, flt) -> NPType:
"""Thermal emission from the planet.
Parameters
----------
tp : float or ndarray
Equilibrium temperature of the planet [K]
tstar : float or ndarray
Effective temperature of the star [K]
k : float or ndarray
Planet-star radius ratio
flt: Filter
Passband transmission
Returns
-------
float or ndarray
References
"""
wl = linspace(flt.wl_min, flt.wl_max, 100)
tm = flt(wl)
return k**2 * (summed_planck(tp, wl, tm) / summed_planck(tstar, wl, tm))
# Doppler boosting
# ================
[docs]def doppler_beaming_factor(teff: float, flt, dataset: str = 'BT-Settl'):
"""The photon weighted bandpass-integrated Doppler beaming factor.
Parameters
----------
teff
Effective temperature of the star [K]
flt
Passband transmission
dataset
Either 'BT-Settl' or 'Husser2013'
Returns
-------
float
The photon weighted bandpass-integrated Doppler beaming factor.
Notes
------
See Bloemen et al., MNRAS, 410 (2010) and Claret, A. et al., A&A, 641 (2020).
"""
if dataset == 'Husser2013':
ip = create_husser2013_interpolator()
else:
ip = create_bt_settl_interpolator()
wl_nm = ip.grid[1][1:-1]
wl = 1e-9 * wl_nm
fl = ip(vstack((tile(teff, wl_nm.size), wl_nm)).T)
b = 5. + diff(log(fl)) / diff(log(wl))
w1 = flt(wl_nm)[1:] * wl[1:] * fl[1:]
w2 = flt(wl_nm)[:-1] * wl[:-1] * fl[:-1]
w = 0.5 * (w1 + w2)
m = isfinite(b)
return sum(w[m] * b[m]) / sum(w[m])
[docs]def doppler_beaming_amplitude(mp: NPType, ms: NPType, period: NPType, b: NPType) -> NPType:
"""The amplitude of the Doppler beaming signal.
Calculates the amplitude of the Doppler beaming signal following the approach
described by Loeb & Gaudi in [Loeb2003]_ . Note that you need to pre-calculate the photon-weighted bandpass-integrated
doppler Beaming factor [Bloemen2010]_ [Barclay2012]_ for the star and the instrument using ``doppler_beaming_factor``.
Parameters
----------
mp : float or ndarray
Planetary mass [MJup]
ms : float or ndarray
Stellar mass [MSun]
period : float or ndarray
Orbital period [d]
b : float or ndarray
Photon-weighted bandpass-integrated Doppler beaming factor
Returns
-------
float or ndarray
Doppler boosting signal amplitude
References
----------
.. [Loeb2003] Loeb, A. & Gaudi, B. S. Periodic Flux Variability of Stars due to the Reflex Doppler Effect Induced by
Planetary Companions. Astrophys. J. 588, L117–L120 (2003).
.. [Bloemen2010] Bloemen, S. et al. Kepler observations of the beaming binary KPD 1946+4340. MNRAS 410, (2010).
.. [Barclay2012] Barclay, T. et al. PHOTOMETRICALLY DERIVED MASSES AND RADII OF THE PLANET AND STAR IN THE TrES-2
SYSTEM. AspJ 761, 53 (2012).
"""
return b / c *(2*pi*G/(d2s*period))**(1/3) * ((mp*mj2kg)/(ms*ms2kg)**(2/3))
# Ellipsoidal variations
# ======================
def ellipsoidal_variation(f: NPType, theta: NPType, mp: float, ms: float, a: float, i: float, e: float, u: float, g: float):
raise NotImplementedError
#return ellipsoidal_variation_amplitude(mp, ms, a, i, u, g)
[docs]def ellipsoidal_variation_signal(f: NPType, theta: NPType, e: float) -> NPType:
"""
Parameters
----------
f
True anomaly [rad]
theta
Angle between the line-of-sight and the star-planet direction
e
Eccentricity
Returns
-------
"""
return -((1 + e * cos(f)) / (1-e**2))**3 * cos(2*theta)
[docs]def ellipsoidal_variation_amplitude(mp: NPType, ms: NPType, a: NPType, i: NPType, u: NPType, g: NPType) -> NPType:
"""The amplitude of the ellipsoidal variation signal.
Calculates the amplitude of the ellipsoidal variation signal following the approach described by
Lillo-Box et al. in [Lillo-Box2014]_, page 11.
Parameters
----------
mp : float or ndarray
Planetary mass [MJup]
ms : float or ndarray
Stellar mass [MSun]
a : float or ndarray
Semi-major axis of the orbit divided by the stellar radius
i : float or ndarray
Orbital inclination [rad]
u : float or ndarray
Linear limb darkening coefficient
g : float or ndarray
Gravity darkening coefficient
Returns
-------
ev_amplitude: float or ndarray
The amplitude of the ellipsoidal variation signal
References
----------
.. [Lillo-Box2014] Lillo-Box, J. et al. Kepler-91b: a planet at the end of its life. A&A 562, A109 (2014).
"""
ae = 0.15 * (15 + u) * (1 + g) / (3 - g)
return ae * (mp*mj2kg)/(ms*ms2kg) * a**-3 * sin(i)**2
[docs]def reflected_fr(a: NPType, ab: NPType, r: NPType = 1.5) -> NPType:
"""Reflected flux ratio per projected area element.
Parameters
----------
a
Scaled semi-major axis [Rsun]
ab
Bond albedo
r
Inverse of the phase integral
Returns
-------
fr : float
Reflected flux ratio
"""
return r * ab / a ** 2
[docs]def flux_ratio(tstar: NPType, a: NPType, f: NPType, ab: NPType, l: NPType, r: NPType = 1.5, ti: NPType = 0) -> NPType:
"""Total flux ratio per projected area element.
Parameters
----------
tstar
Effective stellar temperature [K]
a
Scaled semi-major axis [Rs]
f
Redistribution factor
ab
Bond albedo
l
Wavelength [m]
r
Inverse of the phase integral
ti
Temperature [K]
Returns
-------
fr: float
Total flux ratio
"""
return reflected_fr(a, ab, r) + thermal_fr(tstar, a, f, ab, l, ti)
[docs]def solve_teq(fr, tstar, a, ab, l, r=1.5, ti=0):
"""Solve the equilibrium temperature.
Parameters
----------
fr
Flux ratio
tstar
Effective stellar temperature [K]
a
Scaled semi-major axis [Rs]
ab
Bond albedo
l
Wavelength [m]
r
Inverse of the phase integral
ti
Temperature [K]
Returns
-------
Teq : float or ndarray
Equilibrium temperature
"""
Bs = planck(tstar, l)
try:
return brentq(lambda Teq: reflected_fr(a, ab, r) + planck(Teq + ti, l) / Bs - fr, 5, tstar)
except ValueError:
return nan
[docs]def solve_ab(fr, tstar, a, f, l, r=1.5, ti=0):
"""Solve the Bond albedo.
Parameters
----------
fr :
Flux ratio [-]
tstar :
Effective stellar temperature [K]
a :
Scaled semi-major axis [Rs]
A :
Bond albedo [-]
l :
Wavelength [m]
r :
Inverse of the phase integral [-]
ti :
Temperature [K]
Returns
-------
A : Bond albedo
"""
try:
return brentq(lambda ab: reflected_fr(a, ab, r) + thermal_fr(tstar, a, f, ab, l, ti) - fr, 0, 0.3)
except ValueError:
return nan
[docs]def solve_redistribution(fr, tstar, a, ab, l):
"""Solve the redistribution factor.
Parameters
----------
fr :
Flux ratio [-]
tstar :
Effective stellar temperature [K]
a :
Scaled semi-major axis [Rs]
ab :
Bond albedo [-]
l :
Wavelength [m]
r :
Inverse of the phase integral [-]
Ti : t
Temperature [K]
Returns
-------
f : Redistribution factor
"""
Teqs = solve_teq(fr, tstar, l)
return brentq(lambda f: equilibrium_temperature(tstar, a, f, ab) - Teqs, 0.25, 15)