Source code for pytransit.contamination.contamination

#  PyTransit: fast and easy exoplanet transit modelling in Python.
#  Copyright (C) 2010-2019  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/>.

"""Module to model flux contamination in transit light curves.

Module to model flux contamination in transit light curves. The `pytransit.contamination` module contains two flux
contamination models: a physical model based on simulated stellar spectra (`SMContamination`), and a simple one
based on black body approximation (`BBContamination`). The models share an common interface.
"""

from os.path import join
from typing import Union, Iterable

import pandas as pd
import xarray as xa

from numba import njit
from matplotlib.pyplot import subplots, setp
from numpy import transpose, newaxis, uint32, ndarray, asarray, zeros_like, log, exp, ceil, linspace, array, nan
from pandas import DataFrame
from scipy.interpolate import interp1d, RegularGridInterpolator

from .instrument import Instrument
from ..stars import read_bt_settl_table, read_husser2013_table
from ..utils.phasecurves import planck


[docs]@njit def contaminate_light_curve(flux: ndarray, contamination: ndarray, pbids: ndarray) -> ndarray: """Contaminates a transit light curve. Contaminates a transit light curve with npb passbands. Parameters ---------- flux: 1d array-like Transit light curve with npb passbands. contamination: 1d array-like Array of per-passband contamination values. pbids: 1d array-like Passband indices that map each light curve element to a single passband. Returns ------- Contaminated transit light curve """ npt = flux.size contaminated_flux = zeros_like(flux) for i in range(npt): contaminated_flux[i] = contamination[pbids[i]] + (1.0-contamination[pbids[i]])*flux[i] return contaminated_flux
class _BaseContamination: """A base class to model flux contamination (blending) in transit light curves. A base class to model flux contamination (blending) in transit light curves. Notes ----- This class does not implement all the necessary functionality, and is not meant to be used as-is. """ def __init__(self, instrument: Instrument, ref_pb: str) -> None: """ Parameters ---------- instrument Instrument configuration ref_pb name of the reference passband """ self.instrument = instrument self._ri = self.instrument.pb_names.index(ref_pb) self._rpb = ref_pb def relative_fluxes(self, teff: Union[float, ndarray]): raise NotImplementedError def relative_flux_mixture(self, teffs: Iterable, fractions: Iterable): raise NotImplementedError def contamination(self, cr: float, teff1: float, teff2: float): raise NotImplementedError def exposure_times(self, teff, rtime, rflux=1.0, tflux=1.0): """Exposure times that give equal flux as in the reference passband Parameters ---------- teff : float Effective stellar temperature [K] rtime : float Exposure time in the reference passband rflux : float, optional Flux in the reference passband with the reference exposure time tflux : float, optional Target flux in the reference passband """ return rtime * tflux / rflux / self.relative_fluxes(teff) def c_as_pandas(self, ci, teff1, teff2): """Contamination as a pandas DataFrame.""" return pd.DataFrame(self.contamination(ci, teff1, teff2), columns=pd.Index(self.instrument.pb_names, name='passband'), index=pd.Index(uint32(teff2), name='teff')) def c_as_xarray(self, ci, teff1, teff2): """Contamination as an xarray DataArray.""" from xarray import DataArray return DataArray(self.contamination(ci, teff1, teff2), name='contamination', dims=['teff', 'passband'], coords=[uint32(teff2), self.instrument.pb_names], attrs={'TEff_1': teff1}) def plot(self, teff1: float, teff2: float, cref, wlref=600, wlmin=305, wlmax=995, nwl=500, figsize=None, axs=None): """Plots the contamination model. Parameters ---------- teff1 teff2 cref wlref wlmin wlmax nwl figsize Returns ------- """ wl = linspace(wlmin, wlmax, nwl) ft = (1 - cref) * self.relative_flux(teff1, wl, wlref) fc = cref * self.relative_flux(teff2, wl, wlref) if axs is None: fig, axs = subplots(2, 1, figsize=figsize, sharex='all', constrained_layout=True) else: fig = None axs[0].plot(wl, ft, label='host') axs[0].plot(wl, fc, label='contaminant') axs[1].plot(wl, fc / (ft + fc), 'k') for ax in axs: ax.axhline(cref, ls='--', c='0.75') ax.axvline(wlref, ls='--', c='0.75') axs[0].legend(fontsize='small') setp(axs, xlim=(wlmin, wlmax)) setp(axs[0], ylabel='Relative flux') setp(axs[1], ylim=(-0.05, 1.05), xlabel='Wavelength [nm]', ylabel='Contamination') #fig.tight_layout() return fig
[docs]class BBContamination(_BaseContamination): """Third light contamination based on black body approximation. This class offers a simple black-body model for flux contamination in which the target star and the contaminant(s) are approximated as black bodies with effective temperatures Tt, Tc1, Tc2, ..., Tcn. """ def __init__(self, instrument: Instrument, ref_pb: str, delta_l: float = 10): """ Parameters ---------- instrument Instrument configuration. ref_pb Reference passband name. delta_l Wavelength grid spacing in nm. """ super().__init__(instrument, ref_pb) self._delta_l = delta_l self._wl_grids = [] for f in self.instrument.filters: nwl = int(ceil((f.wl_max - f.wl_min) / self._delta_l)) self._wl_grids.append(linspace(f.wl_min, f.wl_max, nwl))
[docs] @staticmethod def absolute_flux(teff: float, wl: Union[float, Iterable]) -> ndarray: """The absolute flux given an effective temperature and wavelength. Parameters ---------- teff The effective temperature in K wl The wavelength (or an array of) in nm Returns ------- Black body spectral radiance. """ return planck(teff, 1e-9 * wl)
[docs] @staticmethod def relative_flux(teff: float, wl: Union[float, ndarray], wlref: float) -> ndarray: """The black body flux normalized to a given reference wavelength. Parameters ---------- teff The effective temperature of the radiating body [K] wl The wavelength [nm] wlref The reference wavelength [nm] Returns ------- The black body flux normalized to a given reference wavelength """ return planck(teff, 1e-9 * wl) / planck(teff, 1e-9 * wlref)
[docs] def absolute_fluxes(self, teff: float) -> ndarray: """Calculates the integrated absolute fluxes for all filters for a star with the given effective temperature Parameters ---------- teff The effective temperature of the radiating body [K] Returns ------- The integrated absolute fluxes for the filters in the instrument. """ return array( [(self.absolute_flux(teff, g) * f(g)).mean() for f, g in zip(self.instrument.filters, self._wl_grids)])
[docs] def relative_fluxes(self, teff: Union[float, ndarray]) -> ndarray: """Calculates the integrated fluxes for all filters normalized to the reference passband. Parameters ---------- teff The effective temperature of the radiating body [K] Returns ------- The integrated fluxes for all filters normalized to the reference passband. """ fluxes = self.absolute_fluxes(teff) return fluxes / fluxes[self._ri]
[docs] def contamination(self, cref: float, teff1: float, teff2: float) -> ndarray: """Calculates the contamination factors for all the filters given the contamination in the reference passband. Parameters ---------- cref Reference passband contamination teff1 Host star effective temperature teff2 Contaminant effective temperature Returns ------- Contamination factors for all the filters. """ ft = (1.0 - cref) * self.relative_fluxes(teff1) fc = cref * self.relative_fluxes(teff2) return fc / (ft + fc)
[docs]class SMContamination(_BaseContamination): """A class that models flux contamination based on stellar spectrum models. """ def __init__(self, instrument: Instrument, ref_pb: str = None, data: str = 'BT-SETTL') -> None: """ Parameters ---------- instrument Instrument configuration ref_pb name of the reference passband """ super().__init__(instrument, ref_pb) self.data = data if data.lower() == 'bt-settl': self._spectra: DataFrame = read_bt_settl_table() elif data.lower() == 'husser2013': self._spectra: DataFrame = read_husser2013_table() else: raise ValueError('The dataset has to be either "BT-SETTL" or "Husser2013"') self.teff = teff = self._spectra.index.values self.wl = wl = self._spectra.columns.values # 2D interpolator for the spectrum table # -------------------------------------- self._rgi = RegularGridInterpolator((self.teff, self.wl), self._spectra.values) # Dataframe indices # ----------------- self.ipb = pd.Index(self.instrument.pb_names, name='passband') self.iteff = pd.Index(teff, name='teff') # Per-passband fluxes # ------------------- self._compute_relative_flux_tables(0, ref_pb) def _integrate_fluxes(self, rdc=None): i, sp, wl = self.instrument, self._spectra, self.wl int_fluxes = transpose([(sp.values * f(wl) * qe(wl)).sum(1) for f, qe in zip(i.filters, i.qes)]) self._af = pd.DataFrame(int_fluxes, columns=self.ipb, index=self.iteff) def _compute_relative_flux_tables(self, rdc=None, rpb=None): self._integrate_fluxes(rdc) af = self._af.values if rpb is not None: self._rpb = rpb self._ri = self.instrument.pb_names.index(rpb) self._rf = pd.DataFrame(af / af[:, self._ri][:, newaxis], columns=self.ipb, index=self.iteff) self._ip = interp1d(self.iteff.values, self._rf.values.T, bounds_error=False, fill_value=nan) self._ip_abs = interp1d(self.iteff.values, self._af.values.T, bounds_error=False, fill_value=nan)
[docs] def absolute_flux(self, teff: float, wl: Union[float, Iterable]) -> ndarray: """The absolute flux given an effective temperature and a set of wavelength. Parameters ---------- teff: float The effective temperature [K] wl: array-like The wavelengths to calculate the flux in [nm] Returns ------- Spectral radiance. """ return self._rgi((teff, wl))
[docs] def relative_flux(self, teff: float, wl: Union[float, ndarray], wlref: float) -> ndarray: """The stellar flux normalized to a given reference wavelength. Parameters ---------- teff: float The effective temperature of the radiating body [K] wl: array-like The wavelengths to calculate the flux in [nm] wlref: float The reference wavelength [nm] Returns ------- The flux normalized to a given reference wavelength """ return self.absolute_flux(teff, wl) / self.absolute_flux(teff, wlref)
[docs] def relative_fluxes(self, teff: Union[float, ndarray], rdc=None, rpb=None): if (rpb is not None and rpb != self._rpb) or (rdc is not None): self._compute_relative_flux_tables(rdc, rpb) return self._ip(teff)
[docs] def relative_flux_mixture(self, teffs, fractions, rdc=None): teffs = asarray(teffs) fractions = asarray(fractions) x = zeros_like(teffs) assert len(teffs) == len(fractions) + 1 assert 0.0 < fractions.sum() < 1.0 x[1:] = fractions x[0] = 1. - x[1:].sum() return (self.relative_fluxes(teffs, rdc) * x).sum(1)
[docs] def contamination(self, cref: Union[float, ndarray], teff1: Union[float, ndarray], teff2: Union[float, ndarray], absolute: bool = False): """Contamination given reference contamination, host TEff, and contaminant TEff(s) Per-passband contamination given the contamination in the reference passband and TEffs of the two stars. Parameters ---------- cref : float contamination in the reference passband teff1 Effective stellar temperature [K] teff2 Effective stellar temperature [K] Returns ------- Per-passband contamination """ if absolute: teff1, teff2 = asarray(teff1), asarray(teff2) a, b = self._ip_abs(teff1), self._ip_abs(teff2) else: cref, teff1, teff2 = asarray(cref), asarray(teff1), asarray(teff2) a = (1.0 - cref) * self._ip(teff1) b = cref * self._ip(teff2) return (b / (a + b)).T