Multiband emission

mlpoppyns.simulator.multiband_emission.emission_radio

Models for the pulsars' radio beam geometry and luminosity.

Authors:

Michele Ronchi (ronchi@ice.csic.es)
Celsa Pardo Araujo (pardo@ice.csic.es)

beam_aperture(P)

Half angular aperture in [rad] of the radio beam as a function of the spin period.

For the "standard_period_cone" model, the aperture is derived by assuming that the radio beam width extends into the open field line region around the magnetic poles of the star. See eq. (3.29) in Lorimer and Kramer (2005) and eq. (2) in Johnston et al. (2020) for a derivation and discussion of this model. For the "power-law_period_cone" model, we adopted a generic power-law of the spin period (see Maciesiak et al. 2012).

Parameters:

Name	Type	Description	Default
P	ndarray	Array of spin periods of the pulsars in [s].	required

Returns:

Type	Description
ndarray	Half angular aperture of the radio beam in [rad].

Source code in mlpoppyns/simulator/multiband_emission/emission_radio.py

def beam_aperture(P: np.ndarray) -> np.ndarray:
    """
    Half angular aperture in [rad] of the radio beam as a function of the spin period.

    For the "standard_period_cone" model, the aperture is derived by assuming that the radio beam width extends
    into the open field line region around the magnetic poles of the star.
    See eq. (3.29) in Lorimer and Kramer (2005) and eq. (2) in Johnston et al. (2020)
    for a derivation and discussion of this model.
    For the "power-law_period_cone" model, we adopted a generic power-law of the spin period (see Maciesiak et al. 2012).

    Args:
        P (np.ndarray): Array of spin periods of the pulsars in [s].

    Returns:
        (np.ndarray): Half angular aperture of the radio beam in [rad].
    """

    beam_model = cfg["radio_beam_model"]

    if beam_model == "standard_period_cone":
        rho_b = np.sqrt(9.0 * np.pi * cfg["r_em"] / (2.0 * const.C * P))
    elif beam_model == "power-law_period_cone":
        rho_b = cfg["rho_b_0"] * P ** cfg["a_beam"]
        # Convert the half angular aperture from [deg] to [rad].
        rho_b = rho_b * np.pi / 180.0
    else:
        raise ValueError(
            "The radio beam model does not exist. Choose between standard_period_cone or power-law_period_cone."
        )

    return rho_b

beam_fraction(chi, rho_b)

Fraction of solid angle covered by the radio beam in an entire pulsar rotation as a function of the inclination angle chi and the radio beam aperture. This is also the probability that the radio beam intercepts our line of sight. See Appendix in Emmering and Chevalier (1989).

Parameters:

Name	Type	Description	Default
chi	ndarray	Array of inclination angles of the magnetic axis with respect to the rotation axis of the pulsars in [rad].	required
rho_b	ndarray	Array of angular apertures of the radio beam of the pulsars in [rad].	required

Returns:

Type	Description
ndarray	Fraction of solid angle swept by the radio beam.

Source code in mlpoppyns/simulator/multiband_emission/emission_radio.py

def beam_fraction(chi: np.ndarray, rho_b: np.ndarray) -> np.ndarray:
    """
    Fraction of solid angle covered by the radio beam in an entire pulsar rotation
    as a function of the inclination angle chi and the radio beam aperture.
    This is also the probability that the radio beam intercepts our line of sight.
    See Appendix in Emmering and Chevalier (1989).

    Args:
        chi (np.ndarray): Array of inclination angles of the magnetic axis
            with respect to the rotation axis of the pulsars in [rad].
        rho_b (np.ndarray): Array of angular apertures of the radio beam of the pulsars in [rad].

    Returns:
        (np.ndarray): Fraction of solid angle swept by the radio beam.
    """

    beam_frac = np.cos(np.maximum(0, (chi - rho_b))) - np.cos(
        np.minimum(np.pi / 2.0, (chi + rho_b))
    )

    return beam_frac

calculate_radio_emission(P, P_dot, chi, dist)

Compute the radio beam geometry, the intrinsic bolometric radio flux and the DM. Note that the luminosity and DM are computed only for those pulsars whose beams cross our line of sight. This function is used in the simulate_population_magrot_det.py script.

Parameters:

Name	Type	Description	Default
P	ndarray	Array of spin periods of the pulsars in [s].	required
P_dot	ndarray	Array of neutron star spin period derivatives in [s s^-1].	required
chi	ndarray	Array of inclination angles of the pulsars in [rad].	required
dist	ndarray	Array of distances from the ICRS origin in [kpc].	required

Returns:

Type	Description
Tuple[ndarray, ndarray, ndarray, ndarray, ndarray, ndarray, ndarray]	A Tuple
ndarray	containing the following information: Boolean mask to select the pulsars whose radio beam crosses our line of sight. Intrinsic pulse widths in [s]. Bolometric radio luminosity [erg s^-1]. Bolometric radio flux in [erg s^-1 cm^-2]. Spectral index.

Source code in mlpoppyns/simulator/multiband_emission/emission_radio.py

def calculate_radio_emission(
    P: np.ndarray,
    P_dot: np.ndarray,
    chi: np.ndarray,
    dist: np.ndarray,
) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray,]:
    """
    Compute the radio beam geometry, the intrinsic bolometric radio flux and the DM.
    Note that the luminosity and DM are computed only for those pulsars whose beams cross our line of sight.
    This function is used in the simulate_population_magrot_det.py script.

    Args:
        P (np.ndarray): Array of spin periods of the pulsars in [s].
        P_dot (np.ndarray): Array of neutron star spin period derivatives in [s s^-1].
        chi (np.ndarray): Array of inclination angles of the pulsars in [rad].
        dist (np.ndarray): Array of distances from the ICRS origin in [kpc].

    Returns:
        (Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray]): A Tuple
        containing the following information:

            - Boolean mask to select the pulsars whose radio beam crosses our line of sight.
            - Intrinsic pulse widths in [s].
            - Bolometric radio luminosity [erg s^-1].
            - Bolometric radio flux in [erg s^-1 cm^-2].
            - Spectral index.
    """

    # Determining the bolometric radio luminosity.
    # Choose one of the two implementations either based on the P and Pdot or the Edot dependence.
    # See line 251 in the config_simulator.py file.
    radio_luminosity_model = cfg["radio_luminosity_model"]

    if radio_luminosity_model == "lum_radio_ppdot":
        L_radio_bol = pdf_luminosity_radio_ppdot(P, P_dot)
    elif radio_luminosity_model == "lum_radio_edot":
        L_radio_bol = pdf_luminosity_radio_edot(P, P_dot)
    else:
        raise ValueError(
            "The radio luminosity model does not exist. Choose between lum_radio_ppdot or lum_radio_edot."
        )

    # Determining the radio beam angular aperture.
    rho_beam = beam_aperture(P)

    # Determining the solid angle covered by the two radio beams.
    solid_angle_beam = solid_angle_radio_beams(rho_beam)

    # Drawing a random angular intercept for the line of sight.
    # Note that since we assume symmetry between the northern and southern hemisphere of the star
    # we only need to consider one hemisphere, e.g., the northern one.
    los_grid = np.linspace(0.0, np.pi / 2, cfg["resolution"])
    los_rand = rs.random_from_pdf(los_grid, np.sin, len(P))

    # Determining if the pulsar's radio beam intercepts our line of sight.
    intercepted_radio = los_intercept(
        chi,
        rho_beam,
        los_rand,
    )

    # Computing the intrinsic pulse width of the radio pulse.
    w_int = np.zeros(len(P))
    w_int[intercepted_radio] = pulse_width(
        chi[intercepted_radio],
        rho_beam[intercepted_radio],
        los_rand[intercepted_radio],
    )

    # Computing the intrinsic bolometric radio flux.
    S_radio_bol = np.zeros(len(P))
    S_radio_bol[intercepted_radio] = flux_radio(
        L_radio_bol[intercepted_radio],
        dist[intercepted_radio],
        solid_angle_beam[intercepted_radio],
    )

    # Convert pulse width from [rad] to [s].
    w_int_s = w_int * P / (2.0 * np.pi)

    # Computing the spectral index for the radio power-law spectrum.
    spectral_index = compute_spectral_index(
        cfg["mean_spectral_index"],
        cfg["std_spectral_index"],
        len(P),
    )

    return (
        intercepted_radio,
        w_int_s,
        L_radio_bol,
        S_radio_bol,
        spectral_index,
    )

compute_spectral_index(mean, sigma, NS_number)

Draw a random spectral index from a Gaussian distribution (see Posselt et al. 2023).

Parameters:

Name	Type	Description	Default
mean	float	Mean spectral index for the Gaussian distribution.	required
sigma	float	Standard deviation for the Gaussian distribution.	required
NS_number	int	Number of neutron stars for which sampling the spectral index.	required

Returns:

Type	Description
ndarray	Array of spectral indices.

Source code in mlpoppyns/simulator/multiband_emission/emission_radio.py

def compute_spectral_index(
    mean: float, sigma: float, NS_number: int
) -> np.ndarray:
    """
    Draw a random spectral index from a Gaussian distribution (see Posselt et al. 2023).

    Args:
        mean (float): Mean spectral index for the Gaussian distribution.
        sigma (float): Standard deviation for the Gaussian distribution.
        NS_number (int): Number of neutron stars for which sampling the spectral index.

    Returns:
        (np.ndarray): Array of spectral indices.
    """

    spectral_index = np.random.normal(
        mean,
        sigma,
        NS_number,
    )

    return spectral_index

flux_density_radio(S_radio_bol, spectral_index, f, f_min=10000000.0, f_max=100000000000.0)

Compute the radio flux density at a given frequency f assuming a power law spectral shape for the radio emission.

Parameters:

Name	Type	Description	Default
S_radio_bol	ndarray	pulsar bolometric radio flux in [erg s^(-1) cm^(-2)].	required
spectral_index	float	spectral index of the radio emission, assuming a power-law spectrum.	required
f	float	frequency in [Hz] at which the radio luminosity has to be computed.	required
f_min	float	frequency lower limit of the radio emission spectrum [Hz].	10000000.0
f_max	float	frequency upper limit of the radio emission spectrum [Hz].	100000000000.0

Returns:

Type	Description
ndarray	Intrinsic pulsar radio flux density in [Jy] at the frequency f.

Source code in mlpoppyns/simulator/multiband_emission/emission_radio.py

def flux_density_radio(
    S_radio_bol: np.ndarray,
    spectral_index: np.ndarray,
    f: float,
    f_min: float = 1.0e7,
    f_max: float = 1.0e11,
) -> np.ndarray:
    """
    Compute the radio flux density at a given frequency f assuming a power law spectral shape for the radio emission.

    Args:
        S_radio_bol (np.ndarray): pulsar bolometric radio flux in [erg s^(-1) cm^(-2)].
        spectral_index (float): spectral index of the radio emission,
            assuming a power-law spectrum.
        f (float): frequency in [Hz] at which the radio luminosity has to be computed.
        f_min (float): frequency lower limit of the radio emission spectrum [Hz].
        f_max (float): frequency upper limit of the radio emission spectrum [Hz].

    Returns:
        (np.ndarray): Intrinsic pulsar radio flux density in [Jy] at the frequency f.
    """

    S_radio_f = (
        (spectral_index + 1)
        * S_radio_bol
        / (f_max ** (spectral_index + 1) - f_min ** (spectral_index + 1))
        * f**spectral_index
    ) / const.JY_TO_ERG

    return S_radio_f

flux_radio(L_radio, d, solid_angle)

Compute the intrinsic bolometric radio flux for each pulsars in [erg s^(-1) cm^(-2)].

Parameters:

Name	Type	Description	Default
L_radio	ndarray	Pulsar bolometric radio luminosity in [erg s^(-1)].	required
d	ndarray	Distance to the pulsar in [kpc].	required
solid_angle	ndarray	Solid angle covered by the radio beams in [sr].	required

Returns:

Type	Description
ndarray	Intrinsic pulsar bolometric radio flux in [erg s^(-1) cm^(-2)].

Source code in mlpoppyns/simulator/multiband_emission/emission_radio.py

def flux_radio(
    L_radio: np.ndarray,
    d: np.ndarray,
    solid_angle: np.ndarray,
) -> np.ndarray:
    """
    Compute the intrinsic bolometric radio flux for each pulsars in [erg s^(-1) cm^(-2)].

    Args:
        L_radio (np.ndarray): Pulsar bolometric radio luminosity in [erg s^(-1)].
        d (np.ndarray): Distance to the pulsar in [kpc].
        solid_angle (np.ndarray): Solid angle covered by the radio beams in [sr].

    Returns:
        (np.ndarray): Intrinsic pulsar bolometric radio flux in [erg s^(-1) cm^(-2)].
    """

    # Convert distance from [kpc] to [cm].
    d_cm = d * const.KPC_TO_CM

    S_radio = L_radio / (solid_angle * d_cm**2)

    return S_radio

los_intercept(chi, rho_b, los)

Evaluate if the radio beam intercepts the line of sight (LOS), assuming random orientation of the LOS with respect to the rotation axis.

Parameters:

Name	Type	Description	Default
chi	ndarray	Array of inclination angles of the pulsars in [rad].	required
rho_b	ndarray	Array of angular apertures of the radio beam of the pulsars in [rad].	required
los	ndarray	Polar angle of the line of sight intercept computed with respect to the rotation axis of the star [rad].	required

Returns:

Type	Description
ndarray	Array of boolean variables: true if the radio beam intercepts the LOS and false if not.

Source code in mlpoppyns/simulator/multiband_emission/emission_radio.py

def los_intercept(
    chi: np.ndarray, rho_b: np.ndarray, los: np.ndarray
) -> np.ndarray:
    """
    Evaluate if the radio beam intercepts the line of sight (LOS), assuming random
    orientation of the LOS with respect to the rotation axis.

    Args:
        chi (np.ndarray): Array of inclination angles of the pulsars in [rad].
        rho_b (np.ndarray): Array of angular apertures of the radio beam of the pulsars in [rad].
        los (np.ndarray): Polar angle of the line of sight intercept computed with respect
            to the rotation axis of the star [rad].

    Returns:
        (np.ndarray): Array of boolean variables: true if the radio beam intercepts the LOS and false if not.
    """

    condition = (los > chi - rho_b) & (los < chi + rho_b)
    intercepted = np.array(condition)

    return intercepted

loss_rotational_energy(P, P_dot)

Compute the loss of the rotational energy given the period and period derivative (see, e.g., eq. (3.5) in Lorimer and Kramer, 2004).

Parameters:

Name	Type	Description	Default
P	ndarray	array of spin periods of the pulsars in [s].	required
P_dot	ndarray	array of spin period derivatives of the pulsars in [s/s].	required

Returns:

Type	Description
ndarray	Loss of the rotational energy [erg s^(-1)].

Source code in mlpoppyns/simulator/multiband_emission/emission_radio.py

def loss_rotational_energy(P: np.ndarray, P_dot: np.ndarray) -> np.ndarray:
    """
    Compute the loss of the rotational energy given the period and period derivative
    (see, e.g., eq. (3.5) in Lorimer and Kramer, 2004).

    Args:
        P (np.ndarray): array of spin periods of the pulsars in [s].
        P_dot (np.ndarray): array of spin period derivatives of the pulsars in [s/s].

    Returns:
        (np.ndarray): Loss of the rotational energy [erg s^(-1)].
    """
    NS_inertia = 2.0 / 5.0 * cfg["NS_mass"] * cfg["NS_radius"] ** 2
    Erot_dot = NS_inertia * (2.0 * np.pi) ** 2 * (P_dot / (P**3))

    return Erot_dot

pdf_luminosity_radio_edot(P, P_dot)

Draw random bolometric radio luminosities from a distribution that depends on the loss of the rotational energy. We assume a random log-normal spread for the normalization constant L_0.

Parameters:

Name	Type	Description	Default
P	ndarray	array of spin periods of the pulsars in [s].	required
P_dot	ndarray	array of spin period derivatives of the pulsars in [s/s].	required

Returns:

Type	Description
ndarray	pulsar radio luminosity [erg s^(-1)] drawn from a log-normal distribution.

Source code in mlpoppyns/simulator/multiband_emission/emission_radio.py

def pdf_luminosity_radio_edot(P: np.ndarray, P_dot: np.ndarray) -> np.ndarray:
    """
    Draw random bolometric radio luminosities from a distribution that depends on the loss of the rotational energy.
    We assume a random log-normal spread for the normalization constant L_0.

    Args:
        P (np.ndarray): array of spin periods of the pulsars in [s].
        P_dot (np.ndarray): array of spin period derivatives of the pulsars in [s/s].

    Returns:
        (np.ndarray): pulsar radio luminosity [erg s^(-1)] drawn from a log-normal distribution.
    """
    NS_number = len(P)
    L_0 = 10 ** np.random.normal(
        cfg["L_radio_log10_mean"],
        cfg["L_radio_log10_sigma"],
        NS_number,
    )
    Erot_dot = loss_rotational_energy(P, P_dot)
    L_radio = L_0 * (Erot_dot / cfg["Erot_dot_0"]) ** cfg["epsilon_L"]

    return L_radio

pdf_luminosity_radio_ppdot(P, P_dot)

Draw random bolometric radio luminosities from a distribution that depends on the spin period and spin period derivative. We assume a random log-normal spread for the normalization constant L_0.

Parameters:

Name	Type	Description	Default
P	ndarray	Array of spin periods of the pulsars in [s].	required
P_dot	ndarray	Array of spin period derivatives of the pulsars in [s/s].	required

Returns:

Type	Description
ndarray	Pulsar radio luminosity [erg s^(-1)] drawn from a log-normal distribution.

Source code in mlpoppyns/simulator/multiband_emission/emission_radio.py

def pdf_luminosity_radio_ppdot(P: np.ndarray, P_dot: np.ndarray) -> np.ndarray:
    """
    Draw random bolometric radio luminosities from a distribution that depends on the spin period
    and spin period derivative. We assume a random log-normal spread for the normalization constant L_0.

    Args:
        P (np.ndarray): Array of spin periods of the pulsars in [s].
        P_dot (np.ndarray): Array of spin period derivatives of the pulsars in [s/s].

    Returns:
        (np.ndarray): Pulsar radio luminosity [erg s^(-1)] drawn from a log-normal distribution.
    """

    NS_number = len(P)

    L_0 = 10 ** np.random.normal(
        cfg["L_radio_log10_mean"],
        cfg["L_radio_log10_sigma"],
        NS_number,
    )
    L_radio = L_0 * (P ** (-3) * P_dot) ** cfg["epsilon_L"]

    return L_radio

pulse_width(chi, rho_b, los)

Formula to evaluate the intrinsic pulse width in [rad] from the radio beam angular aperture. This assumes that the line of sight intercepts the radio beam with an angle beta with respect to the center of the beam, so that -rho_b < beta < rho_b, with rho_b the angular aperture of the beam. See eq. (1) in Maciesiak et al. (2011a) and eq. (3.26) in Lorimer and Kramer (2004).

Parameters:

Name	Type	Description	Default
chi	ndarray	Array of inclination angles between the magnetic axis and the rotation axis [rad].	required
rho_b	ndarray	Array of angular apertures of the radio beam of the pulsars in [rad].	required
los	ndarray	Polar angles of the line of sight intercept computed with respect to the rotation axis of the star [rad].	required

Returns:

Type	Description
ndarray	Array of intrinsic pulse widths in [rad].

Source code in mlpoppyns/simulator/multiband_emission/emission_radio.py

def pulse_width(
    chi: np.ndarray, rho_b: np.ndarray, los: np.ndarray
) -> np.ndarray:
    """
    Formula to evaluate the intrinsic pulse width in [rad] from the radio beam angular aperture.
    This assumes that the line of sight intercepts the radio beam with an angle beta with
    respect to the center of the beam, so that -rho_b < beta < rho_b, with rho_b the angular
    aperture of the beam. See eq. (1) in Maciesiak et al. (2011a) and eq. (3.26) in Lorimer and Kramer (2004).

    Args:
        chi (np.ndarray): Array of inclination angles between the magnetic axis and the rotation axis [rad].
        rho_b (np.ndarray): Array of angular apertures of the radio beam of the pulsars in [rad].
        los (np.ndarray): Polar angles of the line of sight intercept computed with respect
            to the rotation axis of the star [rad].

    Returns:
        (np.ndarray): Array of intrinsic pulse widths in [rad].
    """

    # Compute the impact parameter, i.e., the angular distance between the LOS intercept and the center of the beam.
    beta = los - chi

    # Compute the quantity sin(w/4) where w is the pulse width.
    sin_w_4 = np.array(
        np.sqrt(
            (np.sin(rho_b / 2.0) ** 2 - np.sin(beta / 2) ** 2)
            / (np.sin(chi + beta) * np.sin(chi))
        )
    )
    # If sin(w/4) > 1 set w = 2*np.pi.
    sin_w_4[sin_w_4 > 1] = 1.0

    w_int = 4.0 * np.arcsin(sin_w_4)

    return w_int

radio_population_intercepted(dict_pop, coverage_radio, full_population=False)

Filter and compute properties of a population of neutron stars whose radio beams intercept our line of sight and that fall in the survey sky coverage.

Parameters:

Name	Type	Description	Default
dict_pop	dict	Dictionary containing the properties of a neutron star population.	required
coverage_radio	ndarray	A boolean mask to filter only stars in the sky coverage of radio surveys.	required
full_population	bool	If True, return arrays of size cfg["NS_number"] (all stars, filling with zeros where not intercepted). If False, return arrays only for intercepted stars. Defaults to False.	False

Returns:

Type	Description
dict	A dictionary containing properties of the neutron stars whose radio beams intercept our line of sight.

Source code in mlpoppyns/simulator/multiband_emission/emission_radio.py

def radio_population_intercepted(
    dict_pop: dict,
    coverage_radio: np.ndarray,
    full_population: bool = False,
) -> dict:
    """
    Filter and compute properties of a population of neutron stars whose radio beams intercept our line of sight
    and that fall in the survey sky coverage.

    Args:
        dict_pop (dict): Dictionary containing the properties of a neutron star population.
        coverage_radio (np.ndarray): A boolean mask to filter only stars in the sky coverage of radio surveys.
        full_population (bool, optional): If True, return arrays of size `cfg["NS_number"]` (all stars, filling
            with zeros where not intercepted). If False, return arrays only for intercepted stars.
            Defaults to False.

    Returns:
        (dict): A dictionary containing properties of the neutron stars whose radio beams intercept our line of sight.
    """
    if full_population:
        # Find the pulsars whose radio beam intercepts our line of sight and compute the intrinsic properties
        # of their radio emission.
        (
            intercepted_radio,
            w_int_s,
            L_radio_bol,
            S_radio_bol,
            spectral_index,
        ) = calculate_radio_emission(
            dict_pop["P"],
            dict_pop["P_dot"],
            dict_pop["chi"],
            dict_pop["dist"],
        )

        # Determine which stars could in principle be detected.
        detectable_radio = intercepted_radio & coverage_radio

        dictionary_radio_pop = {key: value for key, value in dict_pop.items()}

        dictionary_radio_pop["intercepted_radio"] = intercepted_radio

        # Computing the DM for the stars that fall into the surveys' sky coverage and whose
        # radio beam intercepts our line of sight.
        DM = np.zeros(cfg["NS_number"])
        DM[detectable_radio] = edm.compute_DM(
            dictionary_radio_pop["l"][detectable_radio],
            dictionary_radio_pop["b"][detectable_radio],
            dictionary_radio_pop["dist"][detectable_radio],
            cfg["ed_model"],
        )

        # Computing the scattering timescale at 327 MHz of these stars.
        tau_sc = np.zeros(cfg["NS_number"])
        tau_sc[DM != 0] = edm.compute_tau_sc_327(DM[DM != 0])

    else:
        # Select only the stars that can, in principle, be detected in radio, i.e., those that they lie within the
        # observed region.
        dict_pop_filtered = {
            key: value[coverage_radio] for key, value in dict_pop.items()
        }

        # Find the pulsars whose radio beam intercepts our line of sight and compute the intrinsic properties
        # of their radio emission.
        (
            intercepted_radio,
            w_int_s,
            L_radio_bol,
            S_radio_bol,
            spectral_index,
        ) = calculate_radio_emission(
            dict_pop_filtered["P"],
            dict_pop_filtered["P_dot"],
            dict_pop_filtered["chi"],
            dict_pop_filtered["dist"],
        )

        # Select only the stars that can be detected in radio by the considered surveys.
        dictionary_radio_pop = {
            key: value[intercepted_radio]
            for key, value in dict_pop_filtered.items()
        }

        w_int_s = w_int_s[intercepted_radio]
        L_radio_bol = L_radio_bol[intercepted_radio]
        S_radio_bol = S_radio_bol[intercepted_radio]
        spectral_index = spectral_index[intercepted_radio]

        # Computing the DM for the stars that fall into the surveys' sky coverage and whose
        # radio beam intercepts our line of sight.
        DM = edm.compute_DM(
            dictionary_radio_pop["l"],
            dictionary_radio_pop["b"],
            dictionary_radio_pop["dist"],
            cfg["ed_model"],
        )

        # Computing the scattering timescale at 327 MHz of these stars.
        tau_sc = edm.compute_tau_sc_327(DM)

    dictionary_radio_pop["w_int"] = w_int_s
    dictionary_radio_pop["L_radio_bol"] = L_radio_bol
    dictionary_radio_pop["S_radio_bol"] = S_radio_bol
    dictionary_radio_pop["spectral_index"] = spectral_index
    dictionary_radio_pop["DM"] = DM
    dictionary_radio_pop["tau_sc"] = tau_sc

    return dictionary_radio_pop

solid_angle_radio_beams(rho_b)

Total solid angle covered by the two radio beams as a function of the radio beam aperture.

Parameters:

Name	Type	Description	Default
rho_b	ndarray	Array of angular apertures of the radio beam of the pulsars in [rad].	required

Returns:

Type	Description
ndarray	Total solid angle covered by the radio beams.

Source code in mlpoppyns/simulator/multiband_emission/emission_radio.py

def solid_angle_radio_beams(rho_b: np.ndarray) -> np.ndarray:
    """
    Total solid angle covered by the two radio beams as a function of the radio beam aperture.

    Args:
        rho_b (np.ndarray): Array of angular apertures of the radio beam of the pulsars in [rad].

    Returns:
        (np.ndarray): Total solid angle covered by the radio beams.
    """

    solid_angle = 4 * np.pi * (1 - np.cos(rho_b))

    return solid_angle