from typing import Optional, Union
import astropy.units as u
import numpy as np
from astropy.units import Quantity
from tqdm.auto import tqdm
from pytmosph3r.util import pt
from pytmosph3r.util.util import merge_attrs, to_SI
from .emission import Emission
[docs]
class Phasecurve(Emission):
"""The Phasecurve module relies on the :class:`~pytmosph3r.emission.Emission` module: we iterate over
:attr:`n_phases` observer longitudes and scale the flux using the projection of the surface onto that direction.
You can use the same parameters as :class:`~pytmosph3r.emission.Emission` for this module.
IMPORTANT NOTE: :attr:`top_flux_from_star` and :attr:`mu_from_obs` will add a considerable computing time since
the emission flux has to be re-computed at EACH phase.
Args:
n_phases (int) : Number of phases for the phase curve. Defaults to 100.
start (float) : Phase at which to start the phasecurve, in SI (radians) or astropy units.
end (float) : Phase at which to end the phasecurve, in SI (radians) or astropy units.
phases (array) : List of phases (ignores :attr:`n_phases`, :attr:`start` and :attr:`end`). Assumed to be in
SI, or astropy units.
kwargs (dict): Parameters passed to the :class:`~pytmosph3r.emission.Emission` module.
Returns:
(array): A series of :attr:`n_phases` Emission fluxes.
"""
def __init__(self,
n_phases: Optional[int] = None,
start: Union[None, float, Quantity[pt.angle]] = None,
end: Union[None, float, Quantity[pt.angle]] = None,
phases: Union[None, np.ndarray, Quantity[pt.angle]] = None,
wns=None, *args, **kwargs):
"""
We provide two methods to select the phases, the first one (1) provide an array with the phases wanted.
The second one (2) create the phases array from a number of phase and the bounds.
Args:
n_phases (Optional[int]): (2)
start (None, float, Quantity[pt.angle]): (2)
end (None, float, Quantity[pt.angle]): (2)
phases (Union[None, np.ndarray, Quantity[pt.angle]]): Typed array (1)
wns ():
*args ():
**kwargs ():
"""
super().__init__(*args, **kwargs)
self.n_phases: int = n_phases
"""Number of phases in the phase curve."""
self.start : float = to_SI(start, u.rad)
self.end: float = to_SI(end, u.rad)
self.phases: np.ndarray = to_SI(phases, u.rad)
self.wns: np.ndarray = wns
self.store_raw_flux: bool = True # INFO: always activates this in phasecurves
[docs]
def build(self, model):
self.model = model
# Get emission parameters...
e = model.emission
if e is not None: # ... from emission object if computed
self = merge_attrs(self, e)
else: # ... by recomputing them otherwise
super().build(model)
# Set default parameters
if self.phases is not None:
if isinstance(self.phases, (int, float)):
self.phases=[self.phases] # handle one phase case
assert hasattr(self.phases, "__len__"), "Phases should be given as a list or array."
self.warning("You have set a list of 'phases', so 'n_phases' is ignored.")
self.n_phases = len(self.phases)
if self.n_phases is None:
self.n_phases = 100
self.warning("Number of phases set to 100.")
self.n_phases = int(self.n_phases)
if self.start is None:
self.start = to_SI(-180*u.deg, u.rad)
if self.end is None:
self.end = to_SI(180*u.deg, u.rad)
if self.phases is None:
self.phases = to_SI(np.linspace(self.start, self.end, self.n_phases)*u.rad, u.rad)
[docs]
def compute(self, model):
"""Computation of a phase curve, along :attr:`n_phases` observer longitudes. This function is called by :func:`compute`."""
self.info(f"Computing phase curve with {self.n_phases} points...")
self.flux = np.zeros((self.n_phases,self.model.opacity.k_data.Nw))
initial_longitude = self.model.observer.longitude # save longitude since we're going to overwrite it during the phase curve
if self.wns is None:
self.wns = model.opacity.k_data.wns
# TODO pick wns for calculations, as in lightcurve module
# re-calculate raw flux if possible
# Check the case where we only need to perform one call to `compute`.
# - Planet is tidally locked
# and - mu_from_obs is False
# Then the position of the observer does not affect the `raw_flux` computed.
if self.model.orbit.is_tidal_locked and not self.mu_from_obs:
try:
self.raw_flux = model.emission.raw_flux
# self.flux = model.emission.flux
self.factor = model.emission.factor
self.factor_normalized = model.emission.factor_normalized
except:
super().compute(model) # compute self.flux
for i, phase in enumerate(tqdm(self.phases)):
self.model.observer.position_from_phase(phase)
if not self.model.orbit.is_tidal_locked or self.mu_from_obs:
# `raw_flux` need to be computed for each phase.
super().compute(model) # re-compute self.flux # very costly
else:
# Only need to update the projection
self.compute_projection()
self.flux[i] = np.sum(self.raw_flux * self.projected_surface[..., None], axis=(0,1))
self.model.observer.longitude = initial_longitude # resurrect initial setup
if self.planet_to_star_flux_ratio is True:
self.flux_normalized = np.copy(self.flux)
try:
self.flux_normalized /= self.factor_normalized.value
except:
self.flux_normalized /= self.factor_normalized
try:
self.flux /= self.factor.value
except:
self.flux /= self.factor # float case
[docs]
def outputs(self):
return ['wns', 'flux', 'flux_normalized']