# -*- coding: utf-8 -*-
import astropy.units as u
import numpy as np
import scipy
KM2AU = (1 * u.km).to_value(u.AU)
[docs]
def deltaMag(p, Rp, d, Phi):
"""Calculates delta magnitudes for a set of planets, based on their albedo,
radius, and position with respect to host star.
Args:
p (ndarray):
Planet albedo
Rp (astropy Quantity array):
Planet radius in units of km
d (astropy Quantity array):
Planet-star distance in units of AU
Phi (ndarray):
Planet phase function
Returns:
~numpy.ndarray:
Planet delta magnitudes
"""
dMag = -2.5 * np.log10(
p * (Rp.to_value(u.km) * KM2AU / d.to_value(u.AU)) ** 2 * Phi
)
return dMag
[docs]
def betaStar_Lambert():
"""Compute the Lambert phase function deltaMag-maximizing phase angle
Args:
None
Returns:
float:
Value of beta^* in radians
"""
betastarexpr = (
lambda beta: -(np.pi - beta) * np.sin(beta) ** 3 / np.pi
+ 2
* ((np.pi - beta) * np.cos(beta) + np.sin(beta))
* np.sin(beta)
* np.cos(beta)
/ np.pi
)
betastar = scipy.optimize.fsolve(betastarexpr, 63 * np.pi / 180)[0]
return betastar
[docs]
def min_deltaMag_Lambert(Completeness, s=None):
"""Calculate the minimum deltaMag at given separation(s) assuming a Lambert phase
function
Args:
Completeness (BrownCompleteness):
BrownCompleteness object
s (float or ~numpy.ndarray, optional):
Projected separations (in AU) to compute minimum delta mag at.
If None (default) then uses Completeness.xnew
Returns:
~numpy.ndarray:
Minimum deltaMag values
"""
# this is the output of betaStar_Lambert (no need to recompute every time)
betastar = 1.1047288186445432
# if no s input supplied, use the full array of separations from Completeness
if s is None:
s = Completeness.xnew
# allocate output
dmagmin = np.zeros(s.size)
# idenitfy breakpoints
Ppop = Completeness.PlanetPopulation
bp1 = Ppop.rrange.min().to(u.AU).value * np.sin(betastar)
bp2 = Ppop.rrange.max().to(u.AU).value * np.sin(betastar)
# compute minimum delta mags
dmagmin[s < bp1] = -2.5 * np.log10(
Ppop.prange.max()
* ((Ppop.Rprange.max() / Ppop.rrange.min()).decompose().value) ** 2
* Completeness.PlanetPhysicalModel.calc_Phi(
(np.arcsin(s[s < bp1] / Ppop.rrange.min().value)) * u.rad
)
)
inds = (s >= bp1) & (s < bp2)
dmagmin[inds] = -2.5 * np.log10(
Ppop.prange.max()
* ((Ppop.Rprange.max().to(u.AU).value / s[inds])) ** 2
* Completeness.PlanetPhysicalModel.calc_Phi(betastar * u.rad)
* np.sin(betastar) ** 2
)
dmagmin[s >= bp2] = -2.5 * np.log10(
Ppop.prange.max()
* ((Ppop.Rprange.max() / Ppop.rrange.max()).decompose().value) ** 2
* Completeness.PlanetPhysicalModel.calc_Phi(
(np.arcsin(s[s >= bp2] / Ppop.rrange.max().value)) * u.rad
)
)
return dmagmin
[docs]
def max_deltaMag_Lambert(Completeness, s=None):
"""Calculate the maximum deltaMag at given separation(s) assuming a Lambert phase
function
Args:
Completeness (BrownCompleteness):
BrownCompleteness object
s (float or ~numpy.ndarray, optional):
Projected separations (in AU) to compute minimum delta mag at.
If None (default) then uses Completeness.xnew
Returns:
~numpy.ndarray:
Maximum deltaMag values
"""
Ppop = Completeness.PlanetPopulation
dmagmax = -2.5 * np.log10(
Ppop.prange.min()
* ((Ppop.Rprange.min() / Ppop.rrange.max()).decompose().value) ** 2
* Completeness.PlanetPhysicalModel.calc_Phi(
(np.pi - np.arcsin(s / Ppop.rrange.max().value)) * u.rad
)
)
return dmagmax