<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Tarea 2, parte 2
<hr>
Pregunta 1
Según <a href="https
Step1:  <hr>
Pregunta 3
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.
Step2:  <hr>
Pregunta 4
Se calculan los coeficientes según se detalla en en <a href="https
Step3:  Se ejecuta batman como se explica en la documentación, entregando como parámetros los valores obtenidos a lo largo de este trabajo.
<hr>
Pregunta Bonus
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.
<ASSISTANT_TASK:>
Python Code:
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()
Explanation: Tarea 2, parte 2
<hr>
Pregunta 1
Según <a href="https://ui.adsabs.harvard.edu/#abs/2003ApJ...585.1038S/abstract">Seager & Mallén-Ornelas (2003)</a>,
La geometría del tránsito puede ser descrita por (se presentan imágenes del trabajo de Seager & Mallén-Ornelas):
<img src="transitshape.png">
y la duración total como:
<img src="transittotaltime.png">
<img src="fig1.png">
<em>a</em> es el semieje mayor.
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:
<ol>
<li><strong>La razón entre el radio planetario y el estelar</strong> $$R_P/R_* = \sqrt{\Delta F}$$ directamente de definir $\Delta F \equiv (F_{no transit}-F_{transit})/F_{no transit} = (R_P/R*)^2$</li>
<li><strong>El parámetro de impacto</strong> 
<img src="b.png">
definido como la distancia proyectada entre el los centros del planeta y de la estrella durante la mitad del tránsito, en términos de $R_*$. Se deriva directamente de la ecuación anterior y de la ecuación de la forma del tránsito.</li>
<li><strong>La razón entre el semieje mayor de la órbita y el radio estelar</strong>
<img src="ar.png">
a partir de la ecuación de la duración del tránsito.</li>
<li><strong>La densidad estelar</strong>
<img src="rho.png">
a partir de la ecuación anterior y de la Tercera Ley de Kepler cuando $M_P \ll M_*$. Depende del parámetro de impacto puesto que $b$ afecta la duración del tránsito.
</ol>
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.
<hr>
Pregunta 2
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 <em>median filter</em> 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.
End of explanation
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'
Explanation: <hr>
Pregunta 3
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.
End of explanation
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
import batman
params = batman.TransitParams()       #object to store transit parameters
params.t0 = 0.2                     #time of inferior conjunction
params.per = p                       #orbital period
params.rp = ratio                     #planet radius (in units of stellar radii)
params.a = 180                 #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(-720, 720, 6000)  #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
Explanation: <hr>
Pregunta 4
Se calculan los coeficientes según se detalla en en <a href="https://arxiv.org/pdf/1503.07020v3.pdf">Espinoza & Jordán 2015</a>, a partir de una versión modificada del código disponible en <a href="https://github.com/nespinoza/limb-darkening">https://github.com/nespinoza/limb-darkening</a>.
End of explanation
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()
Explanation: Se ejecuta batman como se explica en la documentación, entregando como parámetros los valores obtenidos a lo largo de este trabajo.
<hr>
Pregunta Bonus
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.
End of explanation
<END_TASK>