HW 3 #1-5

1)
All 3 algorithms used the same step size (0.0000243) and range (0.37026 to 324.22235) for this comparison.

Box-Least-Squares

 

It is not clear that there is any signal of a transiting exoplanet here, with the peaks not standing out from the noise. (A peak power of ~25 does not seem meaningful when the average is ~7 to 15, though it corresponds to a peak in both others)

 

Lomb-Scargle

 

This algorithm shows what seem to be actual signals, notably the one at ~4.2 days. My p-value calculations are ambiguous, with the most significant 3 tending to hit zero with all the software I tried. (eg: libreoffice)

Plavchan

 

This algorithm also shows what appear to be actual signals, though with a different pattern. Again, there is a large peak around 4.2 days. The smaller ones to to have periods that are multiples of it. Outside of the universal step-size and interval tunings, this algorithm also lets one focus on the most extreme points as they are the most likely to have actual planets.

 

2)
Using systemic, I get a planet at 4.2308 days with a mass of 0.47814 Mj and eccentricity of 0.01409. This is close to exoplanets.org's values, with differences ranging from rounding error/displayed digits to almost 10%. These values are at/within the 1-sigma errors at exoplanet.org, however. (4.230785 ± 3.6×10-5 days, 0.461 ± 0.0164 Mj*sin(i), eccentricity 0.013 ± 0.012)

The 4.23 day signal has a p-value of 1.58e-179, with 4 of the next 5 signals near it. The residual periodogram generally shows p-values around 1, with the only exceptions being relatively suspicious (1 day, and 12 synodic lunar cycles).

3)


4)
It is not clear that true anomaly actually only goes between 0 and 2 pi radians, but it should have no effect on orbital distance. Link to code: https://github.com/pdn4kd/freezing-tyrion/blob/master/HW3-4.py

5)
Solving the error propagation formula in general:

For power laws, this means:

Plugging in the transit depth formula:

Plugging in the semi-amplitude formula from in in-class:

For a typical star (∂R/R and ∂M/M ~ 0.1), this implies errors of ~10% in radius and ~7% in mass.

Homework 3

https://drive.google.com/open?id=0By4_OxQdf-GpMEcyWnNRX2pYRTg

Homework 3

Link to Drive

Homework II: Problem 5

import matplotlib.pyplot as plt
import math
import numpy as np

Ts = 5493                                                             #temperature of the star
sigma = 5.67*10E-8                                             #stefan-boltzmann constant
rp = 1.65*69E6                                                    #radius of the planet
rs = 0.872*696E6                                                #radius of the star
Ls = sigma*(Ts**4)*4*math.pi*(rs**2)                #luminosity of the star
Mstar = 0.764*1.989*10E30                               # mass of the star
a = 0.0769*1.4960E11                                        #semi-major axis
G = 6.674*10E-11                                                #gravitational constant
P = ((4*((math.pi)**2)*a**3)/(G*Mstar))**(0.5)   #period of orbit in seconds
w = (2*math.pi)/P                                                #orbital frequency
t_steps = 5000                                                     #number of time steps
s_steps = 500                                                       #number of spacial steps
L = np.zeros(t_steps)                                           #Luminosity array
x = np.zeros(t_steps)                                            #x array
t = np.arange (0, P, P/t_steps)                              #time array

length = 2*rs + 4*rp                                            #length of the box
area = np.zeros(t_steps)                                      #array of area

for i in range (0, t_steps):
        yp = 0
        xp = a*math.cos(w*t[i])                              #position of the planet at time t[i]
        x[i]=xp
        if abs(xp) > (rp + rs):
                L[i] = Ls
        else:
                countstar = 0                                      #counts inside the star
                countplanet = 0                                  #counts inside the planet
                for j in range (0, (s_steps)-1):
                        for k in range (0, (s_steps)-1):
                                xtestj = (j*length) / s_steps - length/2.0
                                ytestk = (k*length) / s_steps - length/2.0
                                ls = ((xtestj**2) + (ytestk**2))**(0.5)
                                if ls < rs:
                                        countstar = countstar + 1
                                        lp = (((xtestj - xp)**2) + ((ytestk - yp)**2))**(0.5)
                                        if lp < rp:
                                                countplanet = countplanet + 1
                area[i] = 1.0*countplanet/ countstar
                L[i] = (1-area[i])*Ls
plt.plot (t, L)
plt.xlabel('t')
plt.ylabel('L')
plt.title('Kepler 447b Transit with i=90 and e=0')
plt.show()

This code is to model the transit of Kepler 447b (parameters for this system were taken from wikipedia-https://en.wikipedia.org/wiki/List_of_exoplanets). 

Of course, this python code is not as robust as it could be. It assumes (among other things) that the inclination, i, is 90 degrees, and that the eccentricity, e, is zero. Also, it does not accurately model the secondary transit.

This code does take into account the luminosity of the host star, the host star and planet radii, the host star mass, the semi-major axis of the planet, the period of the planet's orbit, the impact parameter equal to zero, and it models the ingress and egress of the transit. 

Assignment 2 #5

LINK TO CODE