Según Seager & Mallén-Ornelas (2003),
La geometría del tránsito puede ser descrita por (se presentan imágenes del trabajo de Seager & Mallén-Ornelas):
A partir de las características de la curva de luz de un tránsito, teniendo en cuenta la geometría del evento y la Tercera Ley de Kepler, se pueden obtener cuatro parámetros derivables de observables del sistema:
Si se considera la relación masa-radio para una estrella $$R_* = kM_*^x $$ para algunas constantes $k,x$ se obtiene un sistema de cinco ecuaciones y cinco incógnitas, y se pueden derivar las cantidades físicas una por una.
En este paso se reliza una estandarización de los cuartos de modo que queden en torno al promedio de las medianas. Luego, se realizan dos pasos de median filter con ventanas de 11 valores y se normaliza el flujo. Para el período sugerido (704.2 días), se grafica el flujo en función de la fase. En verde se incluyen las barras de error asociadas al flujo.
In [1]:
import numpy as np
from scipy.signal import medfilt
import matplotlib.pyplot as plt
import kplr
%matplotlib inline
client = kplr.API()
koi = client.koi(1274.01)
lcs = koi.get_light_curves(short_cadence=True)
p = 704.2
time, flux, ferr, med = [], [], [], []
for lc in lcs:
with lc.open() as f:
# The lightcurve data are in the first FITS HDU.
hdu_data = f[1].data
time.append(hdu_data["time"][~np.isnan(hdu_data["pdcsap_flux"])])
flux.append(hdu_data["pdcsap_flux"][~np.isnan(hdu_data["pdcsap_flux"])])
ferr.append(hdu_data["pdcsap_flux_err"][~np.isnan(hdu_data["pdcsap_flux"])])
# Ignora los NaN al hacer append
normFlux, normFerr, phase = flux, ferr, time
for i in range(0,len(flux)):
med.append(np.median(flux[i]))
prom = np.mean(med)
for i in range(0,len(flux)):
normFlux[i] = normFlux[i] - (med[i] - prom)
normFlux[i] = medfilt(normFlux[i], 11)
fig, ax = plt.subplots(2,1,figsize=(15,20))
for i in range(0,len(ax)):
ax[i].set_ylim(0.996,1.0007)
ax[i].set_title('KOI-1274.01')
ax[i].set_xlabel('Phase',size=16)
ax[i].set_ylabel('Normalized flux',size=14)
for i in range(0,len(normFlux)):
normFlux[i] = normFlux[i]/prom
normFerr[i] = normFerr[i]/prom
phase[i] = time[i]/p %1
ax[0].errorbar(phase[i], normFlux[i], normFerr[i], fmt='g.', ecolor='green', ms = 3)
ax[0].plot(phase[i], normFlux[i],'k.')
ax[1].errorbar(phase[i], normFlux[i], normFerr[i], fmt='g.', ecolor='green', ms = 3)
ax[1].plot(phase[i], normFlux[i],'k--', alpha=.2)
ax[1].set_xlim(0.699,0.7005)
plt.show()
plt.close()
A simple vista, aproximadamente, se tiene
$T_T = 0.7$ d
$T_F = 0.4$ d
$\Delta F = 0.003$
Con esto, se calculan los cuatro parámetros a partir de observables (sin considerar en este punto una ley de potencias que relacione masa y radio).
Como las ecuaciones son las del trabajo de Seager & Mallén-Ornelas, las suposiciones necesarias para determinar las cantidades son las que allí se hacen.
In [2]:
df = 0.003
tt = 0.7
tf = 0.4
sintf = np.sin(tf*np.pi/p)**2 # un par de variables auxiliares
sintt = np.sin(tt*np.pi/p)**2
ratio = np.sqrt(df) #Rp/R*
b = np.sqrt( ((1-ratio)**2 - (sintf)/(sintt) *(1+ratio)**2) /(1-(sintf/sintt)) )
aR = np.sqrt( ((1+ratio)**2 - b**2 *(1-sintt)) /sintt )
i = np.arccos(b/aR)
i = np.degrees(i)
rho = aR**3 * 365.25**2 / 215**3 / p**2
print 'Rp/R* \t = \t' + repr(ratio)
print 'b \t = \t' + repr(b)
print 'a/R* \t = \t' + repr(aR)
print 'i \t = \t' + repr(i)
print 'rho \t = \t' + repr(rho) + ' densidades solares'
Se calculan los coeficientes según se detalla en en Espinoza & Jordán 2015, a partir de una versión modificada del código disponible en https://github.com/nespinoza/limb-darkening.
In [3]:
from scipy.optimize import leastsq
from scipy.interpolate import UnivariateSpline
import scipy.integrate as integrate
w, r = np.loadtxt('kepler_response_hires1.txt', unpack=True)
w = 10*w
S = UnivariateSpline(w,r,s=0,k=1)
min_w = min(w)
max_w = max(w)
idx = np.where((w>min_w)&(w<max_w))[0]
S_wav = np.append(np.append(min_w,w[idx]),max_w)
S_res = np.append(np.append(S(min_w),r[idx]),S(max_w))
I = np.array([])
wavelengths = np.array([])
f = open('grav_4.5_lh_1.25.dat','r')
counter = 0
while(True):
l = f.readline()
if(l==''):
break
# If no jump of line or comment, save the intensities:
if(l[0]!='#' and l[:3]!='\n'):
splitted = l.split('\t')
if(len(splitted)==18):
splitted[-1] = (splitted[-1])[:-1] # The last one always has a jump of line (\n), so erase it.
wavelength = np.double(splitted[0])*10 # Convert wavelengths, which are in nanometers, to angstroms.
intensities = np.double(np.array(splitted[1:])) # Get the intensities.
ndigits = len(str(int(intensities[1])))
# Only if I(1) is different from zero, fit the LDs:
if(intensities[0]!=0.0):
intensities[1:] = intensities[1:]/1e5 # Kurucz doesn't put points on his files (e.g.: 0.8013 is 8013).
intensities[1:] = intensities[1:]*intensities[0] # All the rest of the intensities are normalized w/r to the center one.
if(counter == 0):
I = intensities
else:
I = np.vstack((I,intensities))
wavelengths = np.append(wavelengths,wavelength)
counter = counter + 1
f.close()
mu = np.array([1.0,0.9,0.8,0.7,0.6,0.5,0.4,0.3,0.25,0.2,0.15,0.125,0.1,0.075,0.05,0.025,0.01])
# Define the number of mu angles at which we will perform the integrations:
nmus = len(mu)
# Now integrate intensity through each angle:
I_l = np.array([])
for i in range(nmus):
# Interpolate the intensities:
Ifunc = UnivariateSpline(wavelengths,I[:,i],s=0,k=1)
integrand = S_res*Ifunc(S_wav)
integration_results = np.trapz(integrand, x=S_wav)
I_l = np.append(I_l,integration_results)
I0 = I_l/(I_l[0]) # Normalize profile with respect to I(mu = 1):
# Define A matrix for the linear system:
A = np.zeros([2,2])
# Define b vector for the linear system:
b = np.zeros(2)
# Obtain the alpha_n_k and beta_k that fill the A matrix and b vector:
for n in range(1,3,1):
for k in range(1,3,1):
A[n-1,k-1] = sum(((1.0-mu)**n)*((1.0-mu)**k))
b[n-1] = sum(((1.0-mu)**n)*(1.0-I0))
u = list(np.linalg.solve(A,b))
print u
Se ejecuta batman como se explica en la documentación, entregando como parámetros los valores obtenidos a lo largo de este trabajo.
In [4]:
import batman
params = batman.TransitParams() #object to store transit parameters
params.t0 = 0. #time of inferior conjunction
params.per = p #orbital period
params.rp = ratio #planet radius (in units of stellar radii)
params.a = aR #semi-major axis (in units of stellar radii)
params.inc = i #orbital inclination (in degrees)
params.ecc = 0. #eccentricity
params.w = 90. #longitude of periastron (in degrees)
params.limb_dark = "quadratic" #limb darkening model
params.u = u #limb darkening coefficients
t = np.linspace(-0.025, 0.025, 100) #times at which to calculate light curve
m = batman.TransitModel(params, t) #initializes model
fluxBatman = m.light_curve(params) #calculates light curve
plt.plot(t, fluxBatman)
plt.xlabel("Time from central transit")
plt.ylabel("Relative flux")
plt.show()
##############
#oc
Se repiten los pasos para el planeta anterior. Dado que los valores del flujo son más dispersos, se realiza el paso de median filter con ventanas de 25 valores y se centra el gráfico en un posible tránsito, que corresponde al sector alrededor de la zona de mínimo flujo (obviando un punto, seguramente espuria, en phase ~ 0.35). Si el planeta fuera del radio de la Tierra y su estrella fuera del radio del Sol, entonces $$R_p/R_* = 0.009 \rightarrow \Delta F = 0,000081$$
Esto se acerca a lo que muestra el gráfico.
In [5]:
koi = client.koi(7016.01)
lcs = koi.get_light_curves(short_cadence=True)
p = koi.koi_period
time, flux, ferr, med = [], [], [], []
for lc in lcs:
with lc.open() as f:
# The lightcurve data are in the first FITS HDU.
hdu_data = f[1].data
time.append(hdu_data["time"][~np.isnan(hdu_data["pdcsap_flux"])])
flux.append(hdu_data["pdcsap_flux"][~np.isnan(hdu_data["pdcsap_flux"])])
ferr.append(hdu_data["pdcsap_flux_err"][~np.isnan(hdu_data["pdcsap_flux"])])
# Ignora los NaN al hacer append
normFlux, normFerr, phase = flux, ferr, time
for i in range(0,len(flux)):
med.append(np.median(flux[i]))
prom = np.mean(med)
for i in range(0,len(flux)):
normFlux[i] = normFlux[i] - (med[i] - prom)
normFlux[i] = medfilt(normFlux[i], 25)
fig, ax = plt.subplots(2,1,figsize=(15,20))
for i in range(0,len(ax)):
ax[i].set_ylim(0.996,1.0007)
ax[i].set_title('KOI-7016.01')
ax[i].set_xlabel('Phase',size=16)
ax[i].set_ylabel('Normalized flux',size=14)
for i in range(0,len(normFlux)):
normFlux[i] = normFlux[i]/prom
normFerr[i] = normFerr[i]/prom
phase[i] = time[i]/p %1
ax[0].errorbar(phase[i], normFlux[i], normFerr[i], fmt='g.', ecolor='green', ms = 3)
ax[0].plot(phase[i], normFlux[i],'k.')
ax[1].errorbar(phase[i], normFlux[i], normFerr[i], fmt='g.', ecolor='green', ms = 3)
ax[1].plot(phase[i], normFlux[i],'k.', alpha=.2)
ax[1].set_xlim(0.762,0.782)
ax[1].set_ylim(0.9985,1.001)
plt.show()
plt.close()