# AST3220, Cosmo I, project 2 -- Big Bang Nucleosynthesis
# Jakob Borg, UiO, 2025
# Module providing the reaction rates for the BBN calculations.
# From "ON THE SYNTHESIS OF ELEMENTS AT VERY HIGH TEMPERATURES"
# by R.V.Wagoner et al, 1967

# Imports:
...


class ReactionRates:
    def __init__(self, *args, **kwargs):
        """
        Initialize the reaction rates for the BBN calculations. These does not
        necessarily require any input parameters, but can be extended if needed.
        """
        # Define constants and parameters:
        self.tau = ...
        self.q = ...
        ...


    # -----------------------------------------------------------------------
    # Reaction for the 2 first species (n and p):
    # -----------------------------------------------------------------------
    def get_weak_force(self, T9):
        """
        Computes decay rates for the neutron proton weak force interactions
        Weak a. 1-3 \\
        n + nu_e <-> p + e \\
        n + e+ <-> p + nu_e_anti \\
        n <-> p + e + nu_e_anti

        Parameters
        ----------
        T9 : float or np.ndarray
            The temperature in units of 10^9 Kelvin.

        Returns
        -------
        float or np.ndarray
            The weak interaction rate in 1/s.
        """
        rate_n_to_p = ...
        rate_p_to_n = ...
        return rate_n_to_p, rate_p_to_n
    
    # -----------------------------------------------------------------------
    # Reaction for up to the 3 first species (n, p and D):
    # -----------------------------------------------------------------------
    def get_np_to_D(self, T9, rho_b):
        """
        Strong and Electromagnetic b.1 \\
        n + p <-> D + gamma

        Parameters
        ----------
        T9 : float or np.ndarray
            The temperature in units of 10^9 Kelvin.
        rho_b : float or np.ndarray
            The baryon density in g/cm^3.

        Returns
        -------
        float or np.ndarray
            The rate for the reaction in 1/s.
        """
        rate_np = ...
        rate_D = ...
        return rate_np, rate_D
    
    # -----------------------------------------------------------------------
    # Reaction for up to the 4 first species (n, p, D and T):
    # -----------------------------------------------------------------------
    def get_nD_to_T(self, T9, rho_b):
        """
        Strong and Electromagnetic b.3 \\
        n + D <-> T + gamma
        """
        rate_nD = ...
        rate_T = ...
        return rate_nD, rate_T
    
    def get_DD_to_pT(self, T9, rho_b):
        """
        # Strong and Electromagnetic b.8 \\
        # D + D <-> p + T
        """
        rate_DD = ...
        rate_pT = ...
        return rate_DD, rate_pT
    
    # -----------------------------------------------------------------------
    # Reaction for up to the 5 first species (n, p, D, T and He3):
    # -----------------------------------------------------------------------
    # Fill in the missing reaction rates here
    
    ...