Fitting Program

First, we were to incorporate limb-darkening into our transit model code. The figure below plots a transit without limb-darkening on the top and with limb-darkening on the bottom. A link to the code is given below, as well. You will have to download the code once you travel to the following link.

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

Second, we were to create a program that fits the transit data for HD 209458b. I used a brute force method to fit the transit data. After looking up values for the host star's radius, the semi-major axis of the planets orbit, and the inclination of the orbital plane, my python code cyclically tries values for the period, radius of the planet, and time of periastron while calculating a chi-squared statistic for each combination of parameters. It records the values for period, planet radius, and time of periastron that minimize the chi-squared statistic and writes them to a file called HD209458b_parameters.txt.

I had problems with my limb-darkening implementation with my fitting program so in an effort to get reasonable results, limb-darkening was not taken into account during the fitting procedure.

The following is a link to the code (fitting_transit.py), the input data (100binned.txt), a holistic plot of the fit (HD209458b_fit.png), a zoomed in plot of the fit of one transit (HD209458b_close_fit.png), and a file with the final best-fit parameters for the chi-squared value, the time of periastron, the planet radius, and the period of the orbit (HD209458b_parameters.txt).

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

The resultant fit for one transit is shown below.

The values that my fitting program calculate are listed below, as well as the literature values provided on the Extrasolar Planet Encyclopedia (http://exoplanet.eu/catalog/hd_209458_b/).

My values
Time_periastron = 2459490.74074 BJD

Period = 3.51863425926 days
Radius of planet = 1.34985507246 Jupiter radii

Literature values
Time_periastron = 2452968.399 JD
Period = 3.52472 days
Radius of planet = 1.38 Jupiter radii

Comparing these, we see that the period is comparable within %1 and the radius is comparable within %3. However, the time of periastron is much less reliable where our value is incorrect by over ~6,500 days. If this difference were some integer multiple of the period then the discrepancy could be explained, however, it is not. To obtain more reliable results for our parameters, we would need to develop a better method for determining a "best fit" as well as correctly incorporating limb-darkening (ie: MCMC, etc.). 

HW3 #6

Transit generating code as extended to allow surface brightness to vary across the disk of the star, and output the results as a CSV.:
https://github.com/pdn4kd/freezing-tyrion/blob/master/HW3-6.py

Example Jupiter/Sun transits with varying limb darkening parameters (no brightness variation, variation consistent with a Lambertian radiator, and an 80/20 split between the two):

---

This is a very simple (brute force) fitting demonstration for transit data from HD 209458 b. At present, it searches the entire specified parameter space at the user-specified resolution, and returns the results with the lowest reduced chi-squared. The only difference in transit calculations between this and the previous portion are that epoch is assumed to be at middle eclipse time, and a correction is included to prevent the planet from eclipsing the star when the star should be in front of the planet:
https://github.com/pdn4kd/freezing-tyrion/blob/master/HW3-6b.py

My use of unicode in the final graph means that Python 3 may be required for the plots to display correctly. Each individual test seems to take somewhat less than 1 second, so be aware of the large number of minutes/hours possible if a large parameter space is used.

example plot for HD 209458 b using 10x binning:
min chisquared:  0.0193337980195
r_planet (Jupiters):  1.4
p_planet (days):  3.52
sma (AU): 0.0471784512364
b (stellar radii):  0.5
Epoch:  2453235.5
Depth, duration, i: 0.123408960809 0.0357690898906 0.87758256189

And 100x binning:
min chisquared:  0.00087832414439
r_planet (Jupiters):  1.4
p_planet (days):  3.52
sma (AU): 0.0471784512364
b (stellar radii):  0.5
Epoch:  2453235.5
Depth, duration, i: 0.123408960809 0.0357690898906 0.87758256189

(unbinned data did not show the transits!)

Transit Fitting Program

Link

Run fit_transit.pro with no parameters. Utilizes mpfitfun.

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.