from typing import Tuple
import numpy as np
from solrat.engine.functions.decorators import log_method
[docs]
class RadiativeTransferCoefficients:
"""
This is a container class for all radiative transfer coefficients.
"""
def __init__(
self,
eta_rho_aI,
eta_rho_aQ,
eta_rho_aU,
eta_rho_aV,
eta_rho_sI,
eta_rho_sQ,
eta_rho_sU,
eta_rho_sV,
epsilonI,
epsilonQ,
epsilonU,
epsilonV,
):
self.eta_rho_aI = eta_rho_aI
self.eta_rho_aQ = eta_rho_aQ
self.eta_rho_aU = eta_rho_aU
self.eta_rho_aV = eta_rho_aV
self.eta_rho_sI = eta_rho_sI
self.eta_rho_sQ = eta_rho_sQ
self.eta_rho_sU = eta_rho_sU
self.eta_rho_sV = eta_rho_sV
self.epsilonI = epsilonI
self.epsilonQ = epsilonQ
self.epsilonU = epsilonU
self.epsilonV = epsilonV
self._eta_tau_scale = self._compute_eta_tau_scale()
[docs]
def get_eta_I(self):
return np.real(self.eta_rho_aI - self.eta_rho_sI)
[docs]
def get_eta_Q(self):
return np.real(self.eta_rho_aQ - self.eta_rho_sQ)
[docs]
def get_eta_U(self):
return np.real(self.eta_rho_aU - self.eta_rho_sU)
[docs]
def get_eta_V(self):
return np.real(self.eta_rho_aV - self.eta_rho_sV)
[docs]
def get_rho_I(self):
return np.imag(self.eta_rho_aI - self.eta_rho_sI)
[docs]
def get_rho_Q(self):
return np.imag(self.eta_rho_aQ - self.eta_rho_sQ)
[docs]
def get_rho_U(self):
return np.imag(self.eta_rho_aU - self.eta_rho_sU)
[docs]
def get_rho_V(self):
return np.imag(self.eta_rho_aV - self.eta_rho_sV)
[docs]
def get_epsilon_I(self):
return np.real(self.epsilonI)
[docs]
def get_epsilon_Q(self):
return np.real(self.epsilonQ)
[docs]
def get_epsilon_U(self):
return np.real(self.epsilonU)
[docs]
def get_epsilon_V(self):
return np.real(self.epsilonV)
[docs]
@log_method
def split_eta_rho(
self,
) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
"""
Split real and imaginary parts of eta_rho for constructing the propagation matrix K
"""
etaI = np.real(self.eta_rho_aI - self.eta_rho_sI)
etaQ = np.real(self.eta_rho_aQ - self.eta_rho_sQ)
etaU = np.real(self.eta_rho_aU - self.eta_rho_sU)
etaV = np.real(self.eta_rho_aV - self.eta_rho_sV)
rhoQ = np.imag(self.eta_rho_aQ - self.eta_rho_sQ)
rhoU = np.imag(self.eta_rho_aU - self.eta_rho_sU)
rhoV = np.imag(self.eta_rho_aV - self.eta_rho_sV)
return etaI, etaQ, etaU, etaV, rhoQ, rhoU, rhoV
[docs]
@log_method
def K_z(self) -> np.ndarray:
r"""
RT matrix K related to the z-propagation equation:
.. math::
dStokes/dz = - K * Stokes + epsilon - K_continuum Stokes + epsilon_continuum
K = K_A - K_S
If N = 1 was used (no atom concentration provided), then the RTE code computes primed Kʹ and epsilonʹ:
.. math::
Kʹ = K / N
epsilonʹ = epsilon / N
dStokes/dz = - Kʹ * N * Stokes + epsilonʹ * N [ - K_continuum Stokes + epsilon_continuum ]
But N value does not impact the K_tau kernel (see below)
Reference: (6.83-6.85)
"""
etaI, etaQ, etaU, etaV, rhoQ, rhoU, rhoV = self.split_eta_rho()
# shape: [len(nu), 4, 4]
K = np.empty((len(etaI), 4, 4), dtype=np.float64)
K[:, 0, 0] = etaI
K[:, 0, 1] = etaQ
K[:, 0, 2] = etaU
K[:, 0, 3] = etaV
K[:, 1, 0] = etaQ
K[:, 1, 1] = etaI
K[:, 1, 2] = rhoV
K[:, 1, 3] = -rhoU
K[:, 2, 0] = etaU
K[:, 2, 1] = -rhoV
K[:, 2, 2] = etaI
K[:, 2, 3] = rhoQ
K[:, 3, 0] = etaV
K[:, 3, 1] = rhoU
K[:, 3, 2] = -rhoQ
K[:, 3, 3] = etaI
return K
[docs]
@log_method
def K_tau(self) -> np.ndarray:
r"""
RT matrix K related to the tau-propagation equation:
.. math::
dtau_line = - eta_line * dz
dStokes/dtau_line = K_tau_line * Stokes - epsilon_tau_line + Stokes / eta_LC - BP(T) eI / eta_LC
eta_LC = eta_line / eta_continuum = dtau_line / dtau_continuum
Stokes[tau->+inf] -> BP(T0)
Reference: modified (9.36), (9.42)
"""
K = self.K_z()
return K / self._eta_tau_scale
def _compute_eta_tau_scale(self) -> np.ndarray:
"""
Line-center eta, assuming a single spectral line
"""
scale = np.max(np.abs(np.real(self.eta_rho_aI - self.eta_rho_sI)))
if scale == 0:
raise ValueError(
"eta_tau_scale is zero: the spectral line has no opacity over the supplied frequency grid. "
"Check that the frequency grid is centered on the actual transition wavelength."
)
return scale
[docs]
@log_method
def epsilon_z(self) -> np.ndarray:
"""
RT matrix emission coefficient epsilon related to the z-propagation equation: see the comments for K_z()
Reference: (6.83-6.85)
"""
# shape: [len(nu), 4, 1]
epsilon = np.zeros((len(self.eta_rho_sV), 4, 1), dtype=np.float64)
epsilon[:, 0, 0] = self.epsilonI
epsilon[:, 1, 0] = self.epsilonQ
epsilon[:, 2, 0] = self.epsilonU
epsilon[:, 3, 0] = self.epsilonV
return epsilon
[docs]
@log_method
def epsilon_tau(self) -> np.ndarray:
"""
RT matrix emission coefficient epsilon related to the tau-propagation equation: see the comments for K_tau()
Reference: modified (9.36)
"""
return self.epsilon_z() / self._eta_tau_scale