Friday, December 16, 2016

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?

Wednesday, November 30, 2016

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.). 

Tuesday, November 29, 2016

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!)

Wednesday, November 16, 2016