Homework 4

1)
6.

7. Earth does not appear to be at any risk of a particular catastrophe. Semi-major axis is randomly moving between the minimum detectable changes, and there is a slow tradeoff in parameters from eccentricity to inclination.

2)
Suitably apocalyptic results are obtained when adding a planet with the parameters of :

NIBIRU     m=1.40000000000000000E-02 r=3.d0 d=28.0
3.14159265358979323E+01 -2.71828182845904528E+00 -1.01001000100001000E-01
7.36084560000000000E-04  7.43519750000000000E-07 -2.00460016972063023E-11
0. 0. 0.

A number of them simply disappear from the simulation after a bit.

3)
a. The two properties of a protoplanetary disk that determine what planets will form are mass, and metallicity.
*b. The three types of planetary migration are named Type I, II, and III. Type I are from the planet causing overdensities in the disk inside and outside of its orbit. The density profile of the disk results in the planet transferring angular momentum to the overdensities, and moving inwards. In Type II migration, the planet is sufficiently massive to open up a gap in the disk. (something about viscous forces and drag). Type III is runaway...  Seager(?) describes the timescales for Type I and II as being much the same.
c. Terrestrial planets can survive giant planet migration (including the formation of hot Jupiters) by outwards scattering, or (more likely) forming exterior to the giant planet and migrating inwards.

4)
a. An "inversion layer" is a layer of a planet's atmosphere where the temperature profile is inverted; instead of decreasing, the temperature increases with altitude. This is important becauseof how it affects wind patterns and heat flow, with atmospheres partially stratifying into layers based on temperature inversions.
b.
Peverest = exp(-8.8/8) ~= 0.33287108369807955
Alternatively:
ln(Peverest) = -alt/scale
scale*8.8/8 = alt

For an isothermal atmosphere, scale = T*k/(M*g)
k = 1.38065582e-23

M (molecular weight of the atmosphere) is assumed to be the same as for Earth. Approximating as 80% Nitrogen, 20% Oxygen, M = (0.8*28+0.2*32)/6.022e26 = 4.782464297575557e-26 kg

(exoplanet properties: a = 0.2 AU, m = 5 M_earth, r = 1.2 R_earth, P_0 = 1.0 atm)
Simplifying the Stephan-Boltzmann equation for blackbody equilibrium:
T^4 = c/a^2; T = (c/a^2)^0.25
Earth: a = 1 AU, T = 287 K, c = 300^4
This planet:
T = 287/sqrt(0.2) = 641.75150954243963 K
g = 9.81*5/1.2^2 = 34.0625 m/s/s

scale = 641.75150954243963*1.38064852e-23/(28.8/6.022e26*34.0625) = 5439.03400930283 m
Checking scale height with assumed Earth parameters: 287*1.38064852e-23/((28*0.8+32*0.2)/6.022e26*9.81) = 8445.867900509964 m, a bit more than 5% larger than the given 8 km.

5439.03400930283*8.8/8 = 5982.937410233113 m
At the actual level of precision, ~6.0 km altitude.

5)
a. The primary cause of gaps in the asteroid belt is resonances with Jupiter causing the asteroids to be perturbed into increasingly eccentric orbits. The ultimate result is a close encounter (or coll
ision!) with another planet dramatically altering the asteroid's orbit.
b. The most surprising features of early exoplanet discoveries were the planets existing at all where they were found. Neither pulsars, nor <10 day orbits around FGK stars were expected.
c. Exoplanet letter designations are assigned in order of discovery, with ties being broken by whichever is closest to its parent star.
d. Kepler has demonstrated that smaller planets are more common than larger ones, at least down to where it begins to suffer from poor completeness (in the sub-neptune/super-earth range).
e. The "habitable" zone is the range of distances from a star where we expect that at least some kinds of rocky planets could have liquid water on their surfaces.

6)
a. The maximum phase amplitude will be just before secondary eclipse. The angle will be:
90-arccos(695700e3/0.723332/1459597870700) = 0.36837 degrees.
Lighting is ambiguous, and assumed to be cos(theta/2)**2 since that gives a plausible shape to the light curve. This gives a phase of 0.99998966614598372. (about 1 - 10^-6)
b. With this inclined orbit, maximum illumination will have the planet 40 degrees from conjuction. This is conveniently also peristron. Using the above guess at illumination, the phase is 0.88302222155948906. Apparent brightness at the same phase angle is proportional to 1/r^2, so the relative distance is:
sqrt(0.88302222155948906/0.99998966614598372) = 0.93969747614672106
0.9397 is the periastron for an orbit with a semi-major axis of 1.
p = a*(1-e)
e = 1-(p/a)
e = 1-(0.93969747614672106/1) = 0.060302523853278944
e ~= 0.06
(Venus' eccentricity is approximated as 0)

7)
a. Mean anomaly is where the body would be along its orbit if it were in a 0 eccentricity orbit. True anomaly is where along the orbit the body actually is. Eccentric anomaly is an auxiliary angle to facilitate converting between the easily found mean anomaly and more difficult true anomaly.

b.

Given the difficulties in hand-drawn diagrams, perhaps electronic ones would be better?

Homework 4

Link

Final Project: Confidence intervals for the KS-test of KIC 11558725

Here is the link to my dual final project for Dr. Sakidja and Dr. Plavchan.

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




Have a great Christmas break!

Final Homework

Here is the link to my final homework for AST 613.

https://drive.google.com/file/d/0By4_OxQdf-GpdG1udlctMmxLN1U/view




Merry Christmas!

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. 

Assignment 2 #5

LINK TO CODE

Homework 2

1)
As the planet mass is ~13% that of the star, it matters considerably.
velocity:

r is a fraction of the semi-major axis inversely proportional to relative mass:

Plugging in the numbers, V is 98.75 m/s. A speed of almost 100 m/s seems quite feasible.


2)
The planet orbiting VB10 was discovered by looking for periodicities in astrometric data. It was refuted by the suggested parameters (and a large portion of the physically plausible parameters) being ruled out by improvements in the radial velocity technique. The RV refutation was helped by the astrometric orbital parameters indicating a near edge-on inclination.

3)
a.
For the astrometric method, planet mass (as a fraction of stellar mass) shows an inverse relationship with semi-major axis. For the radial velocity method, it increases with the square root of stellar mass and semi-major axis.

b.
The cross-over depends strongly on stellar mass and distance, ranging from ~1.5e11 m (1 AU) for the plausible near/low mass case to ~3e12 m (20 AU) for a more distant and likely overly high mass star. (The upper left portion of the second graph should probably be considered unphysical)

4)
 Starting with the Planck equation, a spherical star, and a spherical grain in thermal equilibrium (like the in-class derivation, but with variable opacity):

Evaluating the integral in Wolfram Alpha yields:

Substituting:

(The equation given has 2 instead of sqrt(2), but I think it was described as a typesetting error?)

5)
https://github.com/pdn4kd/freezing-tyrion/blob/master/HW2-5.py

Homework 2

Link to assigment 2

Assignment 2

Link to Assignment 2

HW1

1)
a.
Estimating the maximum angular separation of Earth and Jupiter from the Sun at a distance of 10     parsecs: Earth's Aphelion is 1.01673 AU, Jupiter's is 5.45492 AU.1 AU at 1 pc is 1", so scaling proportionally, these are (to excessive significant figures) 0.101673" and 0.545492" from the Sun. Earth is questionable with the nominal HST resolution of 0.05" (the star and planet would be adjacent), but Jupiter (some 20 pixels away) could be distinguished.
b.
Using the planet-to-star flux ratio approximation R_p²/a²,

Planet	Radius (m)	Semi-Major Axis (m)	Flux Ratio	Magnitudes Fainter
Earth	6.371e6	1.49598023e11	1.814e-9	21.9
Jupiter	6.9911e7	7.78299e11	8.069e-9	20.2

2)
The problems of defining a planet at the high mass end are how they become increasingly star-like. Newly formed brown dwarfs are notable for having effective temperatures (few thousand kelvin), compositions (hydrogen-rich), and densities very close to still-forming red dwarfs. Possible solutions include looking at the spectra for deuterium, helium-3, and lithium abundances (higher in brown dwarfs than stars, at least once the stars have had time to burn those elements/isotopes).

3)
a.
Using Kepler's 3rd law in its full Newtonian glory:
T² == 4*pi²/(G*(M1+M2))*a³, or T == 2*pi*a^1.5/sqrt(G(M1+M2)).
For the stars (a == 1200 AU, M1 == M2 == 0.9 M_sol), T == 9.78e11 s, or 31000 (Julian) years.
b.
For the planet, T == 1.004e7 s, or 0.318 (Julian years), so Tstar/Tplanet ~= 97400 orbits

4)
Using Kepler's 3rd law in simplified form:
T² == a³,  a == T^2/3 == (120/365.25)^2/3 == 0.47613 AU
Periastron == a(1-e)
e == 1 - Periastron/a == 0.36992 ~= 0.37

5)

 I would not expect the decreasing density vague power law guess (mass ~ radius^2.5) to hold for higher masses (topping out at perhaps 2x Jupiter's radius even with >10x Jupiter's mass), given predictions for brown dwarf sizes and the relatively small radii of low mass red dwarfs.

Homework 1

https://drive.google.com/drive/folders/0B23b_IF-Qh6TMmxzaGpub1ZTQXc

Homework 1

1a) If the distance (d) to the Sun from Earth is ~ 1 AU and the distance to the Sun from Jupiter is ~ 5.3 AU, then the maximum angular separation of both planets from a distance (D) of 10 parsecs is given by the equation,

θ = d/D.

Plugging in the correct values for Earth and Jupiter we get,

θ_Earth= 4.84813681 × 10^-7 radians = 0.1” (arcsecs)

and

θ_Jupiter= 2.56951251 × 10^-6 radians = 0.53” (arcsecs)

Considering the stated resolution for the Hubble Space Telescope is 0.05” (arcsecs), we can resolve the separation of either planet from the Sun at a distance of 10 parsecs away.

1b) If the radius (r) of Earth is 6.371 x 10^6 meters and the radius of Jupiter is 69.911 x 10^6 meters, then the planet-to-star flux ratio is approximately given by

flux ratio ≈ r^2 / d^2.

So plugging in the correct values for Earth and Jupiter we get,

Earth-Sun flux ratio ≈ 1.81774025 × 10^-9

and

Jupiter-Sun flux ratio ≈ 2.28383156 × 10^-7

2) Problems defining the planets on the high-mass end is, in part, due to the over-lap of high-mass planets and low-mass brown dwarfs. Low-end properties of mass and temperature for brown dwarfs have been observed to be ~13 Jupiter masses and ~300K-750K. Also, the high-end properties of mass and temperature for exoplanets have been observed to be ~28 Jupiter masses and ~7500K. From this, it is obvious that there is over-lap between the characteristics of brown dwarfs and exoplanets. Direct spectrometry of the objects could help clarify the line between these two. However, this is difficult considering, from 1a, the resolving power of the HST is just small enough to resolve a planet at 1 AU from its host star. Even though the HST has the angular resolution to resolve this planet perhaps, with a flux ratio (from 1b) as small as ~10^-9, the planet's reflected light would still be washed out by the light from its host star.

3a) The orbital period of a binary star system is given by

So plugging in the correct values for HD 80606 and HD 80607 and solving for T we get,

T ≈ 30,984 years

3b) The orbital period for a planet-star system is given by the same formula as above where M1 >> M2. So plugging in the correct values for HD 80606 and HD 80606b we get,

T ≈ 0.318 years

Meaning that HD 80606b makes ~97,433 orbits for every one binary orbit of
HD 80606 and HD 80607.

4) If the orbital period is known then we can solve for the length of the semi-major axis using,

Plugging in the orbital period of 120 days and using a M1 = 1 solar mass >> M2 we get,

a ≈ 0.48 AU

Then, knowing that the closest approach of the planet is 0.3 AU we can solve for the eccentricity using

Radius of periastron = a(1-e)

Plugging in the values for the radius of periastron and the semi-major axis we get that

e ≈ 0.375.

5)

Assuming that the relationship between mass and radius is a power law given by

m^n = r ; n is any real number

Then the average of n for all the planets is

n ≈ 0.162.

Therefore the power law relationship between mass and radius is

m^0.162 ≈ r.

Planets larger than Jupiter would follow this relationship until their densities become so high that electron degeneracy pressures would overcome the classical concepts of pressure and temperature as they relate to ideal gases.