Photochemistry

Module related to loading in photodisocciation data.

CrossSection

Loads cross-section data for molecule.

Source code in src/freckll/reactions/photo.py

class CrossSection:
    """Loads cross-section data for molecule."""

    def __init__(self, molecule: SpeciesFormula, wavelength: u.Quantity, cross_section: u.Quantity) -> None:
        """Initialize and load the cross-section data.

        Args:
            molecule: The molecule for which the cross-section is loaded.
            wavelength: The wavelength at which the cross-section is measured.
            cross_section: The cross-section value at the given wavelength.

        """
        self.wavelength = wavelength
        self.cross_section = cross_section
        self.molecule = molecule

    def interp_to(self, wavelength: u.Quantity, temperature: u.Quantity, pressure: u.Quantity) -> "CrossSection":
        """Interpolate the cross section to the given wavelength, temperature and pressure

        Args:
            wavelength: The wavelength to which the cross section is interpolated.
            temperature: The temperature at which the cross section is measured.
            pressure: The pressure at which the cross section is measured.

        Returns:
            CrossSection: A new CrossSection object with the interpolated values.


        """

        new_wavelength = self.wavelength.to(wavelength.unit, equivalencies=u.spectral())

        return CrossSection(
            self.molecule,
            wavelength,
            np.interp(wavelength.value, new_wavelength.value, self.cross_section.value, left=0.0, right=0.0)
            << self.cross_section.unit,
        )

    def __add__(self, other: "CrossSection") -> "CrossSection":
        """Add two cross sections together.

        Args:
            other: The other cross section to add.

        Returns:
            CrossSection: A new CrossSection object with the summed values.

        """
        if not isinstance(other, CrossSection):
            raise TypeError("Can only add CrossSection objects")

        if self.molecule != other.molecule:
            raise ValueError("Cannot add cross sections of different molecules")

        return CrossSection(self.molecule, self.wavelength, self.cross_section + other.cross_section)

__add__(other)

Add two cross sections together.

Parameters:

Name	Type	Description	Default
other	CrossSection	The other cross section to add.	required

Returns:

Name	Type	Description
CrossSection	CrossSection	A new CrossSection object with the summed values.

Source code in src/freckll/reactions/photo.py

def __add__(self, other: "CrossSection") -> "CrossSection":
    """Add two cross sections together.

    Args:
        other: The other cross section to add.

    Returns:
        CrossSection: A new CrossSection object with the summed values.

    """
    if not isinstance(other, CrossSection):
        raise TypeError("Can only add CrossSection objects")

    if self.molecule != other.molecule:
        raise ValueError("Cannot add cross sections of different molecules")

    return CrossSection(self.molecule, self.wavelength, self.cross_section + other.cross_section)

__init__(molecule, wavelength, cross_section)

Initialize and load the cross-section data.

Parameters:

Name	Type	Description	Default
molecule	SpeciesFormula	The molecule for which the cross-section is loaded.	required
wavelength	Quantity	The wavelength at which the cross-section is measured.	required
cross_section	Quantity	The cross-section value at the given wavelength.	required

Source code in src/freckll/reactions/photo.py

def __init__(self, molecule: SpeciesFormula, wavelength: u.Quantity, cross_section: u.Quantity) -> None:
    """Initialize and load the cross-section data.

    Args:
        molecule: The molecule for which the cross-section is loaded.
        wavelength: The wavelength at which the cross-section is measured.
        cross_section: The cross-section value at the given wavelength.

    """
    self.wavelength = wavelength
    self.cross_section = cross_section
    self.molecule = molecule

interp_to(wavelength, temperature, pressure)

Interpolate the cross section to the given wavelength, temperature and pressure

Parameters:

Name	Type	Description	Default
wavelength	Quantity	The wavelength to which the cross section is interpolated.	required
temperature	Quantity	The temperature at which the cross section is measured.	required
pressure	Quantity	The pressure at which the cross section is measured.	required

Returns:

Name	Type	Description
CrossSection	CrossSection	A new CrossSection object with the interpolated values.

Source code in src/freckll/reactions/photo.py

def interp_to(self, wavelength: u.Quantity, temperature: u.Quantity, pressure: u.Quantity) -> "CrossSection":
    """Interpolate the cross section to the given wavelength, temperature and pressure

    Args:
        wavelength: The wavelength to which the cross section is interpolated.
        temperature: The temperature at which the cross section is measured.
        pressure: The pressure at which the cross section is measured.

    Returns:
        CrossSection: A new CrossSection object with the interpolated values.


    """

    new_wavelength = self.wavelength.to(wavelength.unit, equivalencies=u.spectral())

    return CrossSection(
        self.molecule,
        wavelength,
        np.interp(wavelength.value, new_wavelength.value, self.cross_section.value, left=0.0, right=0.0)
        << self.cross_section.unit,
    )

PhotoMolecule

Represents a molecule used in photodissociation reactions.

This class contains the cross-section data and quantum yields for the molecule.

Source code in src/freckll/reactions/photo.py

class PhotoMolecule:
    """Represents a molecule used in photodissociation reactions.

    This class contains the cross-section data and quantum yields for the molecule.

    """

    def __init__(self, molecule: SpeciesFormula, cross_section: CrossSection):
        """Initialize the PhotoMolecule with a molecule and its cross-section.

        Args:
            molecule: The molecule for which the cross-section is loaded.
            cross_section: The cross-section data for the molecule.



        """
        self.molecule = molecule
        self.cross_section = cross_section
        self.quantum_yields: dict[str | int, QuantumYield] = {}

    def add_quantum_yield(self, branch_id: int | str, quantum_yield: QuantumYield):
        """Add a quantum yield to the molecule.

        Will interpolate the quantum yield to the cross-section wavelength.

        Args:
            branch_id: The ID of the branching ratio.
            quantum_yield: The quantum yield data to be added.


        """
        self.quantum_yields[branch_id] = quantum_yield.interp_to(self.cross_section.wavelength)

    def get_quantum_yield(self, branch_id: int | str) -> QuantumYield:
        """Get the quantum yield for a given branch.

        If the quantum yield is not found, a default quantum yield of 1.0 is returned.

        Args:
            branch_id: The ID of the branching ratio.
        Returns:
            QuantumYield: The quantum yield data for the given branch.


        """

        if branch_id not in self.quantum_yields:
            return QuantumYield(
                self.molecule, branch_id, self.cross_section.wavelength, np.ones(self.cross_section.wavelength.shape)
            )

        return self.quantum_yields[branch_id]

    def interp_to(self, wavelength: u.Quantity, temperature: u.Quantity, pressure: u.Quantity) -> "PhotoMolecule":
        """Interpolate the cross section and quantum yields to the given wavelength.

        Args:
            wavelength: The wavelength to which the cross section and quantum yields are interpolated.

        Returns:
            PhotoMolecule: A new PhotoMolecule object with the interpolated values.


        """

        new_cross_section = self.cross_section.interp_to(wavelength, temperature, pressure)
        new_molecule = PhotoMolecule(self.molecule, new_cross_section)

        for branch_id, qy in self.quantum_yields.items():
            new_molecule.add_quantum_yield(branch_id, qy.interp_to(wavelength))

        return new_molecule

    def reaction_rate(self, branch_id: int | str, flux: u.Quantity) -> list[FreckllArray]:
        r"""Compute the reaction rate for a given branch.

        Calculates the integral:

        $$
        R = \int \sigma(\lambda) \cdot QY(\lambda) \cdot F(\lambda) d\lambda
        $$

        where:
        - $R$ is the reaction rate
        - $\sigma(\lambda)$ is the cross-section
        - $QY(\lambda)$ is the quantum yield
        - $F(\lambda)$ is the flux
        The integral is computed over the wavelength range of the cross-section.

        Args:
            branch_id: The ID of the branching ratio.
            flux: The flux at which the reaction rate is computed.

        Returns:
            u.Quantity: The reaction rate for the given branch.


        """
        qy = self.get_quantum_yield(branch_id)
        flux = flux.to(u.photon / u.cm**2 / u.s / u.nm, equivalencies=u.spectral_density(self.cross_section.wavelength))
        reaction_rate = (
            np.trapezoid(qy.qy * self.cross_section.cross_section * flux, self.cross_section.wavelength, axis=-1)
            / u.photon
        )
        return reaction_rate.to(1 / (u.s)).value

__init__(molecule, cross_section)

Initialize the PhotoMolecule with a molecule and its cross-section.

Parameters:

Name	Type	Description	Default
molecule	SpeciesFormula	The molecule for which the cross-section is loaded.	required
cross_section	CrossSection	The cross-section data for the molecule.	required

Source code in src/freckll/reactions/photo.py

def __init__(self, molecule: SpeciesFormula, cross_section: CrossSection):
    """Initialize the PhotoMolecule with a molecule and its cross-section.

    Args:
        molecule: The molecule for which the cross-section is loaded.
        cross_section: The cross-section data for the molecule.



    """
    self.molecule = molecule
    self.cross_section = cross_section
    self.quantum_yields: dict[str | int, QuantumYield] = {}

add_quantum_yield(branch_id, quantum_yield)

Add a quantum yield to the molecule.

Will interpolate the quantum yield to the cross-section wavelength.

Parameters:

Name	Type	Description	Default
branch_id	int | str	The ID of the branching ratio.	required
quantum_yield	QuantumYield	The quantum yield data to be added.	required

Source code in src/freckll/reactions/photo.py

def add_quantum_yield(self, branch_id: int | str, quantum_yield: QuantumYield):
    """Add a quantum yield to the molecule.

    Will interpolate the quantum yield to the cross-section wavelength.

    Args:
        branch_id: The ID of the branching ratio.
        quantum_yield: The quantum yield data to be added.


    """
    self.quantum_yields[branch_id] = quantum_yield.interp_to(self.cross_section.wavelength)

get_quantum_yield(branch_id)

Get the quantum yield for a given branch.

If the quantum yield is not found, a default quantum yield of 1.0 is returned.

Parameters:

Name	Type	Description	Default
branch_id	int | str	The ID of the branching ratio.	required

Returns: QuantumYield: The quantum yield data for the given branch.

Source code in src/freckll/reactions/photo.py

def get_quantum_yield(self, branch_id: int | str) -> QuantumYield:
    """Get the quantum yield for a given branch.

    If the quantum yield is not found, a default quantum yield of 1.0 is returned.

    Args:
        branch_id: The ID of the branching ratio.
    Returns:
        QuantumYield: The quantum yield data for the given branch.


    """

    if branch_id not in self.quantum_yields:
        return QuantumYield(
            self.molecule, branch_id, self.cross_section.wavelength, np.ones(self.cross_section.wavelength.shape)
        )

    return self.quantum_yields[branch_id]

interp_to(wavelength, temperature, pressure)

Interpolate the cross section and quantum yields to the given wavelength.

Parameters:

Name	Type	Description	Default
wavelength	Quantity	The wavelength to which the cross section and quantum yields are interpolated.	required

Returns:

Name	Type	Description
PhotoMolecule	PhotoMolecule	A new PhotoMolecule object with the interpolated values.

Source code in src/freckll/reactions/photo.py

def interp_to(self, wavelength: u.Quantity, temperature: u.Quantity, pressure: u.Quantity) -> "PhotoMolecule":
    """Interpolate the cross section and quantum yields to the given wavelength.

    Args:
        wavelength: The wavelength to which the cross section and quantum yields are interpolated.

    Returns:
        PhotoMolecule: A new PhotoMolecule object with the interpolated values.


    """

    new_cross_section = self.cross_section.interp_to(wavelength, temperature, pressure)
    new_molecule = PhotoMolecule(self.molecule, new_cross_section)

    for branch_id, qy in self.quantum_yields.items():
        new_molecule.add_quantum_yield(branch_id, qy.interp_to(wavelength))

    return new_molecule

reaction_rate(branch_id, flux)

Compute the reaction rate for a given branch.

Calculates the integral:

\[ R = \int \sigma(\lambda) \cdot QY(\lambda) \cdot F(\lambda) d\lambda \]

where: - \(R\) is the reaction rate - \(\sigma(\lambda)\) is the cross-section - \(QY(\lambda)\) is the quantum yield - \(F(\lambda)\) is the flux The integral is computed over the wavelength range of the cross-section.

Parameters:

Name	Type	Description	Default
branch_id	int | str	The ID of the branching ratio.	required
flux	Quantity	The flux at which the reaction rate is computed.	required

Returns:

Type	Description
list[FreckllArray]	u.Quantity: The reaction rate for the given branch.

Source code in src/freckll/reactions/photo.py

def reaction_rate(self, branch_id: int | str, flux: u.Quantity) -> list[FreckllArray]:
    r"""Compute the reaction rate for a given branch.

    Calculates the integral:

    $$
    R = \int \sigma(\lambda) \cdot QY(\lambda) \cdot F(\lambda) d\lambda
    $$

    where:
    - $R$ is the reaction rate
    - $\sigma(\lambda)$ is the cross-section
    - $QY(\lambda)$ is the quantum yield
    - $F(\lambda)$ is the flux
    The integral is computed over the wavelength range of the cross-section.

    Args:
        branch_id: The ID of the branching ratio.
        flux: The flux at which the reaction rate is computed.

    Returns:
        u.Quantity: The reaction rate for the given branch.


    """
    qy = self.get_quantum_yield(branch_id)
    flux = flux.to(u.photon / u.cm**2 / u.s / u.nm, equivalencies=u.spectral_density(self.cross_section.wavelength))
    reaction_rate = (
        np.trapezoid(qy.qy * self.cross_section.cross_section * flux, self.cross_section.wavelength, axis=-1)
        / u.photon
    )
    return reaction_rate.to(1 / (u.s)).value

PhotoReactionCall

Source code in src/freckll/reactions/photo.py

class PhotoReactionCall:
    def __init__(
        self,
        reactant: PhotoMolecule,
        products: list[SpeciesFormula],
        branch_id: int | str,
        species_list: t.Optional[list[SpeciesFormula]] = None,
        reactant_index: t.Optional[int] = None,
        product_indices: t.Optional[npt.NDArray[np.integer]] = None,
        tags: t.Optional[list[str]] = None,
    ):
        """Initialize the photodissociation reaction call.

        Args:
            reactant: The reactant of the reaction.
            products: The products of the reaction.
            branch_id: The ID of the branching ratio.
            species_list: The list of species in the network.
            reactant_index: The index of the reactant in the species list.
            product_indices: The indices of the products in the species list.
            tags: The tags associated with the reaction.

        """
        self.reactant = reactant
        self.products = products
        self.tags = tags if tags is not None else ["photodissociation"]

        self.reactant_index = reactant_index
        if self.reactant_index is None:
            self.reactant_index = species_list.index(reactant.molecule)
        self.product_indices = product_indices
        if self.product_indices is None:
            self.product_indices = np.array([species_list.index(p) for p in products], dtype=np.int64)
        self.branch_id = branch_id
        self.tags = list(set(self.tags))

    def interpolate_to(
        self,
        wavelength: u.Quantity,
        temperature: u.Quantity | None = None,
        pressure: u.Quantity | None = None,
    ) -> "PhotoReactionCall":
        """Interpolate the reaction call to the given wavelength.

        Args:
            wavelength: The wavelength to which the reaction call is interpolated.

        Returns:
            PhotoReactionCall: A new PhotoReactionCall object with the interpolated values.


        """
        reactant = self.reactant.interp_to(wavelength, temperature, pressure)

        return PhotoReactionCall(
            reactant,
            self.products,
            self.branch_id,
            reactant_index=self.reactant_index,
            product_indices=self.product_indices,
            tags=self.tags,
        )

    @property
    def molecule(self) -> SpeciesFormula:
        """Return the molecule of the reactant."""
        return self.reactant.molecule

    @property
    def cross_section(self) -> CrossSection:
        """Return the cross-section of the reactant."""
        return self.reactant.cross_section

    def __call__(self, flux: u.Quantity, number_density: FreckllArray) -> Reaction:
        """Call the reaction.

        This method computes the reaction rate and creates a Reaction object.

        Args:
            flux: The flux at which the reaction is computed.
            number_density: The number density of the reactant.


        Returns:
            Reaction: The reaction object representing the photodissociation reaction.

        """
        reaction_rate = self.reactant.reaction_rate(self.branch_id, flux)

        reaction = Reaction(
            reactants=[self.reactant.molecule],
            products=self.products,
            reactants_indices=np.array([self.reactant_index]),
            product_indices=self.product_indices,
            reaction_rate=reaction_rate,
            tags=self.tags,
        )

        reaction.calculate_density_krate(number_density.to(u.cm**-3).value)

        return [reaction]

cross_section property

Return the cross-section of the reactant.

molecule property

Return the molecule of the reactant.

__call__(flux, number_density)

Call the reaction.

This method computes the reaction rate and creates a Reaction object.

Parameters:

Name	Type	Description	Default
flux	Quantity	The flux at which the reaction is computed.	required
number_density	FreckllArray	The number density of the reactant.	required

Returns:

Name	Type	Description
Reaction	Reaction	The reaction object representing the photodissociation reaction.

Source code in src/freckll/reactions/photo.py

def __call__(self, flux: u.Quantity, number_density: FreckllArray) -> Reaction:
    """Call the reaction.

    This method computes the reaction rate and creates a Reaction object.

    Args:
        flux: The flux at which the reaction is computed.
        number_density: The number density of the reactant.


    Returns:
        Reaction: The reaction object representing the photodissociation reaction.

    """
    reaction_rate = self.reactant.reaction_rate(self.branch_id, flux)

    reaction = Reaction(
        reactants=[self.reactant.molecule],
        products=self.products,
        reactants_indices=np.array([self.reactant_index]),
        product_indices=self.product_indices,
        reaction_rate=reaction_rate,
        tags=self.tags,
    )

    reaction.calculate_density_krate(number_density.to(u.cm**-3).value)

    return [reaction]

__init__(reactant, products, branch_id, species_list=None, reactant_index=None, product_indices=None, tags=None)

Initialize the photodissociation reaction call.

Parameters:

Name	Type	Description	Default
reactant	PhotoMolecule	The reactant of the reaction.	required
products	list[SpeciesFormula]	The products of the reaction.	required
branch_id	int | str	The ID of the branching ratio.	required
species_list	Optional[list[SpeciesFormula]]	The list of species in the network.	None
reactant_index	Optional[int]	The index of the reactant in the species list.	None
product_indices	Optional[NDArray[integer]]	The indices of the products in the species list.	None
tags	Optional[list[str]]	The tags associated with the reaction.	None

Source code in src/freckll/reactions/photo.py

def __init__(
    self,
    reactant: PhotoMolecule,
    products: list[SpeciesFormula],
    branch_id: int | str,
    species_list: t.Optional[list[SpeciesFormula]] = None,
    reactant_index: t.Optional[int] = None,
    product_indices: t.Optional[npt.NDArray[np.integer]] = None,
    tags: t.Optional[list[str]] = None,
):
    """Initialize the photodissociation reaction call.

    Args:
        reactant: The reactant of the reaction.
        products: The products of the reaction.
        branch_id: The ID of the branching ratio.
        species_list: The list of species in the network.
        reactant_index: The index of the reactant in the species list.
        product_indices: The indices of the products in the species list.
        tags: The tags associated with the reaction.

    """
    self.reactant = reactant
    self.products = products
    self.tags = tags if tags is not None else ["photodissociation"]

    self.reactant_index = reactant_index
    if self.reactant_index is None:
        self.reactant_index = species_list.index(reactant.molecule)
    self.product_indices = product_indices
    if self.product_indices is None:
        self.product_indices = np.array([species_list.index(p) for p in products], dtype=np.int64)
    self.branch_id = branch_id
    self.tags = list(set(self.tags))

interpolate_to(wavelength, temperature=None, pressure=None)

Interpolate the reaction call to the given wavelength.

Parameters:

Name	Type	Description	Default
wavelength	Quantity	The wavelength to which the reaction call is interpolated.	required

Returns:

Name	Type	Description
PhotoReactionCall	PhotoReactionCall	A new PhotoReactionCall object with the interpolated values.

Source code in src/freckll/reactions/photo.py

def interpolate_to(
    self,
    wavelength: u.Quantity,
    temperature: u.Quantity | None = None,
    pressure: u.Quantity | None = None,
) -> "PhotoReactionCall":
    """Interpolate the reaction call to the given wavelength.

    Args:
        wavelength: The wavelength to which the reaction call is interpolated.

    Returns:
        PhotoReactionCall: A new PhotoReactionCall object with the interpolated values.


    """
    reactant = self.reactant.interp_to(wavelength, temperature, pressure)

    return PhotoReactionCall(
        reactant,
        self.products,
        self.branch_id,
        reactant_index=self.reactant_index,
        product_indices=self.product_indices,
        tags=self.tags,
    )

QuantumYield

Loads quantum yield data for molecule.

This represents the branching ratio of the photodissociation process.

This is generic and it is up to the chemical network to determine how branching ratios are organised

Source code in src/freckll/reactions/photo.py

class QuantumYield:
    """Loads quantum yield data for molecule.

    This represents the branching ratio of the photodissociation process.

    This is generic and it is up to the chemical network to determine how branching ratios
    are organised

    """

    def __init__(
        self, molecule: SpeciesFormula, branch_id: int | str, wavelength: u.Quantity, qy: npt.NDArray[np.float64]
    ):
        """Initialize and load the quantum yield data.

        Args:
            molecule: The molecule for which the quantum yield is loaded.
            branch_id: The ID of the branching ratio.
            wavelength: The wavelength at which the quantum yield is measured.
            qy: The quantum yield value at the given wavelength.
        """
        self.wavelength = wavelength
        self.qy = qy
        self.molecule = molecule
        self.branch_id = branch_id

    def interp_to(self, wavelength: u.Quantity) -> "QuantumYield":
        """Interpolate the quantum yield to the given wavelength.

        Args:
            wavelength: The wavelength to which the quantum yield is interpolated.

        Returns:
            QuantumYield: A new QuantumYield object with the interpolated values.

        """

        new_wavelength = self.wavelength.to(wavelength.unit, equivalencies=u.spectral())

        return QuantumYield(
            self.molecule,
            self.branch_id,
            wavelength,
            np.interp(wavelength.value, new_wavelength.value, self.qy, left=0.0, right=0.0),
        )

__init__(molecule, branch_id, wavelength, qy)

Initialize and load the quantum yield data.

Parameters:

Name	Type	Description	Default
molecule	SpeciesFormula	The molecule for which the quantum yield is loaded.	required
branch_id	int | str	The ID of the branching ratio.	required
wavelength	Quantity	The wavelength at which the quantum yield is measured.	required
qy	NDArray[float64]	The quantum yield value at the given wavelength.	required

Source code in src/freckll/reactions/photo.py

def __init__(
    self, molecule: SpeciesFormula, branch_id: int | str, wavelength: u.Quantity, qy: npt.NDArray[np.float64]
):
    """Initialize and load the quantum yield data.

    Args:
        molecule: The molecule for which the quantum yield is loaded.
        branch_id: The ID of the branching ratio.
        wavelength: The wavelength at which the quantum yield is measured.
        qy: The quantum yield value at the given wavelength.
    """
    self.wavelength = wavelength
    self.qy = qy
    self.molecule = molecule
    self.branch_id = branch_id

interp_to(wavelength)

Interpolate the quantum yield to the given wavelength.

Parameters:

Name	Type	Description	Default
wavelength	Quantity	The wavelength to which the quantum yield is interpolated.	required

Returns:

Name	Type	Description
QuantumYield	QuantumYield	A new QuantumYield object with the interpolated values.

Source code in src/freckll/reactions/photo.py

def interp_to(self, wavelength: u.Quantity) -> "QuantumYield":
    """Interpolate the quantum yield to the given wavelength.

    Args:
        wavelength: The wavelength to which the quantum yield is interpolated.

    Returns:
        QuantumYield: A new QuantumYield object with the interpolated values.

    """

    new_wavelength = self.wavelength.to(wavelength.unit, equivalencies=u.spectral())

    return QuantumYield(
        self.molecule,
        self.branch_id,
        wavelength,
        np.interp(wavelength.value, new_wavelength.value, self.qy, left=0.0, right=0.0),
    )

StarSpectra

Source code in src/freckll/reactions/photo.py

class StarSpectra:
    def __init__(self, wavelength: u.Quantity, flux: u.Quantity, reference_distance: u.Quantity):
        """Initialize the star spectrum with wavelength and flux data.

        Args:
            wavelength: The wavelength of the star spectrum.
            flux: The flux of the star spectrum.
            reference_distance: The reference distance for the flux.

        """
        self.wavelength = wavelength
        self.flux = flux

        self.reference_distance = reference_distance

    def incident_flux(self, distance: u.Quantity) -> u.Quantity:
        r"""Calculate the flux density at the reference distance.


        The flux density is calculated using the formula:
        $$
        F = F_0(\frac{D_0}{d})^2$$
        where:
        - $F$ is the flux density at distance $d$
        - $F_0$ is the flux density at reference distance $D_0$
        - $d$ is the distance from the star
        - $D_0$ is the reference distance


        Returns:
            u.Quantity: The flux density at the reference distance.

        """
        return self.flux * (self.reference_distance / distance) ** 2 / 2

__init__(wavelength, flux, reference_distance)

Initialize the star spectrum with wavelength and flux data.

Parameters:

Name	Type	Description	Default
wavelength	Quantity	The wavelength of the star spectrum.	required
flux	Quantity	The flux of the star spectrum.	required
reference_distance	Quantity	The reference distance for the flux.	required

Source code in src/freckll/reactions/photo.py

def __init__(self, wavelength: u.Quantity, flux: u.Quantity, reference_distance: u.Quantity):
    """Initialize the star spectrum with wavelength and flux data.

    Args:
        wavelength: The wavelength of the star spectrum.
        flux: The flux of the star spectrum.
        reference_distance: The reference distance for the flux.

    """
    self.wavelength = wavelength
    self.flux = flux

    self.reference_distance = reference_distance

incident_flux(distance)

Calculate the flux density at the reference distance.

The flux density is calculated using the formula: $$ F = F_0(\frac{D_0}{d})^2$$ where: - \(F\) is the flux density at distance \(d\) - \(F_0\) is the flux density at reference distance \(D_0\) - \(d\) is the distance from the star - \(D_0\) is the reference distance

Returns:

Type	Description
Quantity	u.Quantity: The flux density at the reference distance.

Source code in src/freckll/reactions/photo.py

def incident_flux(self, distance: u.Quantity) -> u.Quantity:
    r"""Calculate the flux density at the reference distance.


    The flux density is calculated using the formula:
    $$
    F = F_0(\frac{D_0}{d})^2$$
    where:
    - $F$ is the flux density at distance $d$
    - $F_0$ is the flux density at reference distance $D_0$
    - $d$ is the distance from the star
    - $D_0$ is the reference distance


    Returns:
        u.Quantity: The flux density at the reference distance.

    """
    return self.flux * (self.reference_distance / distance) ** 2 / 2

optical_depth(altitude, number_density, cross_sections, cross_section_indices)

Compute the optical depth of the atmosphere.

The optical depth is computed using the formula:

\[ \tau = \int n(z) \sigma(z) dz \]

where:

• \(\tau\) is the optical depth
• \(n(z)\) is the number density of the species at altitude \(z\)
• \(\sigma(z)\) is the cross-section of the species at altitude \(z\)

The integral is computed using the trapezoidal rule.

Source code in src/freckll/reactions/photo.py

def optical_depth(
    altitude: u.Quantity,
    number_density: u.Quantity,
    cross_sections: u.Quantity,
    cross_section_indices: npt.NDArray[np.integer],
) -> FreckllArray:
    r"""Compute the optical depth of the atmosphere.

    The optical depth is computed using the formula:

    $$
    \tau = \int n(z) \sigma(z) dz
    $$

    where:

    - $\tau$ is the optical depth
    - $n(z)$ is the number density of the species at altitude $z$
    - $\sigma(z)$ is the cross-section of the species at altitude $z$

    The integral is computed using the trapezoidal rule.



    """
    dz = np.zeros_like(altitude)
    dz[:-1] = np.diff(altitude)
    dz[-1] = dz[-2]

    number_density_dz = number_density * dz

    cross_section_density = cross_sections[:, None, :] * number_density_dz[cross_section_indices, :, None]

    tau = np.sum(cross_section_density, axis=0).decompose()

    return tau

radiative_transfer(flux_top, optical_depth, incident_angle, albedo=0.0)

Computes the radiative transfer equation.

Computes the radiative transfer equation for a given flux and optical depth.

The flux is then propagated upwards and downwards through the atmosphere. The albedo is taken into account for the upward flux.

Parameters:

Name	Type	Description	Default
flux_top	Quantity	The flux at the top of the atmosphere.	required
optical_depth	FreckllArray	The optical depth of the atmosphere.	required
incident_angle	Quantity	The angle of incidence.	required
albedo	float | NDArray[float64]	The albedo of the surface.	0.0

Returns: u.Quantity: The total flux at each layer.

Source code in src/freckll/reactions/photo.py

def radiative_transfer(
    flux_top: u.Quantity,
    optical_depth: FreckllArray,
    incident_angle: u.Quantity,
    albedo: float | npt.NDArray[np.float64] = 0.0,
) -> u.Quantity:
    r"""Computes the radiative transfer equation.

    Computes the radiative transfer equation for a given flux and optical depth.

    The flux is then propagated upwards and downwards through the atmosphere.
    The albedo is taken into account for the upward flux.


    Args:
        flux_top: The flux at the top of the atmosphere.
        optical_depth: The optical depth of the atmosphere.
        incident_angle: The angle of incidence.
        albedo: The albedo of the surface.
    Returns:
        u.Quantity: The total flux at each layer.

    """

    tau = np.exp(-optical_depth / np.cos(incident_angle))

    flux_down = np.zeros_like(tau.value) << flux_top.unit
    flux_up = np.zeros_like(tau.value) << flux_top.unit

    flux_down[-1] = flux_top * tau[-1]
    num_layers = flux_down.shape[0]
    # Propagate the flux downwards from toa to boa
    for layer in reversed(range(num_layers - 1)):
        flux_down[layer] = flux_down[layer + 1] * tau[layer]

    # If albedo then reflect the flux
    flux_up[0] = flux_down[0] * albedo

    # Propagate the flux upwards from boa to toa
    for layer in range(1, num_layers):
        flux_up[layer] = flux_up[layer - 1] * tau[layer - 1]

    flux_down2 = np.zeros_like(flux_down)
    flux_up2 = np.zeros_like(flux_up)

    # Propagate the flux downwards from toa to boa
    flux_down2[-1] = (0.5 * ((1 - tau[-1]) * (flux_up[-1] + flux_top)) + flux_top) * tau[-1]

    flux_down2[-1] = (0.5 * ((1 - tau[-1]) * (flux_up[-1] + flux_top)) + flux_top) * tau[-1]
    for layer in reversed(range(num_layers)[:-1]):
        flux_down2[layer] = (
            0.5 * ((1 - tau[layer]) * (flux_up[layer] + flux_down[layer + 1])) + flux_down2[layer + 1]
        ) * tau[layer]

    flux_up2[0] = albedo * flux_down2[0]
    for layer in range(num_layers)[1:]:
        flux_up2[layer] = (
            0.5 * ((1 - tau[layer - 1]) * (flux_up[layer - 1] + flux_down[layer])) + flux_up2[layer - 1]
        ) * tau[layer - 1]

    return flux_up2 + flux_down2

rayleigh(spectral_grid, alpha, depolar_factor)

Calculate the Rayleigh scattering cross-section.

The Rayleigh scattering cross-section is given by the formula: $$ \sigma_R = \frac{8 \pi^3}{3} \left(\frac{\alpha}{\lambda}\right)^4 (1 + \delta^2)$$

where:

• \(\sigma_R\) is the Rayleigh scattering cross-section
• \(\alpha\) is the polarizability of the molecule

Parameters:

Name	Type	Description	Default
spectral_grid	Quantity	The wavelength grid.	required
alpha	Quantity	The polarizability of the molecule.	required
depolar_factor	float	The depolarization factor.	required

Returns: u.Quantity: The Rayleigh scattering cross-section.

Source code in src/freckll/reactions/photo.py

def rayleigh(spectral_grid: u.Quantity, alpha: u.Quantity, depolar_factor: float) -> u.Quantity:
    r"""Calculate the Rayleigh scattering cross-section.

    The Rayleigh scattering cross-section is given by the formula:
    $$
    \sigma_R = \frac{8 \pi^3}{3} \left(\frac{\alpha}{\lambda}\right)^4 (1 + \delta^2)$$

    where:

    - $\sigma_R$ is the Rayleigh scattering cross-section
    - $\alpha$ is the polarizability of the molecule

    Args:
        spectral_grid: The wavelength grid.
        alpha: The polarizability of the molecule.
        depolar_factor: The depolarization factor.
    Returns:
        u.Quantity: The Rayleigh scattering cross-section.

    """

    wavelength_cm = spectral_grid.to(u.cm, equivalencies=u.spectral())

    alpha = alpha.to(u.cm**3)

    return 8 * np.pi / 3 * (2 * np.pi / wavelength_cm) ** 4 * alpha**2 * depolar_factor