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.