Multiband surveys

mlpoppyns.simulator.multiband_surveys.survey_radio

Model for the pulsar radio surveys.

We consider the following surveys:

1) PMPS: the Parks Multibeam Pulsar Survey (see Manchester et al. 2001, Lorimer et al. 2006) 2) SMPS: the Swinburne Parkes Multibeam Pulsar Survey (see Edwards et al. 2001, Jacoby et al. 2009) 3) HTRU: the High Time Resolution Universe Survey (see Keith et al. 2010)

Authors:

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

SurveyRadio

Class for any radio survey with a Gaussian telescope beam pattern.

Source code in mlpoppyns/simulator/multiband_surveys/survey_radio.py

class SurveyRadio:
    """
    Class for any radio survey with a Gaussian telescope beam pattern.
    """

    def __init__(self, parameters_path: str) -> None:
        """
        Radio survey initialization.
        The parameters for the survey are imported from a JSON file.

        Args:
            parameters_path (str): Path to the survey_parameter.json file containing
                the parameters of the radio survey. The file must include:

                - deg_factor (float): Degradation factor.
                - G0 (float): Gain at the beam center [KJy ^ (-1)].
                - t_obs (float): Integration time [s].
                - t_samp (float): Sampling time [s].
                - T_sys (float): System temperature [K].
                - nu_central (float): Central frequency of the bandwidth [Hz].
                - BW (float): Frequency bandwidth [Hz].
                - channel_width (float): Width of a single frequency channel [Hz].
                - n_pol (float): Number of polarizations.
                - FWHM (float): FWHM of the beam[arcmin].
                - SNR_th (float): Threshold signal-to-noise ratio.
                - RA_range (np.ndarray): Range of the sky covered by the survey in RA [deg].
                - DEC_range(np.ndarray): Range of the sky covered by the survey in DEC [deg].
                - l_range(np.ndarray): Range of the sky covered by the survey in Galactic longitude l[deg].
                - b_range_abs(np.ndarray): Absolute value of the range of the sky covered
                    by the survey in Galactic latitude b [deg].
                - aperture_config (bool): If True, use the aperture array configuration for the pulsar detection.
        """

        # Load parameters from JSON file.
        with open(parameters_path) as read_file:
            self.parameters = json.load(read_file)

        # Save the parameters.
        self.deg_factor = self.parameters["deg_factor"]
        self.G0 = self.parameters["G0"]
        self.t_obs = self.parameters["t_obs"]
        self.t_samp = self.parameters["t_samp"]
        self.T_sys = self.parameters["T_sys"]
        self.f_central = self.parameters["f_central"]
        self.BW = self.parameters["BW"]
        self.channel_width = self.parameters["channel_width"]
        self.n_pol = self.parameters["n_pol"]
        self.FWHM = self.parameters["FWHM"]
        self.SNR_th = self.parameters["SNR_th"]
        self.RA_range = self.parameters["RA_range"]
        self.DEC_range = self.parameters["DEC_range"]
        self.l_range = self.parameters["l_range"]
        self.b_range_abs = self.parameters["b_range_abs"]
        self.aperture_config = self.parameters["aperture_config"]
        self.FFT_search = self.parameters["FFT_search"]
        self.name = self.parameters["name"]

    def sky_coverage(
        self,
        RA: np.ndarray,
        DEC: np.ndarray,
        l_gal: np.ndarray,
        b_gal: np.ndarray,
    ) -> np.ndarray:
        """
        Determine which neutron stars are in the sky region covered by the survey.

        Args:
            RA (np.ndarray): Right ascension in [deg] defined between [0, 360] deg in ICRS frame.
            DEC (np.ndarray): Declination in [deg] defined between [-90, 90] deg in ICRS frame.
            l_gal (np.ndarray): Galactic longitude in [deg] defined between [-180, 180] deg.
            b_gal (np.ndarray): Galactic latitude in [deg] defined between [-90, 90] deg.

        Returns:
            (np.ndarray): Array of boolean variables: True if the pulsar is in the covered sky region, False if not.
        """
        if self.name == "HTRU high":
            coverage = (
                (RA > self.RA_range[0])
                & (RA < self.RA_range[1])
                & (DEC > self.DEC_range[0])
                & (DEC < self.DEC_range[1])
                & (
                    (
                        (
                            (l_gal > self.l_range[0][0])
                            & (l_gal < self.l_range[0][1])
                        )
                        | (
                            (l_gal > self.l_range[1][0])
                            & (l_gal < self.l_range[1][1])
                        )
                    )
                    | (
                        (np.abs(b_gal) > self.b_range_abs[0])
                        & (np.abs(b_gal) < self.b_range_abs[1])
                    )
                )
            )

        else:
            coverage = (
                (RA > self.RA_range[0])
                & (RA < self.RA_range[1])
                & (DEC > self.DEC_range[0])
                & (DEC < self.DEC_range[1])
                & (l_gal > self.l_range[0])
                & (l_gal < self.l_range[1])
                & (np.abs(b_gal) > self.b_range_abs[0])
                & (np.abs(b_gal) < self.b_range_abs[1])
            )

        return coverage

    def detection_offset(self, n_detection: int) -> np.ndarray:
        """
        Generating a random offset with respect to the beam center for the detections.
        A Gaussian beam pattern is assumed to account for sensitivity decay for off-center detections.
        See Lorimer et al. (1993) and the paragraph following eq. (29) in Bates et al. (2014).

        Args:
            n_detection (int): Number of detections to simulate.

        Returns:
            (np.ndarray): Square of the offset from the beam center for each detection in [arcmin^2] .
        """
        offset2 = np.random.uniform(0.0, self.FWHM**2 / 4.0, n_detection)
        return offset2

    def gain_gaussian_beam(self, offset2: np.ndarray) -> np.ndarray:
        """
        This method simulates the gain pattern of a receiver.
        A Gaussian beam pattern is assumed to account for sensitivity decay for off-center detections.
        See Lorimer et al. (1993) and the paragraph following eq. (29) in Bates et al. (2014).

        Args:
            offset2 (np.ndarray): Squared offset from the beam center in [arcmin^2].

        Returns:
            (np.ndarray): Gain of the telescope for the given offset in [K Jy^(-1)].
        """
        G = self.G0 * np.exp(-2.77 * offset2 / self.FWHM**2)

        return G

    def aperture_array_factor(self, DEC: np.ndarray) -> np.ndarray:
        """
        Computes the aperture array sensitivity correction factor as a function of declination for each pulsar.
        Args:
            DEC (np.ndarray): Declination in [deg] defined between [-90, 90] deg in ICRS frame.
        Returns:
            (np.ndarray): Correction factor emulating the sensitivity of an aperture array for different declinations.
        """

        # Compute the pulsars' offset angles from zenith and convert them to [rad].
        offset_from_zenith = (
            DEC - (self.DEC_range[0] + self.DEC_range[1]) / 2.0
        ) * const.DEG_TO_RAD

        aa_factor = np.cos(offset_from_zenith)

        return aa_factor

    def fft_search_efficiency(self, duty_cycle: np.ndarray) -> np.ndarray:
        """
        Compute the efficiency factor from Morello et al. (2020) (see eq. 44) to account for incoherent FFT search.

        Args:
            duty_cycle (np.ndarray): duty cycle of pulsars.
        Returns:
            (np.ndarray): Correction factor emulating the sensitivity of an incoherent FFT search.
        """

        epsilon = (1.0 + 0.0473 * duty_cycle ** (-0.627)) ** (-1)

        return epsilon

    def radiometer_equation(
        self,
        S_radio_obs_mean: np.ndarray,
        G: np.ndarray,
        w_eff: np.ndarray,
        P: np.ndarray,
        T_sky: np.ndarray,
    ) -> np.ndarray:
        """
        Radiometer equation used to compute the signal-to-noise ratio of each pulsar given the observed
        period-averaged radio flux at a given frequency, the effective pulse width, the spin period and
        the survey parameters (see eq. (A1.22) in Lorimer & Kramer 2005). We are assuming a square pulse
        shape for simplicity with height equal to the observed flux and width equal to the effective width.

        Args:
            S_radio_obs_mean (np.ndarray): Observed period-averaged radio flux density in [Jy].
            G (np.ndarray): Gain of the telescope for the given detection in [K Jy^(-1)].
            w_eff (np.ndarray): Effective pulse width in [s].
            P (np.ndarray): Spin period in [s].
            T_sky (np.ndarray): Sky temperature for every detection in [K].

        Returns:
            (np.ndarray): Signal-to-noise ratio of the detection.
        """
        SNR = np.zeros(len(S_radio_obs_mean))

        # If the effective pulse width is larger than the spin period,
        # emission is continuous and the neutron star cannot be detected as a pulsar.
        cond = w_eff < P

        # Compute the SNR of each detection using the radiometer equation.
        SNR[cond] = (
            S_radio_obs_mean[cond]
            * G[cond]
            * np.sqrt(self.n_pol * self.t_obs * self.BW)
            * np.sqrt((P[cond] - w_eff[cond]) / w_eff[cond])
            / (self.deg_factor * (self.T_sys + T_sky[cond]))
        )

        return SNR

    def simulate_detection(
        self,
        S_radio_obs_mean: np.ndarray,
        l_gal: np.ndarray,
        b_gal: np.ndarray,
        DEC: np.ndarray,
        w_eff: np.ndarray,
        P: np.ndarray,
    ) -> np.ndarray:
        """
        Simulate a detection: If the measured SNR surpasses the threshold SNR_th of the survey
        then the pulsar is detected.

        Args:
            S_radio_obs_mean (np.ndarray): Observed period-averaged radio flux density in [Jy].
            l_gal (np.ndarray): Galactic longitude in [deg] defined between [-180, 180] deg.
            b_gal (np.ndarray): Galactic latitude in [deg] defined between [-90, 90] deg.
            DEC (np.ndarray): Declination in [deg] defined between [-90, 90] deg in ICRS frame.
            w_eff (np.ndarray): Effective pulse width in [s].
            P (np.ndarray): Spin period in [s].

        Returns:
            (np.ndarray): Array of boolean variables: True if the pulsar is detected, False if not.
        """

        # Store the total number of sources.
        n = len(S_radio_obs_mean)

        # Draw a random offset from the telescope beam center.
        offset2 = self.detection_offset(n)
        # Compute the gain corresponding to the offset detections.
        G = self.gain_gaussian_beam(offset2)

        # Compute the sky temperature in the coordinates of each detection at the central frequency of the survey.
        # We choose here to use the refined map from Remazeilles et al. (2015).
        T_sky = sky_temperature_H81refined(l_gal, b_gal, self.f_central)
        SNR_detection = self.radiometer_equation(
            S_radio_obs_mean, G, w_eff, P, T_sky
        )

        if self.aperture_config:
            aa_factor = self.aperture_array_factor(DEC)
            SNR_detection = SNR_detection * aa_factor

        if self.FFT_search:
            # Apply efficiency factor from Morello et al. (2020) (see eq. 44) to account for incoherent FFT search.
            duty_cycle = w_eff / P
            epsilon = self.fft_search_efficiency(duty_cycle)
            SNR_detection = SNR_detection * epsilon

        detected = SNR_detection > self.SNR_th

        return detected

    def detected_radio_population(
        self,
        w_int_s: np.ndarray,
        DM: np.ndarray,
        P: np.ndarray,
        intercepted_radio: np.ndarray,
        coverage: np.ndarray,
        l_gal: np.ndarray,
        b_gal: np.ndarray,
        DEC: np.ndarray,
        S_radio_bol: np.ndarray,
        spectral_index: np.ndarray,
        tau_sc: np.ndarray,
    ) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
        """
        Compute the pulsars detected by each survey.

        Args:
            w_int_s (np.ndarray): Intrinsic pulse widths in [s]
            DM (np.ndarray): Dispersion measure in [pc cm^-3].
            P (np.ndarray): Array of spin periods of the pulsars in [s].
            intercepted_radio: (np.ndarray) Array of boolean variables where True values represent stars whose
                beam crosses our line of sight.
            coverage (np.ndarray): Array of boolean variables indicating the pulsars within the sky coverage.
            l_gal (np.ndarray): Galactic longitude in [deg] defined between [-180, 180] deg.
            b_gal (np.ndarray): Galactic latitude in [deg] defined between [-90, 90] deg.
            DEC (np.ndarray): Declination in [deg] defined between [-90, 90] deg in ICRS frame.
            S_radio_bol (np.ndarray): Pulsar bolometric radio flux in [erg s^(-1) cm^(-2)].
            spectral_index (np.ndarray): Spectral indexes.
            tau_sc (np.ndarray): Scattering timescale in [s].

        Returns:
            (Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]): Tuple consisting of the following arrays:

                - Boolean mask to select the pulsars detected by the survey.
                - Effective pulse width in [s].
                - Period-averaged fluxes at the central frequency of the survey in [Jy]
                - Period-averaged fluxes at the central frequency of 1.429 GHz in [Jy].
        """

        detectable_radio_survey = intercepted_radio & coverage

        # Computing the intrinsic radio flux density in [Jy].
        S_radio_f = np.zeros(len(detectable_radio_survey))
        S_radio_f[detectable_radio_survey] = er.flux_density_radio(
            S_radio_bol[detectable_radio_survey],
            spectral_index[detectable_radio_survey],
            f=self.f_central,
        )
        # Computing the intrinsic radio flux density in [Jy] at a frequency of 1.4 GHz to compare with MeerKAT fluxes.
        S_radio_f_1_4_GHz = np.zeros(len(detectable_radio_survey))
        S_radio_f_1_4_GHz[detectable_radio_survey] = er.flux_density_radio(
            S_radio_bol[detectable_radio_survey],
            spectral_index[detectable_radio_survey],
            f=1.429e9,
        )

        # Compute the effective pulse width in [s] at the survey's central frequency and at 1.4 GHz.
        w_eff = np.zeros(len(detectable_radio_survey))
        w_eff[detectable_radio_survey] = effective_pulse_width(
            w_int_s[detectable_radio_survey],
            DM[detectable_radio_survey],
            self.channel_width,
            self.f_central,
            self.t_samp,
            tau_sc[detectable_radio_survey],
        )
        w_eff_1_4_GHz = np.zeros(len(detectable_radio_survey))
        w_eff_1_4_GHz[detectable_radio_survey] = effective_pulse_width(
            w_int_s[detectable_radio_survey],
            DM[detectable_radio_survey],
            self.channel_width,
            1.429e9,
            self.t_samp,
            tau_sc[detectable_radio_survey],
        )

        # Compute the observed radio flux in [Jy] at the survey's central frequency and at 1.4 GHz.
        S_radio_obs = np.zeros(len(detectable_radio_survey))
        S_radio_obs[detectable_radio_survey] = flux_radio_obs(
            S_radio_f[detectable_radio_survey],
            w_int_s[detectable_radio_survey],
            w_eff[detectable_radio_survey],
        )
        S_radio_obs_1_4_GHz = np.zeros(len(detectable_radio_survey))
        S_radio_obs_1_4_GHz[detectable_radio_survey] = flux_radio_obs(
            S_radio_f_1_4_GHz[detectable_radio_survey],
            w_int_s[detectable_radio_survey],
            w_eff_1_4_GHz[detectable_radio_survey],
        )

        # Compute the period-averaged flux in [Jy] at the survey's central frequency and at 1.4 GHz.
        S_radio_obs_mean = flux_radio_obs_period_average(S_radio_obs, P, w_eff)
        S_radio_obs_mean_1_4GHz = flux_radio_obs_period_average(
            S_radio_obs_1_4_GHz, P, w_eff_1_4_GHz
        )

        detected_radio = np.zeros(len(detectable_radio_survey), dtype=bool)

        detected_radio[detectable_radio_survey] = self.simulate_detection(
            S_radio_obs_mean[detectable_radio_survey],
            l_gal[detectable_radio_survey],
            b_gal[detectable_radio_survey],
            DEC[detectable_radio_survey],
            w_eff[detectable_radio_survey],
            P[detectable_radio_survey],
        )

        return (
            detected_radio,
            w_eff,
            S_radio_obs_mean,
            S_radio_obs_mean_1_4GHz,
        )

__init__(parameters_path)

Radio survey initialization. The parameters for the survey are imported from a JSON file.

Parameters:

Name	Type	Description	Default
parameters_path	str	Path to the survey_parameter.json file containing the parameters of the radio survey. The file must include: deg_factor (float): Degradation factor. G0 (float): Gain at the beam center [KJy ^ (-1)]. t_obs (float): Integration time [s]. t_samp (float): Sampling time [s]. T_sys (float): System temperature [K]. nu_central (float): Central frequency of the bandwidth [Hz]. BW (float): Frequency bandwidth [Hz]. channel_width (float): Width of a single frequency channel [Hz]. n_pol (float): Number of polarizations. FWHM (float): FWHM of the beam[arcmin]. SNR_th (float): Threshold signal-to-noise ratio. RA_range (np.ndarray): Range of the sky covered by the survey in RA [deg]. DEC_range(np.ndarray): Range of the sky covered by the survey in DEC [deg]. l_range(np.ndarray): Range of the sky covered by the survey in Galactic longitude l[deg]. b_range_abs(np.ndarray): Absolute value of the range of the sky covered by the survey in Galactic latitude b [deg]. aperture_config (bool): If True, use the aperture array configuration for the pulsar detection.	required

Source code in mlpoppyns/simulator/multiband_surveys/survey_radio.py

def __init__(self, parameters_path: str) -> None:
    """
    Radio survey initialization.
    The parameters for the survey are imported from a JSON file.

    Args:
        parameters_path (str): Path to the survey_parameter.json file containing
            the parameters of the radio survey. The file must include:

            - deg_factor (float): Degradation factor.
            - G0 (float): Gain at the beam center [KJy ^ (-1)].
            - t_obs (float): Integration time [s].
            - t_samp (float): Sampling time [s].
            - T_sys (float): System temperature [K].
            - nu_central (float): Central frequency of the bandwidth [Hz].
            - BW (float): Frequency bandwidth [Hz].
            - channel_width (float): Width of a single frequency channel [Hz].
            - n_pol (float): Number of polarizations.
            - FWHM (float): FWHM of the beam[arcmin].
            - SNR_th (float): Threshold signal-to-noise ratio.
            - RA_range (np.ndarray): Range of the sky covered by the survey in RA [deg].
            - DEC_range(np.ndarray): Range of the sky covered by the survey in DEC [deg].
            - l_range(np.ndarray): Range of the sky covered by the survey in Galactic longitude l[deg].
            - b_range_abs(np.ndarray): Absolute value of the range of the sky covered
                by the survey in Galactic latitude b [deg].
            - aperture_config (bool): If True, use the aperture array configuration for the pulsar detection.
    """

    # Load parameters from JSON file.
    with open(parameters_path) as read_file:
        self.parameters = json.load(read_file)

    # Save the parameters.
    self.deg_factor = self.parameters["deg_factor"]
    self.G0 = self.parameters["G0"]
    self.t_obs = self.parameters["t_obs"]
    self.t_samp = self.parameters["t_samp"]
    self.T_sys = self.parameters["T_sys"]
    self.f_central = self.parameters["f_central"]
    self.BW = self.parameters["BW"]
    self.channel_width = self.parameters["channel_width"]
    self.n_pol = self.parameters["n_pol"]
    self.FWHM = self.parameters["FWHM"]
    self.SNR_th = self.parameters["SNR_th"]
    self.RA_range = self.parameters["RA_range"]
    self.DEC_range = self.parameters["DEC_range"]
    self.l_range = self.parameters["l_range"]
    self.b_range_abs = self.parameters["b_range_abs"]
    self.aperture_config = self.parameters["aperture_config"]
    self.FFT_search = self.parameters["FFT_search"]
    self.name = self.parameters["name"]

aperture_array_factor(DEC)

Computes the aperture array sensitivity correction factor as a function of declination for each pulsar. Args: DEC (np.ndarray): Declination in [deg] defined between [-90, 90] deg in ICRS frame. Returns: (np.ndarray): Correction factor emulating the sensitivity of an aperture array for different declinations.

Source code in mlpoppyns/simulator/multiband_surveys/survey_radio.py

def aperture_array_factor(self, DEC: np.ndarray) -> np.ndarray:
    """
    Computes the aperture array sensitivity correction factor as a function of declination for each pulsar.
    Args:
        DEC (np.ndarray): Declination in [deg] defined between [-90, 90] deg in ICRS frame.
    Returns:
        (np.ndarray): Correction factor emulating the sensitivity of an aperture array for different declinations.
    """

    # Compute the pulsars' offset angles from zenith and convert them to [rad].
    offset_from_zenith = (
        DEC - (self.DEC_range[0] + self.DEC_range[1]) / 2.0
    ) * const.DEG_TO_RAD

    aa_factor = np.cos(offset_from_zenith)

    return aa_factor

detected_radio_population(w_int_s, DM, P, intercepted_radio, coverage, l_gal, b_gal, DEC, S_radio_bol, spectral_index, tau_sc)

Compute the pulsars detected by each survey.

Parameters:

Name	Type	Description	Default
w_int_s	ndarray	Intrinsic pulse widths in [s]	required
DM	ndarray	Dispersion measure in [pc cm^-3].	required
P	ndarray	Array of spin periods of the pulsars in [s].	required
intercepted_radio	ndarray	(np.ndarray) Array of boolean variables where True values represent stars whose beam crosses our line of sight.	required
coverage	ndarray	Array of boolean variables indicating the pulsars within the sky coverage.	required
l_gal	ndarray	Galactic longitude in [deg] defined between [-180, 180] deg.	required
b_gal	ndarray	Galactic latitude in [deg] defined between [-90, 90] deg.	required
DEC	ndarray	Declination in [deg] defined between [-90, 90] deg in ICRS frame.	required
S_radio_bol	ndarray	Pulsar bolometric radio flux in [erg s^(-1) cm^(-2)].	required
spectral_index	ndarray	Spectral indexes.	required
tau_sc	ndarray	Scattering timescale in [s].	required

Returns:

Type	Description
Tuple[ndarray, ndarray, ndarray, ndarray]	Tuple consisting of the following arrays: Boolean mask to select the pulsars detected by the survey. Effective pulse width in [s]. Period-averaged fluxes at the central frequency of the survey in [Jy] Period-averaged fluxes at the central frequency of 1.429 GHz in [Jy].

Source code in mlpoppyns/simulator/multiband_surveys/survey_radio.py

def detected_radio_population(
    self,
    w_int_s: np.ndarray,
    DM: np.ndarray,
    P: np.ndarray,
    intercepted_radio: np.ndarray,
    coverage: np.ndarray,
    l_gal: np.ndarray,
    b_gal: np.ndarray,
    DEC: np.ndarray,
    S_radio_bol: np.ndarray,
    spectral_index: np.ndarray,
    tau_sc: np.ndarray,
) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
    """
    Compute the pulsars detected by each survey.

    Args:
        w_int_s (np.ndarray): Intrinsic pulse widths in [s]
        DM (np.ndarray): Dispersion measure in [pc cm^-3].
        P (np.ndarray): Array of spin periods of the pulsars in [s].
        intercepted_radio: (np.ndarray) Array of boolean variables where True values represent stars whose
            beam crosses our line of sight.
        coverage (np.ndarray): Array of boolean variables indicating the pulsars within the sky coverage.
        l_gal (np.ndarray): Galactic longitude in [deg] defined between [-180, 180] deg.
        b_gal (np.ndarray): Galactic latitude in [deg] defined between [-90, 90] deg.
        DEC (np.ndarray): Declination in [deg] defined between [-90, 90] deg in ICRS frame.
        S_radio_bol (np.ndarray): Pulsar bolometric radio flux in [erg s^(-1) cm^(-2)].
        spectral_index (np.ndarray): Spectral indexes.
        tau_sc (np.ndarray): Scattering timescale in [s].

    Returns:
        (Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]): Tuple consisting of the following arrays:

            - Boolean mask to select the pulsars detected by the survey.
            - Effective pulse width in [s].
            - Period-averaged fluxes at the central frequency of the survey in [Jy]
            - Period-averaged fluxes at the central frequency of 1.429 GHz in [Jy].
    """

    detectable_radio_survey = intercepted_radio & coverage

    # Computing the intrinsic radio flux density in [Jy].
    S_radio_f = np.zeros(len(detectable_radio_survey))
    S_radio_f[detectable_radio_survey] = er.flux_density_radio(
        S_radio_bol[detectable_radio_survey],
        spectral_index[detectable_radio_survey],
        f=self.f_central,
    )
    # Computing the intrinsic radio flux density in [Jy] at a frequency of 1.4 GHz to compare with MeerKAT fluxes.
    S_radio_f_1_4_GHz = np.zeros(len(detectable_radio_survey))
    S_radio_f_1_4_GHz[detectable_radio_survey] = er.flux_density_radio(
        S_radio_bol[detectable_radio_survey],
        spectral_index[detectable_radio_survey],
        f=1.429e9,
    )

    # Compute the effective pulse width in [s] at the survey's central frequency and at 1.4 GHz.
    w_eff = np.zeros(len(detectable_radio_survey))
    w_eff[detectable_radio_survey] = effective_pulse_width(
        w_int_s[detectable_radio_survey],
        DM[detectable_radio_survey],
        self.channel_width,
        self.f_central,
        self.t_samp,
        tau_sc[detectable_radio_survey],
    )
    w_eff_1_4_GHz = np.zeros(len(detectable_radio_survey))
    w_eff_1_4_GHz[detectable_radio_survey] = effective_pulse_width(
        w_int_s[detectable_radio_survey],
        DM[detectable_radio_survey],
        self.channel_width,
        1.429e9,
        self.t_samp,
        tau_sc[detectable_radio_survey],
    )

    # Compute the observed radio flux in [Jy] at the survey's central frequency and at 1.4 GHz.
    S_radio_obs = np.zeros(len(detectable_radio_survey))
    S_radio_obs[detectable_radio_survey] = flux_radio_obs(
        S_radio_f[detectable_radio_survey],
        w_int_s[detectable_radio_survey],
        w_eff[detectable_radio_survey],
    )
    S_radio_obs_1_4_GHz = np.zeros(len(detectable_radio_survey))
    S_radio_obs_1_4_GHz[detectable_radio_survey] = flux_radio_obs(
        S_radio_f_1_4_GHz[detectable_radio_survey],
        w_int_s[detectable_radio_survey],
        w_eff_1_4_GHz[detectable_radio_survey],
    )

    # Compute the period-averaged flux in [Jy] at the survey's central frequency and at 1.4 GHz.
    S_radio_obs_mean = flux_radio_obs_period_average(S_radio_obs, P, w_eff)
    S_radio_obs_mean_1_4GHz = flux_radio_obs_period_average(
        S_radio_obs_1_4_GHz, P, w_eff_1_4_GHz
    )

    detected_radio = np.zeros(len(detectable_radio_survey), dtype=bool)

    detected_radio[detectable_radio_survey] = self.simulate_detection(
        S_radio_obs_mean[detectable_radio_survey],
        l_gal[detectable_radio_survey],
        b_gal[detectable_radio_survey],
        DEC[detectable_radio_survey],
        w_eff[detectable_radio_survey],
        P[detectable_radio_survey],
    )

    return (
        detected_radio,
        w_eff,
        S_radio_obs_mean,
        S_radio_obs_mean_1_4GHz,
    )

detection_offset(n_detection)

Generating a random offset with respect to the beam center for the detections. A Gaussian beam pattern is assumed to account for sensitivity decay for off-center detections. See Lorimer et al. (1993) and the paragraph following eq. (29) in Bates et al. (2014).

Parameters:

Name	Type	Description	Default
n_detection	int	Number of detections to simulate.	required

Returns:

Type	Description
ndarray	Square of the offset from the beam center for each detection in [arcmin^2] .

Source code in mlpoppyns/simulator/multiband_surveys/survey_radio.py

def detection_offset(self, n_detection: int) -> np.ndarray:
    """
    Generating a random offset with respect to the beam center for the detections.
    A Gaussian beam pattern is assumed to account for sensitivity decay for off-center detections.
    See Lorimer et al. (1993) and the paragraph following eq. (29) in Bates et al. (2014).

    Args:
        n_detection (int): Number of detections to simulate.

    Returns:
        (np.ndarray): Square of the offset from the beam center for each detection in [arcmin^2] .
    """
    offset2 = np.random.uniform(0.0, self.FWHM**2 / 4.0, n_detection)
    return offset2

fft_search_efficiency(duty_cycle)

Compute the efficiency factor from Morello et al. (2020) (see eq. 44) to account for incoherent FFT search.

Parameters:

Name	Type	Description	Default
duty_cycle	ndarray	duty cycle of pulsars.	required

Returns: (np.ndarray): Correction factor emulating the sensitivity of an incoherent FFT search.

Source code in mlpoppyns/simulator/multiband_surveys/survey_radio.py

def fft_search_efficiency(self, duty_cycle: np.ndarray) -> np.ndarray:
    """
    Compute the efficiency factor from Morello et al. (2020) (see eq. 44) to account for incoherent FFT search.

    Args:
        duty_cycle (np.ndarray): duty cycle of pulsars.
    Returns:
        (np.ndarray): Correction factor emulating the sensitivity of an incoherent FFT search.
    """

    epsilon = (1.0 + 0.0473 * duty_cycle ** (-0.627)) ** (-1)

    return epsilon

gain_gaussian_beam(offset2)

This method simulates the gain pattern of a receiver. A Gaussian beam pattern is assumed to account for sensitivity decay for off-center detections. See Lorimer et al. (1993) and the paragraph following eq. (29) in Bates et al. (2014).

Parameters:

Name	Type	Description	Default
offset2	ndarray	Squared offset from the beam center in [arcmin^2].	required

Returns:

Type	Description
ndarray	Gain of the telescope for the given offset in [K Jy^(-1)].

Source code in mlpoppyns/simulator/multiband_surveys/survey_radio.py

def gain_gaussian_beam(self, offset2: np.ndarray) -> np.ndarray:
    """
    This method simulates the gain pattern of a receiver.
    A Gaussian beam pattern is assumed to account for sensitivity decay for off-center detections.
    See Lorimer et al. (1993) and the paragraph following eq. (29) in Bates et al. (2014).

    Args:
        offset2 (np.ndarray): Squared offset from the beam center in [arcmin^2].

    Returns:
        (np.ndarray): Gain of the telescope for the given offset in [K Jy^(-1)].
    """
    G = self.G0 * np.exp(-2.77 * offset2 / self.FWHM**2)

    return G

radiometer_equation(S_radio_obs_mean, G, w_eff, P, T_sky)

Radiometer equation used to compute the signal-to-noise ratio of each pulsar given the observed period-averaged radio flux at a given frequency, the effective pulse width, the spin period and the survey parameters (see eq. (A1.22) in Lorimer & Kramer 2005). We are assuming a square pulse shape for simplicity with height equal to the observed flux and width equal to the effective width.

Parameters:

Name	Type	Description	Default
S_radio_obs_mean	ndarray	Observed period-averaged radio flux density in [Jy].	required
G	ndarray	Gain of the telescope for the given detection in [K Jy^(-1)].	required
w_eff	ndarray	Effective pulse width in [s].	required
P	ndarray	Spin period in [s].	required
T_sky	ndarray	Sky temperature for every detection in [K].	required

Returns:

Type	Description
ndarray	Signal-to-noise ratio of the detection.

Source code in mlpoppyns/simulator/multiband_surveys/survey_radio.py

def radiometer_equation(
    self,
    S_radio_obs_mean: np.ndarray,
    G: np.ndarray,
    w_eff: np.ndarray,
    P: np.ndarray,
    T_sky: np.ndarray,
) -> np.ndarray:
    """
    Radiometer equation used to compute the signal-to-noise ratio of each pulsar given the observed
    period-averaged radio flux at a given frequency, the effective pulse width, the spin period and
    the survey parameters (see eq. (A1.22) in Lorimer & Kramer 2005). We are assuming a square pulse
    shape for simplicity with height equal to the observed flux and width equal to the effective width.

    Args:
        S_radio_obs_mean (np.ndarray): Observed period-averaged radio flux density in [Jy].
        G (np.ndarray): Gain of the telescope for the given detection in [K Jy^(-1)].
        w_eff (np.ndarray): Effective pulse width in [s].
        P (np.ndarray): Spin period in [s].
        T_sky (np.ndarray): Sky temperature for every detection in [K].

    Returns:
        (np.ndarray): Signal-to-noise ratio of the detection.
    """
    SNR = np.zeros(len(S_radio_obs_mean))

    # If the effective pulse width is larger than the spin period,
    # emission is continuous and the neutron star cannot be detected as a pulsar.
    cond = w_eff < P

    # Compute the SNR of each detection using the radiometer equation.
    SNR[cond] = (
        S_radio_obs_mean[cond]
        * G[cond]
        * np.sqrt(self.n_pol * self.t_obs * self.BW)
        * np.sqrt((P[cond] - w_eff[cond]) / w_eff[cond])
        / (self.deg_factor * (self.T_sys + T_sky[cond]))
    )

    return SNR

simulate_detection(S_radio_obs_mean, l_gal, b_gal, DEC, w_eff, P)

Simulate a detection: If the measured SNR surpasses the threshold SNR_th of the survey then the pulsar is detected.

Parameters:

Name	Type	Description	Default
S_radio_obs_mean	ndarray	Observed period-averaged radio flux density in [Jy].	required
l_gal	ndarray	Galactic longitude in [deg] defined between [-180, 180] deg.	required
b_gal	ndarray	Galactic latitude in [deg] defined between [-90, 90] deg.	required
DEC	ndarray	Declination in [deg] defined between [-90, 90] deg in ICRS frame.	required
w_eff	ndarray	Effective pulse width in [s].	required
P	ndarray	Spin period in [s].	required

Returns:

Type	Description
ndarray	Array of boolean variables: True if the pulsar is detected, False if not.

Source code in mlpoppyns/simulator/multiband_surveys/survey_radio.py

def simulate_detection(
    self,
    S_radio_obs_mean: np.ndarray,
    l_gal: np.ndarray,
    b_gal: np.ndarray,
    DEC: np.ndarray,
    w_eff: np.ndarray,
    P: np.ndarray,
) -> np.ndarray:
    """
    Simulate a detection: If the measured SNR surpasses the threshold SNR_th of the survey
    then the pulsar is detected.

    Args:
        S_radio_obs_mean (np.ndarray): Observed period-averaged radio flux density in [Jy].
        l_gal (np.ndarray): Galactic longitude in [deg] defined between [-180, 180] deg.
        b_gal (np.ndarray): Galactic latitude in [deg] defined between [-90, 90] deg.
        DEC (np.ndarray): Declination in [deg] defined between [-90, 90] deg in ICRS frame.
        w_eff (np.ndarray): Effective pulse width in [s].
        P (np.ndarray): Spin period in [s].

    Returns:
        (np.ndarray): Array of boolean variables: True if the pulsar is detected, False if not.
    """

    # Store the total number of sources.
    n = len(S_radio_obs_mean)

    # Draw a random offset from the telescope beam center.
    offset2 = self.detection_offset(n)
    # Compute the gain corresponding to the offset detections.
    G = self.gain_gaussian_beam(offset2)

    # Compute the sky temperature in the coordinates of each detection at the central frequency of the survey.
    # We choose here to use the refined map from Remazeilles et al. (2015).
    T_sky = sky_temperature_H81refined(l_gal, b_gal, self.f_central)
    SNR_detection = self.radiometer_equation(
        S_radio_obs_mean, G, w_eff, P, T_sky
    )

    if self.aperture_config:
        aa_factor = self.aperture_array_factor(DEC)
        SNR_detection = SNR_detection * aa_factor

    if self.FFT_search:
        # Apply efficiency factor from Morello et al. (2020) (see eq. 44) to account for incoherent FFT search.
        duty_cycle = w_eff / P
        epsilon = self.fft_search_efficiency(duty_cycle)
        SNR_detection = SNR_detection * epsilon

    detected = SNR_detection > self.SNR_th

    return detected

sky_coverage(RA, DEC, l_gal, b_gal)

Determine which neutron stars are in the sky region covered by the survey.

Parameters:

Name	Type	Description	Default
RA	ndarray	Right ascension in [deg] defined between [0, 360] deg in ICRS frame.	required
DEC	ndarray	Declination in [deg] defined between [-90, 90] deg in ICRS frame.	required
l_gal	ndarray	Galactic longitude in [deg] defined between [-180, 180] deg.	required
b_gal	ndarray	Galactic latitude in [deg] defined between [-90, 90] deg.	required

Returns:

Type	Description
ndarray	Array of boolean variables: True if the pulsar is in the covered sky region, False if not.

Source code in mlpoppyns/simulator/multiband_surveys/survey_radio.py

def sky_coverage(
    self,
    RA: np.ndarray,
    DEC: np.ndarray,
    l_gal: np.ndarray,
    b_gal: np.ndarray,
) -> np.ndarray:
    """
    Determine which neutron stars are in the sky region covered by the survey.

    Args:
        RA (np.ndarray): Right ascension in [deg] defined between [0, 360] deg in ICRS frame.
        DEC (np.ndarray): Declination in [deg] defined between [-90, 90] deg in ICRS frame.
        l_gal (np.ndarray): Galactic longitude in [deg] defined between [-180, 180] deg.
        b_gal (np.ndarray): Galactic latitude in [deg] defined between [-90, 90] deg.

    Returns:
        (np.ndarray): Array of boolean variables: True if the pulsar is in the covered sky region, False if not.
    """
    if self.name == "HTRU high":
        coverage = (
            (RA > self.RA_range[0])
            & (RA < self.RA_range[1])
            & (DEC > self.DEC_range[0])
            & (DEC < self.DEC_range[1])
            & (
                (
                    (
                        (l_gal > self.l_range[0][0])
                        & (l_gal < self.l_range[0][1])
                    )
                    | (
                        (l_gal > self.l_range[1][0])
                        & (l_gal < self.l_range[1][1])
                    )
                )
                | (
                    (np.abs(b_gal) > self.b_range_abs[0])
                    & (np.abs(b_gal) < self.b_range_abs[1])
                )
            )
        )

    else:
        coverage = (
            (RA > self.RA_range[0])
            & (RA < self.RA_range[1])
            & (DEC > self.DEC_range[0])
            & (DEC < self.DEC_range[1])
            & (l_gal > self.l_range[0])
            & (l_gal < self.l_range[1])
            & (np.abs(b_gal) > self.b_range_abs[0])
            & (np.abs(b_gal) < self.b_range_abs[1])
        )

    return coverage

effective_pulse_width(w_int, DM, channel_width, f, t_samp, tau_sc)

Measured effective pulse width which is smeared out by the inter-channel dispersion, the scattering with the interstellar medium and by the instrumental sampling time (see eq. (2) in Cordes & McLaughlin 2003).

Parameters:

Name	Type	Description	Default
w_int	ndarray	Intrinsic pulse width in [s].	required
DM	ndarray	Dispersion measure in [pc cm^-3].	required
channel_width	float	Width in frequency of a single frequency channel of the receiver in [Hz].	required
f	float	Central frequency at which the observation is performed [Hz].	required
t_samp	float	Sampling time for the radio survey [s].	required
tau_sc	ndarray	Scattering timescales in [s].	required

Returns: (np.ndarray): Measured effective pulse width in [s].

Source code in mlpoppyns/simulator/multiband_surveys/survey_radio.py

def effective_pulse_width(
    w_int: np.ndarray,
    DM: np.ndarray,
    channel_width: float,
    f: float,
    t_samp: float,
    tau_sc: np.ndarray,
) -> np.ndarray:
    """
    Measured effective pulse width which is smeared out by the inter-channel dispersion, the scattering
    with the interstellar medium and by the instrumental sampling time (see eq. (2) in Cordes & McLaughlin 2003).

    Args:
        w_int (np.ndarray): Intrinsic pulse width in [s].
        DM (np.ndarray): Dispersion measure in [pc cm^-3].
        channel_width (float): Width in frequency of a single frequency channel of the receiver in [Hz].
        f (float): Central frequency at which the observation is performed [Hz].
        t_samp (float): Sampling time for the radio survey [s].
        tau_sc (np.ndarray): Scattering timescales in [s].
    Returns:
        (np.ndarray): Measured effective pulse width in [s].
    """

    tau_DM = smearing_in_channel(DM, channel_width, f)

    # Convert the scattering time to a given observation frequency f assuming a Kolmogorov spectrum.
    tau_sc_f = edm.compute_tau_sc_f(tau_sc, f)

    w_eff = np.sqrt(w_int**2 + tau_sc_f**2 + tau_DM**2 + t_samp**2)

    return w_eff

flux_radio_obs(S_radio_f, w_int, w_eff)

Compute the flux density received by the telescope after taking into account that the pulse has been broadened by the propagation in the interstellar medium. We assume that the total fluence = S_radio_f x w_int is conserved as the pulse propagates in the interstellar medium. Since the pulse is broadened as it propagates, the flux received on Earth is given by S_radio_obs = fluence / w_eff, therefore S_radio_obs < S_radio_f. Also, since in our simulation we are assuming a simple squared pulse shape, the flux density computed here is equal to the peak flux density.

Parameters:

Name	Type	Description	Default
S_radio_f	ndarray	Intrinsic pulsar radio flux in [Jy].	required
w_int	ndarray	Intrinsic pulse width in [rad].	required
w_eff	ndarray	Effective pulse width in [rad].	required

Returns:

Type	Description
ndarray	Observed pulsar radio flux in [Jy].

Source code in mlpoppyns/simulator/multiband_surveys/survey_radio.py

def flux_radio_obs(
    S_radio_f: np.ndarray,
    w_int: np.ndarray,
    w_eff: np.ndarray,
) -> np.ndarray:
    """
    Compute the flux density received by the telescope after taking into account that the pulse has been
    broadened by the propagation in the interstellar medium. We assume that the total fluence = S_radio_f x w_int
    is conserved as the pulse propagates in the interstellar medium. Since the pulse is broadened as it propagates,
    the flux received on Earth is given by S_radio_obs = fluence / w_eff, therefore S_radio_obs < S_radio_f.
    Also, since in our simulation we are assuming a simple squared pulse shape, the flux density computed here is equal
    to the peak flux density.

    Args:
        S_radio_f (np.ndarray): Intrinsic pulsar radio flux in [Jy].
        w_int (np.ndarray): Intrinsic pulse width in [rad].
        w_eff (np.ndarray): Effective pulse width in [rad].

    Returns:
        (np.ndarray): Observed pulsar radio flux in [Jy].
    """

    # Compute the total fluence.
    fluence_f = S_radio_f * w_int

    # Compute the observed radio flux in [Jy].
    S_radio_f_obs = fluence_f / w_eff

    return S_radio_f_obs

flux_radio_obs_period_average(S_radio_f_obs, P, w_eff)

Compute the mean flux density received by the telescope averaged over a spin period. We are assuming a simple squared pulse shape.

Parameters:

Name	Type	Description	Default
S_radio_f_obs	ndarray	Observed pulsar radio flux in [Jy].	required
P	ndarray	Spin period of the pulsar in [s].	required
w_eff	ndarray	Effective pulse width in [s].	required

Returns:

Type	Description
ndarray	Observed pulsar radio flux averaged over a period in [Jy].

Source code in mlpoppyns/simulator/multiband_surveys/survey_radio.py

def flux_radio_obs_period_average(
    S_radio_f_obs: np.ndarray,
    P: np.ndarray,
    w_eff: np.ndarray,
) -> np.ndarray:
    """
    Compute the mean flux density received by the telescope averaged over a spin period.
    We are assuming a simple squared pulse shape.

    Args:
        S_radio_f_obs (np.ndarray): Observed pulsar radio flux in [Jy].
        P (np.ndarray): Spin period of the pulsar in [s].
        w_eff (np.ndarray): Effective pulse width in [s].

    Returns:
        (np.ndarray): Observed pulsar radio flux averaged over a period in [Jy].
    """

    # Compute the observed radio flux in [Jy].
    S_radio_f_obs_mean = S_radio_f_obs * w_eff / P

    return S_radio_f_obs_mean

sky_temperature_H81(l_gal, b_gal, f)

Sky temperature as a function of Galactic longitude and latitude (l, b) and frequency f. This function implements the sky temperature map at 408 MHz from Haslam et al. (1981), downloadable here: https://lambda.gsfc.nasa.gov/product/foreground/haslam_408.cfm.

Parameters:

Name	Type	Description	Default
l_gal	ndarray	Galactic longitude in [deg] defined between [-180, 180] deg.	required
b_gal	ndarray	Galactic latitude in [deg] defined between [-90, 90] deg.	required
f	float	Central frequency at which the observation is performed [Hz].	required

Returns:

Type	Description
ndarray	Measured sky temperature in [K] as a function of the Galactic coordinates at frequency f.

Source code in mlpoppyns/simulator/multiband_surveys/survey_radio.py

def sky_temperature_H81(
    l_gal: np.ndarray, b_gal: np.ndarray, f: float
) -> np.ndarray:
    """
    Sky temperature as a function of Galactic longitude and latitude (l, b) and frequency f.
    This function implements the sky temperature map at 408 MHz from Haslam et al. (1981), downloadable here:
    https://lambda.gsfc.nasa.gov/product/foreground/haslam_408.cfm.

    Args:
        l_gal (np.ndarray): Galactic longitude in [deg] defined between [-180, 180] deg.
        b_gal (np.ndarray): Galactic latitude in [deg] defined between [-90, 90] deg.
        f (float): Central frequency at which the observation is performed [Hz].

    Returns:
        (np.ndarray): Measured sky temperature in [K] as a function of the Galactic coordinates at frequency f.
    """

    # Read the sky temperature map.
    file = pathlib.Path().joinpath(
        cfg["path_to_software"],
        "mlpoppyns/simulator/multiband_surveys/Tsky_map_haslam81.fits",
    )
    hdulist = fits.open(file)
    hdu = hdulist["TEMPERATURE"]
    data = hdu.data

    # Convert coordinates into astropy coordinates object.
    coord = SkyCoord(l_gal, b_gal, frame="galactic", unit="deg")

    # Convert sky coordinates into pixel coordinates and extract the temperatures.
    wcs = WCS(hdu.header)
    x_pixel, y_pixel = wcs.world_to_pixel(coord)
    x_pixel = x_pixel.astype(int)
    y_pixel = y_pixel.astype(int)
    T_sky_400 = data[y_pixel, x_pixel]

    # Rescale to the wanted frequency assuming a sky temperature spectral index of -2.6
    # (see Lawson et al. 1987, Johnston et al. 1992).
    T_sky_f = T_sky_400 * (408.0e6 / f) ** 2.6

    return T_sky_f

sky_temperature_H81refined(l_gal, b_gal, f)

Sky temperature as a function of Galactic longitude and latitude (l, b) and frequency f. This function implements the sky temperature map at 408 MHz from Remazeilles et al. (2015) which is a refinement of the map from Haslam et al. (1981). The map is downloadable here: https://lambda.gsfc.nasa.gov/product/foreground/fg_2014_haslam_408_get.html.

Parameters:

Name	Type	Description	Default
l_gal	ndarray	Galactic longitude in [deg] defined between [-180, 180] deg.	required
b_gal	ndarray	Galactic latitude in [deg] defined between [-90, 90] deg.	required
f	float	Central frequency at which the observation is performed [Hz].	required

Returns:

Type	Description
ndarray	Measured sky temperature in [K] as a function of the Galactic coordinates at frequency f.

Source code in mlpoppyns/simulator/multiband_surveys/survey_radio.py

def sky_temperature_H81refined(
    l_gal: np.ndarray, b_gal: np.ndarray, f: float
) -> np.ndarray:
    """
    Sky temperature as a function of Galactic longitude and latitude (l, b) and frequency f.
    This function implements the sky temperature map at 408 MHz from Remazeilles et al. (2015) which is a
    refinement of the map from Haslam et al. (1981).
    The map is downloadable here:
    https://lambda.gsfc.nasa.gov/product/foreground/fg_2014_haslam_408_get.html.

    Args:
        l_gal (np.ndarray): Galactic longitude in [deg] defined between [-180, 180] deg.
        b_gal (np.ndarray): Galactic latitude in [deg] defined between [-90, 90] deg.
        f (float): Central frequency at which the observation is performed [Hz].

    Returns:
        (np.ndarray): Measured sky temperature in [K] as a function of the Galactic coordinates at frequency f.
    """

    # Read the sky temperature map.
    file = pathlib.Path().joinpath(
        cfg["path_to_software"],
        "mlpoppyns/simulator/multiband_surveys/Tsky_map_haslam81_refined.fits",
    )

    T_sky_map = hp.read_map(file, dtype=np.float64)

    # Convert coordinates into astropy coordinates object.
    coord = SkyCoord(l_gal, b_gal, frame="galactic", unit="deg")

    # Convert sky coordinates into pixel coordinates and extract the temperatures.
    l_g = coord.l.degree
    b_g = coord.b.degree
    vec = np.array(hp.rotator.dir2vec(l_g, b_g, lonlat=True))
    n_side = 512
    pix = hp.pixelfunc.vec2pix(n_side, vec[0], vec[1], vec[2], nest=False)

    T_sky_400 = T_sky_map[pix]

    # Rescale to the wanted frequency assuming a sky temperature spectral index of -2.6
    # (see Lawson et al. 1987, Johnston et al. 1992).
    T_sky_f = np.atleast_1d(T_sky_400 * (408.0e6 / f) ** 2.6)

    return T_sky_f

sky_temperature_approx(l_gal, b_gal, f)

Sky temperature as a function of Galactic longitude and latitude (l, b) and frequency f. This function uses an empirical fit from Narayan (1987) and rescale to the given frequency using a relation from Johnston et al. (1992) (see also eq. (5) in Yusifov & Küçük 2004).

Parameters:

Name	Type	Description	Default
l_gal	ndarray	Galactic longitude in [deg] defined between [-180, 180] deg.	required
b_gal	ndarray	Galactic latitude in [deg] defined between [-90, 90] deg.	required
f	float	Central frequency at which the observation is performed [Hz].	required

Returns:

Type	Description
ndarray	Measured sky temperature in [K] as a function of the Galactic coordinates at frequency f.

Source code in mlpoppyns/simulator/multiband_surveys/survey_radio.py

def sky_temperature_approx(
    l_gal: np.ndarray, b_gal: np.ndarray, f: float
) -> np.ndarray:
    """
    Sky temperature as a function of Galactic longitude and latitude (l, b) and frequency f.
    This function uses an empirical fit from Narayan (1987) and rescale to the given frequency using
    a relation from Johnston et al. (1992) (see also eq. (5) in Yusifov & Küçük 2004).

    Args:
        l_gal (np.ndarray): Galactic longitude in [deg] defined between [-180, 180] deg.
        b_gal (np.ndarray): Galactic latitude in [deg] defined between [-90, 90] deg.
        f (float): Central frequency at which the observation is performed [Hz].

    Returns:
        (np.ndarray): Measured sky temperature in [K] as a function of the Galactic coordinates at frequency f.
    """

    # Sky temperature at 408 MHz from Narayan (1987).
    T_sky_400 = 25.0 + 275.0 / (
        (1.0 + (l_gal / 42.0) ** 2) * (1.0 + (b_gal / 3.0) ** 2)
    )

    # Rescale to the wanted frequency assuming a sky temperature spectral index of -2.6
    # (see Lawson et al. 1987, Johnston et al. 1992).
    T_sky_f = T_sky_400 * (408.0e6 / f) ** 2.6

    return T_sky_f

smearing_in_channel(DM, channel_width, nu)

Dispersive smearing inside a single frequency channel in [s] evaluated for the central frequency of the survey. See eq. (27) in Bates et al. (2014) and appendix A2.4 of Handbook of pulsar astronomy by Lorimer and Kramer (2004).

Parameters:

Name	Type	Description	Default
DM	ndarray	Dispersion measure in [pc cm^-3].	required
channel_width	float	Width in frequency of a single frequency channel of the receiver in [Hz].	required
nu	float	Central frequency at which the observation is performed [Hz].	required

Returns:

Type	Description
ndarray	Intra-channel dispersive smearing in [s].

Source code in mlpoppyns/simulator/multiband_surveys/survey_radio.py

def smearing_in_channel(
    DM: np.ndarray, channel_width: float, nu: float
) -> np.ndarray:
    """
    Dispersive smearing inside a single frequency channel in [s] evaluated for
    the central frequency of the survey. See eq. (27) in Bates et al. (2014)
    and appendix A2.4 of Handbook of pulsar astronomy by Lorimer and Kramer (2004).

    Args:
        DM (np.ndarray): Dispersion measure in [pc cm^-3].
        channel_width (float): Width in frequency of a single frequency channel of the receiver in [Hz].
        nu (float): Central frequency at which the observation is performed [Hz].

    Returns:
        (np.ndarray): Intra-channel dispersive smearing in [s].
    """

    dt = (
        2
        * const.E**2
        / (2 * np.pi * const.M_E * const.C)
        * channel_width
        / nu**3.0
        * DM
        * const.PC_TO_CM
    )

    return dt

mlpoppyns.simulator.multiband_surveys.surveys_wrapper

A wrapper module to include all methods used in the script simulate_population_magrot_det.py to initialize and simulate surveys to detect neutron stars in different wavelengths.

Authors:

Vanessa Graber (graber @ ice.csic.es)
Michele Ronchi (ronchi @ ice.csic.es)
Alberto Garcia-Garcia (garciagarcia @ ice.csic.es)
Celsa Pardo Araujo (pardo @ ice.csic.es)

SurveyData dataclass

A dataclass object to store all the survey data for a simulation.

Attributes:

Name	Type	Description
surveys_cfg	Dict	A dictionary to save all the survey information from config_simulator.
surveys_radio	Dict	A dictionary to store the radio survey class objects.
n_detected_sim	Dict	A dictionary to store the number of stars detected in each survey.
n_detected_complete_sim	Dict	A dictionary to store how many neutron stars are detected in the flux ranges where we assume completeness.
percentage_detected	Dict	A dictionary storing the percentage of detections compared to real detected numbers.
n_created_at_match	Dict	A dictionary storing how many stars we have created to reach the desirable number in each survey.
n_detected_sim_at_match	Dict	A dictionary storing how many stars we have detected when we reach the desirable number in each survey.
batchsize_adjust_flags	Dict	A dictionary storing the flags for adjusting batch size as we approach the target detection number for all surveys.
stop_flags	Dict	A dictionary storing the flags for stopping the simulation as we reach the target detection number for all surveys.
dictionary_detected_radio	Dict	A dictionary storing all the properties of neutron stars detected in the radio surveys.

Source code in mlpoppyns/simulator/multiband_surveys/surveys_wrapper.py

@dataclass
class SurveyData:
    """
    A dataclass object to store all the survey data for a simulation.

    Attributes:
        surveys_cfg (Dict): A dictionary to save all the survey information from config_simulator.
        surveys_radio (Dict): A dictionary to store the radio survey class objects.
        n_detected_sim (Dict): A dictionary to store the number of stars detected in each survey.
        n_detected_complete_sim (Dict): A dictionary to store how many neutron stars are detected in the flux ranges
            where we assume completeness.
        percentage_detected (Dict): A dictionary storing the percentage of detections compared to real detected numbers.
        n_created_at_match (Dict): A dictionary storing how many stars we have created to reach the desirable number in
            each survey.
        n_detected_sim_at_match (Dict): A dictionary storing how many stars we have detected when we reach the desirable
            number in each survey.
        batchsize_adjust_flags (Dict): A dictionary storing the flags for adjusting batch size as we approach the target
            detection number for all surveys.
        stop_flags (Dict): A dictionary storing the flags for stopping the simulation as we reach the target detection
            number for all surveys.
        dictionary_detected_radio (Dict): A dictionary storing all the properties of neutron stars detected in the radio
            surveys.
    """

    surveys_cfg: Dict
    surveys_radio: Dict
    n_detected_sim: Dict
    n_detected_complete_sim: Dict
    percentage_detected: Dict
    n_created_at_match: Dict
    n_detected_sim_at_match: Dict
    batchsize_adjust_flags: Dict
    stop_flags: Dict
    dictionary_detected_radio: Dict

adjust_n_batchsize(survey_data_class)

To speed up the simulation, generate new neutron stars in batches.

The batchsize is adjusted as the synthetic population approaches the observed number of neutron stars in the real surveys. This guarantees a better fine-tuning of the simulated detected numbers.

Parameters:

Name	Type	Description	Default
survey_data_class	SurveyData	The SurveyData dataclass containing the data of all neutron star surveys.	required

Source code in mlpoppyns/simulator/multiband_surveys/surveys_wrapper.py

def adjust_n_batchsize(survey_data_class) -> int:
    """
    To speed up the simulation, generate new neutron stars in batches.

    The batchsize is adjusted as the synthetic population approaches the observed number of neutron stars in the real
    surveys. This guarantees a better fine-tuning of the simulated detected numbers.

    Args:
        survey_data_class (SurveyData): The SurveyData dataclass containing the data of all neutron star surveys.
    """

    surveys_cfg = survey_data_class.surveys_cfg
    n_detected_complete_sim = survey_data_class.n_detected_complete_sim
    percentage_detected = survey_data_class.percentage_detected
    batchsize_adjusting_flags = survey_data_class.batchsize_adjust_flags

    # Evaluate the percentage of neutron stars detected by the simulated
    # surveys with respect to the real surveys and adjust the batch size accordingly.
    for survey in surveys_cfg:
        percentage_detected[survey] = (
            n_detected_complete_sim[survey]
            / surveys_cfg[survey]["detected_real"]
        )

    if (
        all(value > 0.9 for value in percentage_detected.values())
        and not batchsize_adjusting_flags[0.9]
    ):
        batchsize_adjusting_flags[0.9] = True
        n_batchsize = 10000
    elif (
        all(value > 0.95 for value in percentage_detected.values())
        and not batchsize_adjusting_flags[0.95]
    ):
        batchsize_adjusting_flags[0.95] = True
        n_batchsize = 5000

    else:
        n_batchsize = 100000

    return n_batchsize

apply_surveys_coverage_filter(surveys_radio, coverage_dict, pop_dict, idx_remove)

Apply survey coverage criteria to filter a dynamic population dataset based on sky coverage of all surveys and a distance cutoff, and update the indices of entries to be removed.

Parameters:

Name	Type	Description	Default
surveys_radio	dict	A dictionary of radio survey objects, containing the information on the sky coverage.	required
coverage_dict	dict	A dictionary with the sky coverage information for all surveys.	required
pop_dict	dict	A dictionary containing the data of a neutron star population.	required
idx_remove	list	A list of indices of entries to be removed based on the filtering criteria.	required

Returns:

Type	Description
Tuple[dict, list]	A tuple object containing the following attributes: A dictionary containing data for stars that meet the coverage criteria. An updated list of indices of stars that are outside the coverage and should be removed.

Source code in mlpoppyns/simulator/multiband_surveys/surveys_wrapper.py

def apply_surveys_coverage_filter(
    surveys_radio: dict,
    coverage_dict: dict,
    pop_dict: dict,
    idx_remove: list,
) -> Tuple[dict, list]:
    """
    Apply survey coverage criteria to filter a dynamic population dataset based on sky coverage of all surveys and a
    distance cutoff, and update the indices of entries to be removed.

    Args:
        surveys_radio (dict): A dictionary of radio survey objects, containing the information on the sky coverage.
        coverage_dict (dict): A dictionary with the sky coverage information for all surveys.
        pop_dict (dict): A dictionary containing the data of a neutron star population.
        idx_remove (list): A list of indices of entries to be removed based on the filtering criteria.

    Returns:
        (Tuple[dict, list]): A tuple object containing the following attributes:

            - A dictionary containing data for stars that meet the coverage criteria.
            - An updated list of indices of stars that are outside the coverage and should be removed.
    """

    survey_radio_names = list(surveys_radio.keys())

    coverage_tot = coverage_dict["coverage_radio"]

    # Select only neutron stars that fall into the sky region covered by the surveys.
    dictionary_coverage_database = {
        key: value[coverage_tot] for key, value in pop_dict.items()
    }

    # Add coverage for each survey to the dictionary.
    dictionary_coverage_database["coverage_radio"] = coverage_dict[
        "coverage_radio"
    ][coverage_tot]

    for survey_name in survey_radio_names:
        # Add the coverage data for each survey.
        coverage_key = f"coverage_radio_{survey_name}"
        dictionary_coverage_database[coverage_key] = coverage_dict[
            coverage_key
        ][coverage_tot]

    # Remove stars that do not fall into the total sky coverage.
    idx = pop_dict["idx"]
    out_coverage = np.invert(coverage_tot)
    idx_remove += idx[out_coverage].tolist()

    return dictionary_coverage_database, idx_remove

compute_surveys_coverage(surveys_radio, pop_dict, dist_cutoff)

Apply survey coverage criteria to filter a dynamic population dataset based on sky coverage of all surveys and a distance cutoff, and update the indices of entries to be removed.

Parameters:

Name	Type	Description	Default
surveys_radio	dict	A dictionary of radio survey objects, containing the information on the sky coverage.	required
pop_dict	dict	A dictionary containing the sky coordinates of a population of neutron stars.	required
dist_cutoff	float	The maximum heliocentric distance to include in the survey coverage.	required

Returns:

Type	Description
dict	A dictionary with the sky coverage information for all surveys.

Source code in mlpoppyns/simulator/multiband_surveys/surveys_wrapper.py

def compute_surveys_coverage(
    surveys_radio: dict,
    pop_dict: dict,
    dist_cutoff: float,
) -> dict:
    """
    Apply survey coverage criteria to filter a dynamic population dataset based on sky coverage of all surveys and a
    distance cutoff, and update the indices of entries to be removed.

    Args:
        surveys_radio (dict): A dictionary of radio survey objects, containing the information on the sky coverage.
        pop_dict (dict): A dictionary containing the sky coordinates of a population of neutron stars.
        dist_cutoff (float): The maximum heliocentric distance to include in the survey coverage.

    Returns:
        (dict): A dictionary with the sky coverage information for all surveys.
    """

    dist = pop_dict["dist"]
    dist_mask = dist < dist_cutoff

    survey_radio_names = list(surveys_radio.keys())
    coverage_survey_radio = {}

    for name in survey_radio_names:
        # Evaluate the sky coverage for each radio survey.
        coverage_survey_radio[name] = surveys_radio[name].sky_coverage(
            pop_dict["ra"],
            pop_dict["dec"],
            pop_dict["l"],
            pop_dict["b"],
        )

    # Combine the coverage masks of all selected radio surveys into a single mask.
    # It performs a logical OR (|) across all coverage arrays in coverage_survey_radio,
    # for each survey name in survey_radio_names.
    # The result is a single array where a position is True if it is covered by any survey.
    coverage_radio_tot = (
        functools.reduce(
            lambda a, b: a | b,
            (coverage_survey_radio[name] for name in survey_radio_names),
        )
    ) & dist_mask

    # Create a dictionary to save the coverage information.
    coverage_dict = {}
    coverage_dict["coverage_radio"] = coverage_radio_tot

    for survey_name in survey_radio_names:
        # Add the coverage data for each survey.
        coverage_key = f"coverage_radio_{survey_name}"
        coverage_dict[coverage_key] = coverage_survey_radio[survey_name]

    return coverage_dict

create_output_dataframe_surveys(dictionary_detected_radio)

Creates Pandas DataFrames containing the information on the detected neutron stars for each survey.

Parameters:

Name	Type	Description	Default
dictionary_detected_radio	dict	Dictionary containing detected neutron star properties for each radio survey.	required

Returns: (dict): A dictionary of DataFrames, one for each survey containing detected neutron stars' information.

Source code in mlpoppyns/simulator/multiband_surveys/surveys_wrapper.py

def create_output_dataframe_surveys(
    dictionary_detected_radio: dict,
) -> dict:
    """
    Creates Pandas DataFrames containing the information on the detected neutron stars for each survey.

    Args:
        dictionary_detected_radio (dict): Dictionary containing detected neutron star properties for each radio survey.
    Returns:
        (dict): A dictionary of DataFrames, one for each survey containing detected neutron stars' information.
    """

    # Initialize an empty dictionary to store the resulting DataFrames.
    dfs = {}

    # Defining the parameters and units that are common for all radio surveys.
    parameters_radio = [
        "idx",
        "age",
        "ra",
        "dec",
        "l",
        "b",
        "DM",
        "dist",
        "pm_ra",
        "pm_dec",
        "v_ls",
        "B",
        "chi",
        "P",
        "P_dot",
        "L_radio_bol",
        "S_radio_obs_mean",
        "S_radio_obs_mean_1400",
        "w_int",
        "w_eff",
        "tau_sc",
        "spectral_index",
    ]
    units_radio = [
        "",
        "[yr]",
        "[deg]",
        "[deg]",
        "[deg]",
        "[deg]",
        "[pc cm^-3]",
        "[kpc]",
        "[mas yr^-1]",
        "[mas yr^-1]",
        "[km s^-1]",
        "[G]",
        "[rad]",
        "[s]",
        "[s s^-1]",
        "[erg s^-1]",
        "[Jy]",
        "[Jy]",
        "[s]",
        "[s]",
        "[s]",
        "",
    ]

    # Loop over each survey's detected dictionary and generate the corresponding DataFrame.
    for survey_name, survey_data in dictionary_detected_radio.items():
        # Check if the survey is a combined survey like 'HTRU_low_mid'.
        if survey_name == "HTRU_low_mid":
            parameters_survey = parameters_radio + ["HTRU_low", "HTRU_mid"]
            units_survey = units_radio + ["", ""]
        else:
            parameters_survey = parameters_radio
            units_survey = units_radio

        # Build the DataFrame using the appropriate parameters and units.
        df = dfb.build_dataframe(survey_data, parameters_survey, units_survey)
        dfs[
            survey_name
        ] = df  # Store the DataFrame in the dictionary with survey_name as key.

    return dfs

initialize_all_surveys()

Initialize all surveys.

Returns:

Type	Description
SurveyData	A SurveyData object containing all survey data.

Source code in mlpoppyns/simulator/multiband_surveys/surveys_wrapper.py

def initialize_all_surveys() -> SurveyData:
    """
    Initialize all surveys.

    Returns:
        (SurveyData): A SurveyData object containing all survey data.
    """

    # Copy survey configuration data from the configuration file.
    surveys_cfg = cfg["surveys_radio"].copy()

    # Initialize radio surveys.
    surveys_radio, dictionary_detected_radio = initialize_radio_surveys()

    survey_data_class = SurveyData(
        surveys_cfg=surveys_cfg,
        surveys_radio=surveys_radio,
        n_detected_sim={survey: 0 for survey in surveys_cfg},
        n_detected_complete_sim={survey: 0 for survey in surveys_cfg},
        percentage_detected={survey: 0 for survey in surveys_cfg},
        n_created_at_match={survey: 0 for survey in surveys_cfg},
        n_detected_sim_at_match={survey: 0 for survey in surveys_cfg},
        batchsize_adjust_flags={0.9: False, 0.95: False},
        stop_flags={survey: False for survey in surveys_cfg},
        dictionary_detected_radio=dictionary_detected_radio,
    )

    return survey_data_class

initialize_radio_surveys()

Initialize and return radio survey objects and empty detection dictionaries.

Returns:

Type	Description
Tuple[dict, dict]	A tuple of dictionaries containing the following information: A dictionary of initialized radio survey objects, keyed by survey names (e.g., "PMPS", "SMPS", etc.). An empty dictionary for storing detected neutron star data for each of the survey.

Source code in mlpoppyns/simulator/multiband_surveys/surveys_wrapper.py

def initialize_radio_surveys() -> Tuple[dict, dict]:
    """
    Initialize and return radio survey objects and empty detection dictionaries.

    Returns:
        (Tuple[dict, dict]): A tuple of dictionaries containing the following information:

            - A dictionary of initialized radio survey objects, keyed by survey names (e.g., "PMPS", "SMPS", etc.).
            - An empty dictionary for storing detected neutron star data for each of the survey.
    """
    # Get the path to the software directory.
    base_path = pathlib.Path(cfg["path_to_software"])

    # Initialize survey objects.
    radio_surveys = {}

    for survey_name, survey_dict in cfg["surveys_radio"].items():
        if survey_name == "HTRU_low_mid":
            radio_surveys["HTRU_low"] = sr.SurveyRadio(
                str(base_path.joinpath(survey_dict["path_low"]))
            )
            radio_surveys["HTRU_mid"] = sr.SurveyRadio(
                str(base_path.joinpath(survey_dict["path_mid"]))
            )
        else:
            radio_surveys[survey_name] = sr.SurveyRadio(
                str(base_path.joinpath(survey_dict["path"]))
            )

    # Detection dictionary template.
    detection_template = {
        "age": [],
        "ra": [],
        "dec": [],
        "l": [],
        "b": [],
        "DM": [],
        "dist": [],
        "pm_ra": [],
        "pm_dec": [],
        "v_ls": [],
        "B": [],
        "chi": [],
        "P": [],
        "P_dot": [],
        "L_radio_bol": [],
        "S_radio_obs_mean": [],
        "S_radio_obs_mean_1400": [],
        "w_int": [],
        "w_eff": [],
        "tau_sc": [],
        "spectral_index": [],
        "idx": [],
    }

    # Initialize detection dictionaries for each survey.
    detection_dictionaries = {
        survey_name: detection_template.copy()
        for survey_name in cfg["surveys_radio"].keys()
    }

    # Add flags for HTRU low and mid surveys.
    if "HTRU_low_mid" in cfg["surveys_radio"]:
        detection_dictionaries["HTRU_low_mid"].update(
            {"HTRU_low": [], "HTRU_mid": []}
        )

    return radio_surveys, detection_dictionaries

radio_detection(radio_surveys, dictionary_radio, intercepted_radio, logger, full_population=False)

Simulate radio detections for various surveys and update the dictionaries with the properties of detected neutron stars.

Parameters:

Name	Type	Description	Default
radio_surveys	dict	Dictionary containing the radio survey objects.	required
dictionary_radio	dict	Dictionary with properties of radio pulsars.	required
intercepted_radio	ndarray	Boolean mask to select radio pulsars whose beam intercept our line of sight.	required
logger	Logger	Logger object for logging.	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 the properties of detected pulsars for each survey.

Source code in mlpoppyns/simulator/multiband_surveys/surveys_wrapper.py

def radio_detection(
    radio_surveys: dict,
    dictionary_radio: dict,
    intercepted_radio: np.ndarray,
    logger: logging.Logger,
    full_population: bool = False,
) -> dict:
    """
    Simulate radio detections for various surveys and update the dictionaries with the properties
    of detected neutron stars.

    Args:
        radio_surveys (dict): Dictionary containing the radio survey objects.
        dictionary_radio (dict): Dictionary with properties of radio pulsars.
        intercepted_radio (np.ndarray): Boolean mask to select radio pulsars whose beam intercept our line of sight.
        logger (logging.Logger): Logger object for logging.
        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 the properties of detected pulsars for each survey.
    """
    # Initialize variables for HTRU_low and HTRU_mid.
    detected_HTRU_low = np.array([])
    w_eff_low = None
    S_radio_obs_mean_low = None
    S_radio_obs_mean_1400_low = None

    detected_HTRU_mid = np.array([])
    w_eff_mid = None
    S_radio_obs_mean_mid = None
    S_radio_obs_mean_1400_mid = None

    # Process each survey.
    detected_dictionaries = {}

    for survey_name in radio_surveys:
        (
            detected_mask,
            w_eff,
            S_radio_obs_mean,
            S_radio_obs_mean_1400,
        ) = radio_surveys[survey_name].detected_radio_population(
            dictionary_radio["w_int"],
            dictionary_radio["DM"],
            dictionary_radio["P"],
            intercepted_radio,
            dictionary_radio[f"coverage_radio_{survey_name}"],
            dictionary_radio["l"],
            dictionary_radio["b"],
            dictionary_radio["dec"],
            dictionary_radio["S_radio_bol"],
            dictionary_radio["spectral_index"],
            dictionary_radio["tau_sc"],
        )
        detected_dictionaries[survey_name] = update_filtered_dictionary(
            dictionary_radio,
            detected_mask,
            w_eff=w_eff,
            S_radio_obs_mean=S_radio_obs_mean,
            S_radio_obs_mean_1400=S_radio_obs_mean_1400,
        )
        if full_population:
            detected_dictionaries[survey_name].pop("intercepted_radio", None)

            fraction_detected = len(detected_mask[detected_mask]) / len(
                detected_mask
            )
            logger.info(
                f"Fraction of detected pulsars by {survey_name}: {fraction_detected}"
            )

        # Save the properties for the HTRU low and mid surveys separately.
        if survey_name == "HTRU_low":
            detected_HTRU_low = detected_mask
            w_eff_low = w_eff
            S_radio_obs_mean_low = S_radio_obs_mean
            S_radio_obs_mean_1400_low = S_radio_obs_mean_1400

        elif survey_name == "HTRU_mid":
            detected_HTRU_mid = detected_mask
            w_eff_mid = w_eff
            S_radio_obs_mean_mid = S_radio_obs_mean
            S_radio_obs_mean_1400_mid = S_radio_obs_mean_1400

    if (
        "HTRU_low" in radio_surveys.keys()
        and "HTRU_mid" in radio_surveys.keys()
    ):
        # Since the sky coverage of the HTRU mid and low surveys overlap, we remove those stars from the mid
        # survey that are already in the low survey in order to not double count individual objects.
        detected_HTRU_low_mid = detected_HTRU_low | detected_HTRU_mid
        detected_dictionaries["HTRU_low_mid"] = update_filtered_dictionary(
            dictionary_radio,
            detected_HTRU_low_mid,
            w_eff=np.where(detected_HTRU_low, w_eff_low, w_eff_mid),
            S_radio_obs_mean=np.where(
                detected_HTRU_low, S_radio_obs_mean_low, S_radio_obs_mean_mid
            ),
            S_radio_obs_mean_1400=np.where(
                detected_HTRU_low,
                S_radio_obs_mean_1400_low,
                S_radio_obs_mean_1400_mid,
            ),
            HTRU_low=detected_HTRU_low,
            HTRU_mid=detected_HTRU_mid,
            idx=dictionary_radio["idx"],
        )
        if full_population:
            detected_dictionaries["HTRU_low_mid"].pop(
                "intercepted_radio", None
            )

        # Remove the dictionaries containing the results for the individual HTRU low and mid surveys,
        # as we only require the combined detections determined above.
        del detected_dictionaries["HTRU_low"]
        del detected_dictionaries["HTRU_mid"]

    return detected_dictionaries

update_filtered_dictionary(dict_to_update, mask, **kwargs)

This function updates a dictionary to include for each key only the values corresponding to a given Boolean mask.

Parameters:

Name	Type	Description	Default
dict_to_update	dict	The original dictionary containing properties of neutron stars.	required
mask	ndarray	Boolean mask indicating which elements for each key in dict_to_update have to be included.	required
**kwargs	ndarray	Additional property values provided as keyword arguments. These properties will also be filtered using the mask.	{}

Returns:

Type	Description
dict	A new dictionary containing only the filtered values from dict_to_update and combining the keys already present in the original dictionary with the ones provided in kwargs.

Source code in mlpoppyns/simulator/multiband_surveys/surveys_wrapper.py

def update_filtered_dictionary(
    dict_to_update: dict, mask: np.ndarray, **kwargs: np.ndarray
) -> dict:
    """
    This function updates a dictionary to include for each key only the values corresponding to a given Boolean mask.

    Args:
        dict_to_update (dict): The original dictionary containing properties of neutron stars.
        mask (np.ndarray): Boolean mask indicating which elements for each key in `dict_to_update` have to be included.
        **kwargs (np.ndarray): Additional property values provided as keyword arguments.
            These properties will also be filtered using the `mask`.

    Returns:
        (dict): A new dictionary containing only the filtered values from `dict_to_update` and combining
            the keys already present in the original dictionary with the ones provided in `kwargs`.
    """

    # Extract keys from the dictionary that has to be updated.
    properties = list(dict_to_update.keys())

    # Filtering the values in the original dictionary.
    filtered_dict = {
        prop: dict_to_update[prop][mask].tolist() for prop in properties
    }

    # Add and filter the properties specified in `kwargs` to another dictionary.
    additional_filtered_dict = {
        key: value[mask].tolist() for key, value in kwargs.items()
    }

    # Combine the two dictionaries.
    combined_dict = filtered_dict | additional_filtered_dict

    return combined_dict

update_survey_data(survey_data_class, pop_detected_dict_update, survey_type, n_created, idx_remove, logger)

This function updates the survey data dictionaries to include only the properties of neutron stars that are detected in the surveys.

Parameters:

Name	Type	Description	Default
survey_data_class	SurveyData	The SurveyData dataclass containing the data of all neutron star surveys.	required
pop_detected_dict_update	dict	A dictionary containing the properties of neutron stars that have been detected in the surveys.	required
survey_type	str	A string specifying which survey to update, "Radio" or "X-ray".	required
n_created	int	The total number of neutron stars created so far.	required
idx_remove	list	A list of indices of the neutron stars that have been detected and have to be removed from the dynamical database.	required
logger	Logger	A logger instance to log messages.	required

Source code in mlpoppyns/simulator/multiband_surveys/surveys_wrapper.py

def update_survey_data(
    survey_data_class: SurveyData,
    pop_detected_dict_update: dict,
    survey_type: str,
    n_created: int,
    idx_remove: list,
    logger: logging.Logger,
) -> None:
    """
    This function updates the survey data dictionaries to include only the properties of neutron stars that are
    detected in the surveys.

    Args:
        survey_data_class (SurveyData): The SurveyData dataclass containing the data of all neutron star surveys.
        pop_detected_dict_update (dict): A dictionary containing the properties of neutron stars that have been
            detected in the surveys.
        survey_type (str): A string specifying which survey to update, "Radio" or "X-ray".
        n_created (int): The total number of neutron stars created so far.
        idx_remove (list): A list of indices of the neutron stars that have been detected and have to be removed
            from the dynamical database.
        logger (logging.Logger): A logger instance to log messages.
    """

    surveys_cfg = survey_data_class.surveys_cfg
    n_detected_sim = survey_data_class.n_detected_sim
    n_detected_complete_sim = survey_data_class.n_detected_complete_sim
    stop_flags = survey_data_class.stop_flags
    n_detected_sim_at_match = survey_data_class.n_detected_sim_at_match
    n_created_at_match = survey_data_class.n_created_at_match
    dictionary_detected_radio = survey_data_class.dictionary_detected_radio

    for survey in pop_detected_dict_update:
        if survey_type == "radio":
            # Check the number of detected pulsars for each survey.
            n_detected_sim[survey] += len(
                pop_detected_dict_update[survey]["age"]
            )
            n_detected_complete_sim[survey] += len(
                pop_detected_dict_update[survey]["age"]
            )
            logger.info(
                f"Total number of neutron stars detected by {survey}: {n_detected_sim[survey]}"
            )
            # Since we assume that the radio surveys are complete, i.e.,
            # n_detected_complete_sim = n_detected_sim, if the number of simulated detected pulsars
            # matches the real one, store the value of created neutron stars.
            if (
                n_detected_complete_sim[survey]
                >= surveys_cfg[survey]["detected_real"]
                and not stop_flags[survey]
            ):
                n_detected_sim_at_match[survey] = n_detected_sim[survey]
                stop_flags[survey] = True
                n_created_at_match[survey] = n_created

            # Update the detection dictionaries.
            dictionary_detected_radio[survey] = {
                key: value + pop_detected_dict_update[survey][key]
                for key, value in dictionary_detected_radio[survey].items()
            }
            # Update the index list to remove the stars that have been detected from the dynamical database.
            idx_det = list(dictionary_detected_radio[survey]["idx"])
            idx_remove += idx_det

        else:
            logger.error(f"Unknown survey type: {survey_type}")
            sys.exit(1)