Skip to content

Commit

Permalink
adding docstrings
Browse files Browse the repository at this point in the history
  • Loading branch information
@bmorris3
bmorris3 committed Oct 7, 2021
1 parent fb754eb commit 6d6912d
Show file tree
Hide file tree
Showing 6 changed files with 343 additions and 137 deletions.
80 changes: 66 additions & 14 deletions jaxon/continuum.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,27 @@
@jit
def log_hminus_continuum(wavelength_um, temperature, pressure,
volume_mixing_ratio_product, truncation_value=-100):

"""
Compute the H-minus continuum.
Parameters
----------
wavelength_um : numpy.ndarray
Wavelength in microns
temperature : numpy.ndarray
Array of temperatures
pressure : numpy.ndarray
Array of pressures
volume_mixing_ratio_product : float
Product of the VMR of the H- abundance
truncation_value : float
Truncate values below log10 of ``truncation_value``
Returns
-------
absorption_coeff : tensor-like
Log10 absorption coefficient
"""
# first, compute the cross sections (in cm4/dyne)
kappa_bf = bound_free_absorption(wavelength_um, temperature)
kappa_ff = free_free_absorption(wavelength_um, temperature)
Expand All @@ -30,6 +50,21 @@ def log_hminus_continuum(wavelength_um, temperature, pressure,

@jit
def bound_free_absorption(wavelength_um, temperature):
"""
Bound-free absorption.
Parameters
----------
wavelength_um : numpy.ndarray
Wavelength in micron
temperature : numpy.ndarray
Temperature array
Returns
-------
kappa_bf : tensor-like
Bound-free absorption
"""
# Note: alpha has a value of 1.439e4 micron-1 K-1, the value stated in John (1988) is wrong
# here, we express alpha using physical constants
alpha = CONST_C * CONST_H / CONST_K * 10000.0
Expand Down Expand Up @@ -70,6 +105,21 @@ def body_fun(val, x):

@jit
def free_free_absorption(wavelength_um, temperature):
"""
Free-free absorption via John (1988).
Parameters
----------
wavelength_um : numpy.ndarray
Wavelength in micron
temperature : numpy.ndarray
Temperature array
Returns
-------
kappa_ff : tensor-like
Free-free absorption
"""
# coefficients from John (1988)
# to follow his notation (which starts at an index of 1), the 0-index components are 0
# for wavelengths larger than 0.3645 micron
Expand Down Expand Up @@ -123,19 +173,21 @@ def body_fun(val, x):

@jit
def dtauHminusCtm(nus, Tarr, Parr, dParr, volume_mixing_ratio_product, mmw, g):
"""dtau of the H- continuum
Args:
nus: wavenumber matrix (cm-1)
Tarr: temperature array (K)
Parr: temperature array (bar)
dParr: delta temperature array (bar)
volume_mixing_ratio_product: number density for e- times number density for H [N_layer]
mmw: mean molecular weight of atmosphere
g: gravity (cm2/s)
nucia: wavenumber array for CIA
tcia: temperature array for CIA
logac: log10(absorption coefficient of CIA)
"""
dtau of the H- continuum
Parameters
----------
nus: wavenumber matrix (cm-1)
Tarr: temperature array (K)
Parr: temperature array (bar)
dParr: delta temperature array (bar)
volume_mixing_ratio_product: number density for e- times number density for H [N_layer]
mmw: mean molecular weight of atmosphere
g: gravity (cm2/s)
nucia: wavenumber array for CIA
tcia: temperature array for CIA
logac: log10(absorption coefficient of CIA)
Returns:
optical depth matrix [N_layer, N_nus]
Expand Down
54 changes: 29 additions & 25 deletions jaxon/hatp7.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,36 +32,40 @@
spitzer_depth = 1e-2 * np.array([0.161, 0.186])
spitzer_depth_err = 1e-2 * np.array([0.014, 0.008])

# plt.plot(
# central_wl, depths, 'o'
# )
# plt.plot(
# spitzer_wl, spitzer_depth, 'o'
# )

kepler_mean_wl = [0.641] # um
kepler_depth = [19e-6] # eyeballed
kepler_depth_err = [10e-6]
def get_observed_depths():
kepler_mean_wl = [0.641] # um
kepler_depth = [19e-6] # eyeballed
kepler_depth_err = [10e-6]

all_depths = np.concatenate([depths, spitzer_depth])
all_depths_errs = np.concatenate([depths_err, spitzer_depth_err])
all_wavelengths = np.concatenate([central_wl.value, spitzer_wl])
return all_depths, all_depths_errs, all_wavelengths, kepler_mean_wl

all_depths = np.concatenate([depths, spitzer_depth])
all_depths_errs = np.concatenate([depths_err, spitzer_depth_err])
all_wavelengths = np.concatenate([central_wl.value, spitzer_wl])

rprs = float(1.431*u.R_jup / (2.00 * u.R_sun))

g = (G * u.M_jup / u.R_jup**2).to(u.cm/u.s**2).value

t0 = 2454954.357462 # Bonomo 2017
period = 2.204740 # Stassun 2017
rp = 16.9 * u.R_earth # Stassun 2017
rstar = 1.991 * u.R_sun # Berger 2017
a = 4.13 * rstar # Stassun 2017
duration = 4.0398 / 24 # Holczer 2016
b = 0.4960 # Esteves 2015
rho_star = 0.27 * u.g / u.cm ** 3 # Stassun 2017
T_s = 6449 # Berger 2018

a_rs = float(a / rstar)
a_rp = float(a / rp)
rp_rstar = float(rp / rstar)
eclipse_half_dur = duration / period / 2
def get_planet_params():
t0 = 2454954.357462 # Bonomo 2017
period = 2.204740 # Stassun 2017
rp = 16.9 * u.R_earth # Stassun 2017
rstar = 1.991 * u.R_sun # Berger 2017
a = 4.13 * rstar # Stassun 2017
duration = 4.0398 / 24 # Holczer 2016
b = 0.4960 # Esteves 2015
rho_star = 0.27 * u.g / u.cm ** 3 # Stassun 2017
T_s = 6449 # Berger 2018

a_rs = float(a / rstar)
a_rp = float(a / rp)
rp_rstar = float(rp / rstar)
eclipse_half_dur = duration / period / 2

return (
planet_name, a_rs, a_rp, T_s, rprs, t0, period, eclipse_half_dur, b,
rstar, rho_star, rp_rstar
)
166 changes: 105 additions & 61 deletions jaxon/lightcurve.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,70 +8,114 @@
import pymc3_ext as pmx

from .utils import floatX
from .hatp7 import (
planet_name, t0, period, eclipse_half_dur, b, rstar, rho_star, rp_rstar
)
from .hatp7 import get_planet_params

__all__ = [
'eclipse_model'
'eclipse_model',
'get_filter',
'get_light_curve'
]

lcf = search_lightcurve(
planet_name, mission="Kepler", cadence="long"
# quarter=10
).download_all()

slc = lcf.stitch()

phases = ((slc.time.jd - t0) % period) / period
in_eclipse = np.abs(phases - 0.5) < 1.5 * eclipse_half_dur
in_transit = (
(phases < 1.5 * eclipse_half_dur) |
(phases > 1 - 1.5 * eclipse_half_dur)
)
out_of_transit = np.logical_not(in_transit)

slc = slc.flatten(
polyorder=3, break_tolerance=10, window_length=1001, mask=~out_of_transit
).remove_nans()

phases = ((slc.time.jd - t0) % period) / period
in_eclipse = np.abs(phases - 0.5) < 1.5 * eclipse_half_dur
in_transit = (
(phases < 1.5 * eclipse_half_dur) |
(phases > 1 - 1.5 * eclipse_half_dur)
)
out_of_transit = np.logical_not(in_transit)

sc = sigma_clip(
np.ascontiguousarray(slc.flux[out_of_transit], dtype=floatX),
maxiters=100, sigma=8, stdfunc=mad_std
)

phase = np.ascontiguousarray(
phases[out_of_transit][~sc.mask], dtype=floatX
)
time = np.ascontiguousarray(
slc.time.jd[out_of_transit][~sc.mask], dtype=floatX
)

bin_in_eclipse = np.abs(phase - 0.5) < eclipse_half_dur
unbinned_flux_mean = np.mean(sc[~sc.mask].data)

unbinned_flux_mean_ppm = 1e6 * (unbinned_flux_mean - 1)
flux_normed = np.ascontiguousarray(
1e6 * (sc[~sc.mask].data / unbinned_flux_mean - 1.0), dtype=floatX
)
flux_normed_err = np.ascontiguousarray(
1e6 * slc.flux_err[out_of_transit][~sc.mask].value, dtype=floatX
)

filt = Filter.from_name("Kepler")
filt.bin_down(4) # This speeds up integration by orders of magnitude
filt_wavelength, filt_trans = filt.wavelength.to(u.m).value, filt.transmittance


def eclipse_model():
cadence = "long"
cadence_duration = 30 * u.min


def get_light_curve(cadence=cadence):
"""
Parameters
----------
cadence : str {'long', 'short'}
Kepler cadence mode
"""
(planet_name, a_rs, a_rp, T_s, rprs, t0, period, eclipse_half_dur, b,
rstar, rho_star, rp_rstar) = get_planet_params()

lcf = search_lightcurve(
planet_name, mission="Kepler", cadence=cadence
# quarter=10
).download_all()

slc = lcf.stitch()

phases = ((slc.time.jd - t0) % period) / period
in_transit = (
(phases < 1.5 * eclipse_half_dur) |
(phases > 1 - 1.5 * eclipse_half_dur)
)
out_of_transit = np.logical_not(in_transit)

slc = slc.flatten(
polyorder=3, break_tolerance=10, window_length=1001, mask=~out_of_transit
).remove_nans()

phases = ((slc.time.jd - t0) % period) / period
in_transit = (
(phases < 1.5 * eclipse_half_dur) |
(phases > 1 - 1.5 * eclipse_half_dur)
)
out_of_transit = np.logical_not(in_transit)

sc = sigma_clip(
np.ascontiguousarray(slc.flux[out_of_transit], dtype=floatX),
maxiters=100, sigma=8, stdfunc=mad_std
)

phase = np.ascontiguousarray(
phases[out_of_transit][~sc.mask], dtype=floatX
)
time = np.ascontiguousarray(
slc.time.jd[out_of_transit][~sc.mask], dtype=floatX
)

bin_in_eclipse = np.abs(phase - 0.5) < eclipse_half_dur
unbinned_flux_mean = np.mean(sc[~sc.mask].data)

unbinned_flux_mean_ppm = 1e6 * (unbinned_flux_mean - 1)
flux_normed = np.ascontiguousarray(
1e6 * (sc[~sc.mask].data / unbinned_flux_mean - 1.0), dtype=floatX
)
flux_normed_err = np.ascontiguousarray(
1e6 * slc.flux_err[out_of_transit][~sc.mask].value, dtype=floatX
)
return phase, time, flux_normed, flux_normed_err


def get_filter():
"""
Get the Kepler bandpass filter transmittance curve
Returns
-------
filt_wavelength : numpy.ndarray
Wavelengths in the transmittance curve
filt_trans : numpy.ndarray
Transmittances
"""
filt = Filter.from_name("Kepler")
filt.bin_down(4) # This speeds up integration by orders of magnitude
filt_wavelength, filt_trans = filt.wavelength.to(u.m).value, filt.transmittance
return filt_wavelength, filt_trans


def eclipse_model(cadence_duration=cadence_duration):
"""
Compute the (static) eclipse model
Parameters
----------
cadence_duration : astropy.unit.Quantity
Exposure duration
Return
------
eclipse_numpy : numpy.ndarray
Occultation vector normalized to unity out-of-eclipse and zero
in-eclipse.
"""
(planet_name, a_rs, a_rp, T_s, rprs, t0, period, eclipse_half_dur, b,
rstar, rho_star, rp_rstar) = get_planet_params()
phase, time, flux_normed, flux_normed_err = get_light_curve()

with pm.Model():
# Define a Keplerian orbit using `exoplanet`:
orbit = xo.orbits.KeplerianOrbit(
Expand All @@ -83,7 +127,7 @@ def eclipse_model():
eclipse_light_curves = xo.LimbDarkLightCurve([0, 0]).get_light_curve(
orbit=orbit._flip(rp_rstar), r=orbit.r_star,
t=phase * period,
texp=(30 * u.min).to(u.d).value
texp=cadence_duration.to(u.d).value
)

# Normalize the eclipse model to unity out of eclipse and
Expand Down

0 comments on commit 6d6912d

Please sign in to comment.