"""
Phase Functions
See also [Keithly2021]_
"""
import numpy as np
from astropy import units as u
[docs]
def phi_lambert(beta, phiIndex=np.asarray([])):
"""Lambert phase function (isotropic scattering)
See: [Sobolev1975]_
Args:
beta (astropy.units.quantity.Quantity or numpy.ndarray):
phase angle array in radians
phiIndex (numpy.ndarray):
array of indicies of type of exoplanet phase function to use, ints 0-7
Returns:
numpy.ndarray:
Phi, phase function values between 0 and 1
"""
if isinstance(beta, u.Quantity):
beta = beta.to_value(u.rad)
phi = (np.sin(beta) + (np.pi - beta) * np.cos(beta)) / np.pi
return phi
[docs]
def transitionStart(x, a, b):
"""Smoothly transition from one 0 to 1
Args:
x (numpy.ndarray):
x, in deg input value in deg, floats
a (numpy.ndarray):
transition midpoint in deg, floats
b (numpy.ndarray):
transition slope
Returns:
numpy.ndarray:
s, Transition value from 0 to 1, floats
"""
s = 0.5 + 0.5 * np.tanh((x - a) / b)
return s
[docs]
def transitionEnd(x, a, b):
"""Smoothly transition from one 1 to 0
Smaller b is sharper step a is midpoint, s(a)=0.5
Args:
x (numpy.ndarray):
x, in deg input value in deg, floats
a (numpy.ndarray):
a, transition midpoint in deg, floats
b (numpy.ndarray):
transition slope
Returns:
numpy.ndarray:
s, transition value from 1 to 0
"""
s = 0.5 - 0.5 * np.tanh((x - a) / b)
return s
[docs]
def quasiLambertPhaseFunction(beta, phiIndex=np.asarray([])):
"""Quasi Lambert Phase Function from [Agol2007]_
Args:
beta (astropy.units.quantity.Quantity or numpy.ndarray):
planet phase angles in radians
phiIndex (numpy.ndarray):
array of indicies of type of exoplanet phase function to use, ints 0-7
Returns:
numpy.ndarray:
Phi, phase function value
"""
if isinstance(beta, u.Quantity):
beta = beta.to_value(u.rad)
Phi = np.cos(beta / 2.0) ** 4
return Phi
[docs]
def quasiLambertPhaseFunctionInverse(Phi, phiIndex=np.asarray([])):
"""Quasi Lambert Phase Function Inverse
Args:
Phi (numpy.ndarray):
Phi, phase function value, floats
phiIndex (numpy.ndarray):
array of indicies of type of exoplanet phase function to use, ints 0-7
Returns:
numpy.ndarray:
beta, planet phase angles in rad, floats
"""
beta = 2.0 * np.arccos((Phi) ** (1.0 / 4.0))
return beta
[docs]
def hyperbolicTangentPhaseFunc(beta, A, B, C, D, planetName=None):
"""
Optimal Parameters for Earth Phase Function basedon mallama2018 comparison using
mallama2018PlanetProperties.py:
A=1.85908529, B=0.89598952, C=1.04850586, D=-0.08084817
Optimal Parameters for All Solar System Phase Function basedon mallama2018
comparison using mallama2018PlanetProperties.py:
A=0.78415 , B=1.86890455, C=0.5295894 , D=1.07587213
Args:
beta (astropy.units.quantity.Quantity or numpy.ndarray):
Phase Angle in radians
A (float):
Hyperbolic phase function parameter
B (float):
Hyperbolic phase function paramter
C (float):
Hyperbolic phase function parameter
D (float):
Hyperbolic phase function parameter
planetName (string or None):
planet name string all lower case for one of 8 solar system planets
Returns:
numpy.ndarray:
Phi, phase angle in degrees
"""
if planetName is None:
return None # do nothing
elif planetName == "mercury":
A, B, C, D = (
0.9441564,
0.3852919,
2.59291159,
-0.67540991,
) # 0.93940195, 0.40446512, 2.47034733, -0.64361749
elif planetName == "venus":
A, B, C, D = (
1.20324931,
1.57402581,
0.60683886,
0.87010846,
) # 1.26116739, 1.53204409, 0.61961161, 0.84075693
elif planetName == "earth":
A, B, C, D = (
0.78414986,
1.86890464,
0.52958938,
1.07587225,
) # 0.78415 , 1.86890455, 0.5295894 , 1.07587213
elif planetName == "mars":
A, B, C, D = (
1.89881785,
0.48220465,
2.02299497,
-1.02612681,
) # 2.02856459, 0.29590061, 3.32324214, -1.71048535
elif planetName == "jupiter":
A, B, C, D = (
1.3622441,
1.45676529,
0.64490468,
0.7800663,
) # 1.3761512 , 1.45852349, 0.64157352, 0.7983722
elif planetName == "saturn":
A, B, C, D = (
5.02672862,
0.41588155,
1.80205383,
-1.74369974,
) # 5.49410541, 0.37274869, 2.00119662, -2.1551928
elif planetName == "uranus":
A, B, C, D = (
1.54388146,
1.18304642,
0.79972526,
0.37288376,
) # 1.56866334, 1.16284633, 0.81250327, 0.34759469
elif planetName == "neptune":
A, B, C, D = (
1.31369238,
1.41437107,
0.67584636,
0.65077278,
) # 1.37105297, 1.36886173, 0.69506274, 0.609515
if hasattr(beta, "value"):
beta = beta.to("rad").value
Phi = -np.tanh((beta - D) / A) / B + C
return Phi
[docs]
def hyperbolicTangentPhaseFuncInverse(Phi, A, B, C, D, planetName=None):
"""
Optimal Parameters for Earth Phase Function basedon mallama2018 comparison
using mallama2018PlanetProperties.py:
A=1.85908529, B=0.89598952, C=1.04850586, D=-0.08084817
Optimal Parameters for All Solar System Phase Function basedon mallama2018
comparison using mallama2018PlanetProperties.py:
A=0.78415 , B=1.86890455, C=0.5295894 , D=1.07587213
Args:
Phi (numpy.ndarray):
phase angle in degrees
A (float):
Hyperbolic phase function parameter
B (float):
Hyperbolic phase function paramter
C (float):
Hyperbolic phase function parameter
D (float):
Hyperbolic phase function parameter
planetName (string or None):
planet name string all lower case for one of 8 solar system planets
Returns:
numpy.ndarray:
beta, Phase Angle in degrees
"""
if planetName is None:
return None # do nothing
elif planetName == "mercury":
A, B, C, D = (
0.9441564,
0.3852919,
2.59291159,
-0.67540991,
) # 0.93940195, 0.40446512, 2.47034733, -0.64361749
elif planetName == "venus":
A, B, C, D = (
1.20324931,
1.57402581,
0.60683886,
0.87010846,
) # 1.26116739, 1.53204409, 0.61961161, 0.84075693
elif planetName == "earth":
A, B, C, D = (
0.78414986,
1.86890464,
0.52958938,
1.07587225,
) # 0.78415 , 1.86890455, 0.5295894 , 1.07587213
elif planetName == "mars":
A, B, C, D = (
1.89881785,
0.48220465,
2.02299497,
-1.02612681,
) # 2.02856459, 0.29590061, 3.32324214, -1.71048535
elif planetName == "jupiter":
A, B, C, D = (
1.3622441,
1.45676529,
0.64490468,
0.7800663,
) # 1.3761512 , 1.45852349, 0.64157352, 0.7983722
elif planetName == "saturn":
A, B, C, D = (
5.02672862,
0.41588155,
1.80205383,
-1.74369974,
) # 5.49410541, 0.37274869, 2.00119662, -2.1551928
elif planetName == "uranus":
A, B, C, D = (
1.54388146,
1.18304642,
0.79972526,
0.37288376,
) # 1.56866334, 1.16284633, 0.81250327, 0.34759469
elif planetName == "neptune":
A, B, C, D = (
1.31369238,
1.41437107,
0.67584636,
0.65077278,
) # 1.37105297, 1.36886173, 0.69506274, 0.609515
beta = ((A * np.arctanh(-B * (Phi - C)) + D) * u.radian).to("deg").value
return beta
[docs]
def betaFunc(inc, v, w):
"""Calculate the planet phase angle
Args:
inc (float or numpy.ndarray):
planet inclination in rad
v (numpy.ndarray):
planet true anomaly in rad
w (numpy.ndarray):
planet argument of periapsis
Returns:
numpy.ndarray:
beta, planet phase angle
"""
beta = np.arccos(np.sin(inc) * np.sin(v + w))
return beta
[docs]
def phase_Mercury(beta):
"""Mercury phase function
Valid from 0 to 180 deg
Args:
beta (numpy.ndarray):
beta, phase angle in degrees
Returns:
numpy.ndarray:
phase function values
"""
phase = 10.0 ** (
-0.4
* (
6.3280e-02 * beta
- 1.6336e-03 * beta**2.0
+ 3.3644e-05 * beta**3.0
- 3.4265e-07 * beta**4.0
+ 1.6893e-09 * beta**5.0
- 3.0334e-12 * beta**6.0
)
)
return phase
[docs]
def phase_Venus_1(beta):
"""Venus phase function
Valid from 0 to 163.7 deg
Args:
beta (numpy.ndarray):
beta, phase angle in degrees
Returns:
numpy.ndarray:
phase function values
"""
phase = 10.0 ** (
-0.4
* (
-1.044e-03 * beta
+ 3.687e-04 * beta**2.0
- 2.814e-06 * beta**3.0
+ 8.938e-09 * beta**4.0
)
)
return phase
[docs]
def phase_Venus_2(beta):
"""Venus phase function
Valid from 163.7 to 179 deg
Args:
beta (numpy.ndarray):
beta, phase angle in degrees
Returns:
numpy.ndarray:
phase function values
"""
phase = 10.0 ** (-0.4 * (-2.81914e-00 * beta + 8.39034e-03 * beta**2.0))
# 1 Scale Properly
h1 = phase_Venus_1(163.7) - 0.0 # Total height desired over range
h2 = 10.0 ** (
-0.4 * (-2.81914e-00 * 163.7 + 8.39034e-03 * 163.7**2.0)
) - 10.0 ** (-0.4 * (-2.81914e-00 * 179.0 + 8.39034e-03 * 179.0**2.0))
phase = phase * h1 / h2 # Scale so height is proper
# 2 Lateral movement to make two functions line up
difference = phase_Venus_1(163.7) - h1 / h2 * (
10.0 ** (-0.4 * (-2.81914e-00 * 163.7 + 8.39034e-03 * 163.7**2.0))
)
phase = phase + difference
return phase
[docs]
def phase_Venus_melded(beta):
"""
Venus phae function
Valid from 0 to 180 deg
Args:
beta (numpy.ndarray):
beta, phase angle in degrees
Returns:
numpy.ndarray:
phase function values
"""
phase = (
transitionEnd(beta, 163.7, 5.0) * phase_Venus_1(beta)
+ transitionStart(beta, 163.7, 5.0)
* transitionEnd(beta, 179.0, 0.5)
* phase_Venus_2(beta)
+ transitionStart(beta, 179.0, 0.5) * phi_lambert(beta * np.pi / 180.0)
+ 2.766e-04
)
# 2.666e-04 ensures the phase function is entirely positive
# (near 180 deg phase, there is a small region
# where phase goes negative) This small addition fixes this
return phase
[docs]
def phase_Earth(beta):
"""Earth phase function
Valid from 0 to 180 deg
Args:
beta (numpy.ndarray):
beta, phase angle in degrees
Returns:
numpy.ndarray:
phase function values
"""
phase = 10.0 ** (-0.4 * (-1.060e-3 * beta + 2.054e-4 * beta**2.0))
return phase
[docs]
def phase_Mars_1(beta):
"""Mars phase function
Valid from 0 to 50 deg
Args:
beta (numpy.ndarray):
beta, phase angle in degrees
Returns:
numpy.ndarray:
phase function values
"""
phase = 10.0 ** (
-0.4 * (0.02267 * beta - 0.0001302 * beta**2.0 + 0.0 + 0.0)
) # L(λe) + L(LS)
return phase
[docs]
def phase_Mars_2(beta):
"""Mars phase function
Valid from 50 to 180 deg
Args:
beta (numpy.ndarray):
beta, phase angle in degrees
Returns:
numpy.ndarray:
phase function values
"""
phase = (
phase_Mars_1(50.0)
/ 10.0 ** (-0.4 * (-0.02573 * 50.0 + 0.0003445 * 50.0**2.0))
* 10.0 ** (-0.4 * (-0.02573 * beta + 0.0003445 * beta**2.0 + 0.0 + 0.0))
) # L(λe) + L(Ls)
return phase
[docs]
def phase_Mars_melded(beta):
"""Mars phase function
Valid from 0 to 180 degrees
Args:
beta (numpy.ndarray):
beta, phase angle in degrees
Returns:
numpy.ndarray:
phase function values
"""
phase = transitionEnd(beta, 50.0, 5.0) * phase_Mars_1(beta) + transitionStart(
beta, 50.0, 5.0
) * phase_Mars_2(beta)
return phase
[docs]
def phase_Jupiter_1(beta):
"""Jupiter phase function
Valid from 0 to 12 deg
Args:
beta (numpy.ndarray):
beta, phase angle in degrees
Returns:
numpy.ndarray:
phase function values
"""
phase = 10.0 ** (-0.4 * (-3.7e-04 * beta + 6.16e-04 * beta**2.0))
return phase
[docs]
def phase_Jupiter_2(beta):
"""Jupiter phase function
Valid from 12 to 130 deg
Args:
beta (numpy.ndarray):
beta, phase angle in degrees
Returns:
numpy.ndarray:
phase function values
"""
# inds = np.where(beta > 180.0)[0]
# beta[inds] = [180.0] * len(inds)
# assert np.all(
# (
# 1.0
# - 1.507 * (beta / 180.0)
# - 0.363 * (beta / 180.0) ** 2.0
# - 0.062 * (beta / 180.0) ** 3.0
# + 2.809 * (beta / 180.0) ** 4.0
# - 1.876 * (beta / 180.0) ** 5.0
# )
# >= 0.0
# ), "error in beta input"
difference = phase_Jupiter_1(12.0) - 10.0 ** (
-0.4
* (
-2.5
* np.log10(
1.0
- 1.507 * (12.0 / 180.0)
- 0.363 * (12.0 / 180.0) ** 2.0
- 0.062 * (12.0 / 180.0) ** 3.0
+ 2.809 * (12.0 / 180.0) ** 4.0
- 1.876 * (12.0 / 180.0) ** 5.0
)
)
)
phase = difference + 10.0 ** (
-0.4
* (
-2.5
* np.log10(
1.0
- 1.507 * (beta / 180.0)
- 0.363 * (beta / 180.0) ** 2.0
- 0.062 * (beta / 180.0) ** 3.0
+ 2.809 * (beta / 180.0) ** 4.0
- 1.876 * (beta / 180.0) ** 5.0
)
)
)
return phase
[docs]
def phase_Jupiter_melded(beta):
"""Jupiter phase function
Valid from 0 to 130 degrees
Args:
beta (numpy.ndarray):
beta, phase angle in degrees
Returns:
numpy.ndarray:
phase function values
"""
phase = (
transitionEnd(beta, 12.0, 5.0) * phase_Jupiter_1(beta)
+ transitionStart(beta, 12.0, 5.0)
* transitionEnd(beta, 130.0, 5.0)
* phase_Jupiter_2(beta)
+ transitionStart(beta, 130.0, 5.0) * phi_lambert(beta * np.pi / 180.0)
)
return phase
[docs]
def phase_Saturn_2(beta):
"""Saturn phase function (Globe Only Earth Observations)
Valid beta from 0 to 6.5 deg
Args:
beta (numpy.ndarray):
beta, phase angle in degrees
Returns:
numpy.ndarray:
phase function values
"""
phase = 10.0 ** (-0.4 * (-3.7e-04 * beta + 6.16e-04 * beta**2.0))
return phase
[docs]
def phase_Saturn_3(beta):
"""Saturn phase function (Globe Only Pioneer Observations)
Valid beta from 6 to 150. deg
Args:
beta (numpy.ndarray):
beta, phase angle in degrees
Returns:
numpy.ndarray:
phase function values
"""
difference = phase_Saturn_2(6.5) - 10.0 ** (
-0.4
* (
2.446e-4 * 6.5
+ 2.672e-4 * 6.5**2.0
- 1.505e-6 * 6.5**3.0
+ 4.767e-9 * 6.5**4.0
)
)
phase = difference + 10.0 ** (
-0.4
* (
2.446e-4 * beta
+ 2.672e-4 * beta**2.0
- 1.505e-6 * beta**3.0
+ 4.767e-9 * beta**4.0
)
)
return phase
[docs]
def phase_Saturn_melded(beta):
"""
Saturn phase function
Args:
beta (numpy.ndarray):
beta, phase angle in degrees
Returns:
numpy.ndarray:
phase function values
"""
phase = (
transitionEnd(beta, 6.5, 5.0) * phase_Saturn_2(beta)
+ transitionStart(beta, 6.5, 5.0)
* transitionEnd(beta, 150.0, 5.0)
* phase_Saturn_3(beta)
+ transitionStart(beta, 150.0, 5.0) * phi_lambert(beta * np.pi / 180.0)
)
return phase
[docs]
def phiprime_phi(phi):
"""Helper method for Uranus phase function
Valid for phi from -82 to 82 deg
Args:
phi (numpy.ndarray):
phi, planet rotation axis offset in degrees, floats
Returns:
numpy.ndarray:
phiprime, in deg, floats
"""
f = 0.0022927 # flattening of the planet
phiprime = np.arctan2(np.tan(phi * np.pi / 180.0), (1.0 - f) ** 2.0) * 180.0 / np.pi
return phiprime
[docs]
def phase_Uranus(beta, phi=-82.0):
"""Uranus phase function
Valid for beta 0 to 154 deg
Args:
beta (numpy.ndarray):
beta, phase angle in degrees
phi (float or numpy.ndarray):
phi, planet rotation axis offset in degrees, floats
Returns:
numpy.ndarray:
phase function values
"""
phase = 10.0 ** (
-0.4 * (-8.4e-04 * phiprime_phi(phi) + 6.587e-3 * beta + 1.045e-4 * beta**2.0)
)
return phase
[docs]
def phase_Uranus_melded(beta):
"""Uranus phase function
Args:
beta (numpy.ndarray):
beta, phase angle in degrees
Returns:
numpy.ndarray:
phase function values
"""
phase = transitionEnd(beta, 154.0, 5.0) * phase_Uranus(beta) + transitionStart(
beta, 154.0, 5.0
) * phi_lambert(beta * np.pi / 180.0)
return phase
[docs]
def phase_Neptune(beta):
"""Neptune phase function
Valid for beta 0 to 133.14 deg
Args:
beta (numpy.ndarray):
beta, phase angle in degrees
Returns:
numpy.ndarray:
phase function values
"""
phase = 10.0 ** (-0.4 * (7.944e-3 * beta + 9.617e-5 * beta**2.0))
return phase
[docs]
def phase_Neptune_melded(beta):
"""Neptune phase function
Args:
beta (numpy.ndarray):
beta, phase angle in degrees
Returns:
numpy.ndarray:
phase function values
"""
phase = transitionEnd(beta, 133.14, 5.0) * phase_Neptune(beta) + transitionStart(
beta, 133.14, 5.0
) * phi_lambert(beta * np.pi / 180.0)
return phase
[docs]
def realSolarSystemPhaseFunc(beta, phiIndex=np.asarray([])):
"""Uses the phase functions from Mallama 2018 implemented in
mallama2018PlanetProperties.py
Args:
beta (astropy.units.quantity.Quantity or numpy.ndarray):
phase angle array in degrees
phiIndex (numpy.ndarray):
array of indicies of type of exoplanet phase function to use, ints 0-7
Returns:
numpy.ndarray:
Phi, phase function values between 0 and 1
"""
if len(phiIndex) == 0: # Default behavior is to use the lambert phase function
Phi = np.zeros(len(beta)) # instantiate initial array
Phi = phi_lambert(beta)
else:
if hasattr(beta, "unit"):
beta = beta.to("rad").value
beta = beta * 180.0 / np.pi # convert to phase angle in degrees
if not len(phiIndex) == 0 and (len(beta) == 1 or len(beta) == 0):
beta = np.ones(len(phiIndex)) * beta
Phi = np.zeros(len(beta))
# Find indicies of where to use each phase function
mercuryInds = np.where(phiIndex == 0)[0]
venusInds = np.where(phiIndex == 1)[0]
earthInds = np.where(phiIndex == 2)[0]
marsInds = np.where(phiIndex == 3)[0]
jupiterInds = np.where(phiIndex == 4)[0]
saturnInds = np.where(phiIndex == 5)[0]
uranusInds = np.where(phiIndex == 6)[0]
neptuneInds = np.where(phiIndex == 7)[0]
if not len(mercuryInds) == 0:
Phi[mercuryInds] = phase_Mercury(beta[mercuryInds])
if not len(venusInds) == 0:
Phi[venusInds] = phase_Venus_melded(beta[venusInds])
if not len(earthInds) == 0:
Phi[earthInds] = phase_Earth(beta[earthInds])
if not len(marsInds) == 0:
Phi[marsInds] = phase_Mars_melded(beta[marsInds])
if not len(jupiterInds) == 0:
Phi[jupiterInds] = phase_Jupiter_melded(beta[jupiterInds])
if not len(saturnInds) == 0:
Phi[saturnInds] = phase_Saturn_melded(beta[saturnInds])
if not len(uranusInds) == 0:
Phi[uranusInds] = phase_Uranus_melded(beta[uranusInds])
if not len(neptuneInds) == 0:
Phi[neptuneInds] = phase_Neptune_melded(beta[neptuneInds])
return Phi