Here we create an example multicolour transit light curve dataset for the parameter estimation tutorials. The data comes in three versions a) True light curve, b) light curve with white noise, and c) light curve with white noise and systematic noise correlated with a simultaneously observed auxiliary parameter.
In [1]:
import pyfits as pf
import seaborn as sb
from os.path import join
from scipy.ndimage import gaussian_filter1d as GF
from pytransit import MandelAgol
from pytransit.orbits_f import orbits as of
from exotk.utils.orbits import as_from_rhop
warnings.simplefilter("ignore")
random.seed(0)
tc, p, rho, b, k = 1, 2.5, 1.5, 0.5, 0.1
a = as_from_rhop(rho,p)
i = arccos(b/a)
filters = 'u g r i z J2 H2 Ks'.split()
npb = len(filters)
ldcs = [[0.80,0.02],[0.61,0.16],[0.45,0.20],[0.36,0.20],
[0.30,0.20],[0.21,0.22],[0.12,0.25],[0.12,0.19]]
white_noise_std = 0.75e-3
obs_duration = of.duration_circular(p,a,i) + 3./24.
obs_cadence = 1./60./24.
obs_n = int(obs_duration/obs_cadence)
time = linspace(1-0.5*obs_duration,1+0.5*obs_duration, obs_n)
The white noise component is different for each observed passband.
In [2]:
noise_white = [normal(0,white_noise_std, size=obs_n) for ldc in ldcs]
We may still have a residual airmass-dependent baseline effect after dividing by the comparison stars. Here the difference with the next step (systematic noise) is that we know the relationship between the airmass and the effect, and can model it.
$$ F = F_0 \exp(-k\tau) $$
In [3]:
l = 0.05
airmass = 1.2+4*(time-1.05)**2
baseline = exp(-l*airmass)
#plot(time, flux_true[-1]*baseline);
The noise in photometric time series is never truly white, but also contains correlated (red) noise components arising from astrophysical and instrumental sources.
Here we create a systematic signal correlated with a simultaneously observed auxiliary parameter aux_true that varies during the observations more or less randomly. The observations of the auxiliary parameter, aux_obs, also contain a bit of (white) noise. In real life the auxiliary parameter could be the airmass, seeing, telescope rotator angle, air pressure, etc.
The systematic signal is assumed to be wavelength independent to keep things simple. However, the correlation between aux and the signal is nonlinear just to make the analysis a bit more interesting.
In [4]:
aux_true = GF(normal(size=obs_n), 5) + normal(0, 0.05, size=obs_n)
aux_obs = aux_true + normal(0, 0.025, size=obs_n)
noise_sys = aux_true**2 + sin(4*aux_true) + 0.6*sin(1.3+6.6*aux_true)
noise_sys *= 0.0015/noise_sys.std()
In [5]:
with sb.axes_style('white'):
fig, axs = subplots(1,3, figsize=(14,4))
axs[0].plot(time, aux_obs)
axs[1].plot(time, noise_sys)
axs[2].plot(aux_obs, noise_sys, '.')
fig.tight_layout()
In [6]:
tm = MandelAgol()
flux_true = [tm.evaluate(time, k, ldc, tc, p, a, i) for ldc in ldcs]
flux_wn = [flux+noise for flux,noise in zip(flux_true, noise_white)]
flux_sn = [flux+noise_sys for flux in flux_wn]
In [7]:
cols = ([pf.Column('time', 'D', 'BJD', array=time),
pf.Column('aux', 'D', '-', array=aux_obs)] +
[pf.Column('f_tr_{}'.format(filters[j]), 'D', array=flux_true[j]) for j in range(npb)] +
[pf.Column('f_wn_{}'.format(filters[j]), 'D', array=flux_wn[j]) for j in range(npb)] +
[pf.Column('f_sn_{}'.format(filters[j]), 'D', array=flux_sn[j]) for j in range(npb)])
hdul = pf.HDUList([pf.PrimaryHDU(), pf.BinTableHDU.from_columns(cols)])
hdul[1].header['extname'] = 'phot'
hdul.writeto(join('data','obs_data.fits'), clobber=True)
In [8]:
with sb.axes_style('white'):
th, y_offset = time*24-24, 0.005
fig,ax = subplots(1,3,figsize=(12,8))
cs = [sb.desaturate(cm.spectral((j+1)/float(npb+1)), 0.75) for j in range(npb)]
[ax[0].plot(th, flux-j*y_offset, '-', c=cs[j]) for j,flux in enumerate(flux_true)]
[ax[1].plot(th, flux-j*y_offset, '-', c=cs[j]) for j,flux in enumerate(flux_wn)]
[ax[2].plot(th, flux-j*y_offset, '-', c=cs[j]) for j,flux in enumerate(flux_sn)]
setp(ax, xlim=th[[0,-1]], ylim=(0.951,1.004), yticks=[])
fig.tight_layout()