In [1]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
# Import dev version of friedrich:
import sys
sys.path.insert(0, '../')
from friedrich.lightcurve import hat11_params_morris
from friedrich.orientation import planet_position_cartesian
transit_params = hat11_params_morris()
start = 0.5
n_times = 1000
interp_duration = 0.005 # fraction of full transit duration
times1 = np.linspace(transit_params.t0 - transit_params.duration/2 - interp_duration,
transit_params.t0 - transit_params.duration/2 + interp_duration, n_times)
times2 = np.linspace(transit_params.t0 + transit_params.duration/2 - interp_duration,
transit_params.t0 + transit_params.duration/2 + interp_duration, n_times)
times = np.concatenate([times1, times2])
x_center, y_center, z_center = planet_position_cartesian(times, transit_params)
y_plus = y_center + transit_params.rp
y_minus = y_center - transit_params.rp
theta = np.linspace(0, 2*np.pi, 50)
fig, ax = plt.subplots(figsize=(10,10))
ax.plot(np.cos(theta), np.sin(theta), color='k')
ax.plot(x_center, y_center, color='r')
#ax.scatter(x, y)
ax.fill_between(x_center, y_minus, y_plus, color='k', alpha=0.3)
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.set_aspect('equal')
To first order the planet occults a solid band of latitudes in sky-projected coordinates, if we neglect the orbital eccentricity and inclination. Here we're not concerned with the latitudes and longitudes on the stellar surface with respect to the stellar rotation axis, so I refer to latitudes as bands on the stellar surface parallel with the sky-projected x axis.
The limb of the star is defined by the unit circle in sky coordinates $$ x^2 + y^2 = 1 .$$ The cross-sectional radius of the star $\rho$ at a given $y$ is $$ \rho(y) = \sqrt{1 - y^2}.$$ $$ \rho^\prime(y) = y(1 - y^2)^{-1/2}$$ The circumference of the star $C$ at a given $y$ is $$ C(y) = 2\pi\rho(y). $$
We can find the surface area of segments of the star by using the surface of revolution defined by revolving a function $f(y)$ (in this case the unit circle) around the $y$ axis, which has solutions of the form (see derivation) $$ A = \int 2 \pi f(y) \sqrt{ 1 + f^\prime(y)^2 } ~dy. $$
The transit chord approximately extends in $y$ from $b - R_p/R_s$ to $b + R_p/R_s$, and covers only the hemisphere of the star that is visible to the observer (50% of the total area), so the surface area of the star below the transit chord is: $$ A = \int_{b - R_p/R_s}^{b + R_p/R_s} \frac{1}{2} ~ 2\pi \rho(y) \sqrt{ 1 + \rho^\prime(y)^2 } ~dy ~~=~~ \int_{b - R_p/R_s}^{b + R_p/R_s} \pi \rho(y) \sqrt{ 1 + \rho^\prime(y)^2} ~dy, $$
which we can solve numerically:
In [2]:
from scipy.integrate import quad
def rho(y):
"""Radius at a given latitude (latitude = y)"""
return np.sqrt(1 - y**2)
def rho_prime(y):
"""d(rho)/dy"""
return y / np.sqrt(1 - y**2)
def integrand(y):
return np.pi * rho(y) * np.sqrt(1 + rho_prime(y)**2)
y_lower_bound = transit_params.b - transit_params.rp
y_upper_bound = transit_params.b + transit_params.rp
# Integrate the integrand function from y_lower_bound to y_upper_bound:
surface_area = quad(integrand, y_lower_bound, y_upper_bound)[0]
fractional_surface_area_per_transit = surface_area/(4*np.pi)
print("Fractional surface area occulted by one transit of HAT-P-11b: {0}"
.format(fractional_surface_area_per_transit))
n_transits = 5
print("Assuming no overlap, {0:.1f}% of the stellar surface is occulted per stellar rotation ({1} transits)"
.format(n_transits*fractional_surface_area_per_transit*100, n_transits))
In [3]:
y_lower_bound = -1
y_upper_bound = 1
surface_area = quad(integrand, y_lower_bound, y_upper_bound)[0]
print("fractional surface area occulted: {0}".format(surface_area/(4*np.pi)))
If the top half of the visible hemisphere is occulted, we can expect one quarater of the surface area to be occulted:
In [4]:
y_lower_bound = 0
y_upper_bound = 1
surface_area = quad(integrand, y_lower_bound, y_upper_bound)[0]
print("fractional surface area occulted: {0}".format(surface_area/(4*np.pi)))
In [37]:
np.unique(whole_lc.quarters)
Out[37]:
In [41]:
# Get quarters
from friedrich.lightcurve import LightCurve
light_curve_paths = glob('/local/tmp/hat11/*slc.fits')
whole_lc = LightCurve.from_raw_fits(light_curve_paths, name='HAT11')
sort_by_time = np.argsort(whole_lc.times.jd)
all_times = whole_lc.times.jd[sort_by_time]
all_quarters = whole_lc.quarters[sort_by_time]
from friedrich.analysis import MCMCResults
from glob import glob
paths = sorted(glob('/local/tmp/friedrich/hat11/chains*.hdf5'))
# Collect info from all friedrich spot fits:
spots_all_transits = []
quarters_all_transits = []
for i, path in enumerate(paths):
m = MCMCResults(path, hat11_params_morris())
spots_all_transits.append(m.get_spots())
quarters_all_transits.append(all_quarters[np.argmin(np.abs(m.lc.times.jd.mean() -
all_times))])
del m
In [43]:
# Get number of spots with delta BIC > 10
n_spots = np.zeros(len(spots_all_transits))
times = np.zeros(len(spots_all_transits))
quarters = np.zeros(len(spots_all_transits))
for i, transit, q in zip(range(len(spots_all_transits)),
spots_all_transits, quarters_all_transits):
for spot in transit:
if spot.delta_BIC > 10:
n_spots[i] += 1
if times[i] == 0:
times[i] = spot.t0.value
if quarters[i] == 0:
quarters[i] = q
In [7]:
# Plot the distribution of spot numbers:
fig, ax = plt.subplots(1, 2, figsize=(14, 7))
ax[0].plot(times, n_spots, 'o')
ax[0].set(xlabel='JD', ylabel='# spots')
ax[0].set_ylim([-0.1, max(n_spots)+0.1])
ax[1].hist(n_spots, len(np.unique(n_spots)))
ax[1].set(xlabel='# spots')
plt.show()
To estimate the number of detectable spots on the surface of the star, divide the number of spots per transit by the fractional surface area of the star occulted per transit.
In [8]:
import seaborn as sns
sns.set_style("darkgrid")
from astropy.time import Time
n_spots_whole_star = n_spots/fractional_surface_area_per_transit
fig, ax = plt.subplots(1, figsize=(10, 5))
ax.plot(times, n_spots_whole_star, 'o-')
ax.set(xlabel='JD', ylabel='# spots on whole star')
plt.show()
In [9]:
rotation_period = transit_params.per_rot
n_spots_per_rotation = []
times_per_rotation = []
for i in range(3, len(times) - 3):
within_one_rotation = ((times[i] - rotation_period/2 <= times) &
(times[i] + rotation_period/2 >= times))
n_transits_in_window = np.count_nonzero(within_one_rotation)
times_per_rotation.append(np.mean(times[within_one_rotation]))
n_spots_in_window = np.sum(n_spots[within_one_rotation])
n_spots_per_area = n_spots_in_window/(n_transits_in_window*fractional_surface_area_per_transit)
n_spots_per_rotation.append(n_spots_per_area)
times_per_rotation = np.array(times_per_rotation)
n_spots_per_rotation = np.array(n_spots_per_rotation)
fig, ax = plt.subplots(1, figsize=(10, 5))
ax.plot(times_per_rotation, n_spots_per_rotation, 'o-')
ax.set(ylabel='Total number of spots per stellar rotation',
xlabel='Time [JD]')
plt.show()
In [51]:
quarters_unique = np.sort(np.unique(quarters))
n_spots_per_quarter = np.zeros(len(quarters_unique))
n_transits_per_quarter = np.zeros(len(quarters_unique))
fake_uncertainty = np.zeros(len(quarters_unique))
for i, q in enumerate(quarters_unique):
in_this_quarter = (quarters == q)
n_spots_per_quarter[i] += np.sum(n_spots[in_this_quarter])
n_transits_per_quarter[i] = np.count_nonzero(in_this_quarter)
fake_uncertainty[i] = np.std(n_spots[in_this_quarter])/np.sqrt(n_transits_per_quarter[i])
In [65]:
fig, ax = plt.subplots(1, 2, figsize=(12, 6))
ax[0].plot(quarters_unique, n_spots_per_quarter, 'o')
ax[0].set_xlabel('Quarter')
ax[0].set_ylabel('Number of spots')
ax[0].set_xlim([-0.5, 18])
ax[1].plot(quarters_unique, n_transits_per_quarter, 'o')
ax[1].set_xlabel('Quarter')
ax[1].set_ylabel('Number of transits observed')
ax[1].set_xlim([-0.5, 18])
Out[65]:
In [76]:
params, cov = np.polyfit(quarters_unique,
n_spots_per_quarter/n_transits_per_quarter,
1, cov=True)
intercept = params[1]
intercept_uncertainty = np.sqrt(np.diag(cov))[1]
plt.errorbar(quarters_unique, n_spots_per_quarter/n_transits_per_quarter,
fake_uncertainty, fmt='o')
plt.fill_between(quarters_unique, intercept-intercept_uncertainty,
intercept+intercept_uncertainty, alpha=0.2)
plt.axhline(intercept, color='k')
plt.xlabel('Quarter')
plt.ylabel('Mean spots per transit')
plt.xlim([-0.5, 18])
print("Mean number of spots per transit: {0:.2f} +/- {1:.2f}"
.format(intercept, intercept_uncertainty))
In [ ]: