class HealpixFitsMap:
    """ Class to handle fits file containing healpix maps """

    def __init__(
        self,
        input_filename,
        do_pol=True,
        verbose=False,
        fwhm_in=0.0,
        fwhm_in2=None,
        nside_in=16,
        lmax=None,
        map_seed=53543,
        no_ileak=False,
        no_quleak=False,
        compute_derivatives=None,
        derivatives_type="T",
        ext_map_gal=False,
    ):
        """

        Parameters
        ----------
        input_filename : string, or list of strings
            Either fits file containing the sky maps (data will just be
            loaded), or CAMB lensed cl file (.dat) containing lensed
            power spectra with order ell, TT, EE, BB, TE (maps will be
            created on-the-fly), or a list of 3 fits files containing alms
            (maps will be created on-the-fly).
        do_pol : bool, optional
            If True, load temperature and polarisation. Temperature only
            otherwise. Default is True.
        verbose : bool, optional
            If True, print out plenty of useless messages.
        fwhm_in : float, optional
            If input_filename is a CAMB lensed cl file, the generated maps will
            be convolved with a beam having this fwhm_in. In arcmin.
            No effect if you provide maps directly.
        fwhm_in2 : float, optional
            If provided, will generate another set of I, Q, U with this
            resolution (useful for dichroic detectors). Default is None.
            No effect if you provide maps directly.
        nside_in : int, optional
            If input_filename is a CAMB lensed cl file, the maps will be
            generated at a resolution nside_in. No effect if you provide maps
            directly.
        lmax : None or int, optional
            Maximum multipole when creating a map from cl. If none, it
            is set automatically to 2*nside_in
        map_seed : int, optional
            If input_filename is a CAMB lensed cl file, this is the seed used
            to create the maps. No effect if you provide maps directly.
        no_ileak : bool, optional
            If True, load temperature and polarisation, but set the temperature
            to zero to avoid leakages.
        no_quleak : bool, optional
            If True, load temperature and polarisation, but set the
            polarisation to zero to avoid leakages.
        ext_map_gal : bool, optional
            Set it to True if you are reading a map in Galactic coordinate.
            (Planck maps for example).
        compute_derivatives : bool, optional
            If True, return derivatives of relevant Stokes component according
            to the derivatives_type option.
            Be sure that you have enough memory!
        derivatives_type : str, optional
            If 'T' ('P') present in the string return 1st and 2nd derivatives
            of the input temperature (Q and U) map (t=theta, p=phi).
            Note that d/dp is already divided by sin(theta).
            Options:
            - 'T1' will compute 1st derivatives dI/dt, dI/dp alone
            - 'P1' will compute 1st derivatives dQ/dt, dU/dt, dQ/dp, dU/dp.
            Be sure that you have enough memory!
        """
        self.input_filename = input_filename
        self.do_pol = do_pol
        self.verbose = verbose
        self.no_ileak = no_ileak
        self.no_quleak = no_quleak
        self.ext_map_gal = ext_map_gal
        self.fwhm_in = fwhm_in
        self.fwhm_in2 = fwhm_in2
        self.nside_in = nside_in
        if lmax is None:
            self.lmax = 2 * self.nside_in
        else:
            self.lmax = lmax
        self.map_seed = map_seed
        self.compute_derivatives = compute_derivatives
        self.derivatives_type = derivatives_type
        self.I = None
        self.Q = None
        self.U = None

        self.I2 = None
        self.Q2 = None
        self.U2 = None

        fromalms = False
        if type(self.input_filename) == list:
            if self.verbose:
                print("Reading sky maps from alms file...")
            self.load_healpix_fits_map_from_alms()
            fromalms = True
        elif self.input_filename[-4:] == ".dat":
            if self.verbose:
                print("Creating sky maps from cl file...")
            self.create_healpix_fits_map()
        elif self.input_filename[-5:] == ".fits":
            if self.verbose:
                print("Reading sky maps from fits file...")
            self.load_healpix_fits_map()
        else:
            raise IOError(
                """Input file not understood! Should be either a
                fits file containing the sky maps
                (data will be loaded), or a CAMB lensed cl file
                (.dat) containing lensed power spectra with
                order ell, TT, EE, BB, TE
                (maps will be created on-the-fly)."""
            )

        self.set_leakage_to_zero()

        if self.compute_derivatives:
            if "T" in self.derivatives_type:
                self.compute_intensity_derivatives(fromalm=fromalms)
            if "P" in self.derivatives_type:
                self.compute_pol_derivatives(fromalm=fromalms)

    def load_healpix_fits_map(self, force=False):
        """
        Load from disk into memory a sky map.

        Not updated for dichroic for the moment.

        Parameters
        ----------
        force : bool
            If true, force to load the maps in memory even if it is already
            loaded. Default is False.

        Examples
        ----------
        Let's generate fake data
        >>> filename = 'myfits_to_test_.fits'
        >>> write_dummy_map(filename, nside=16)

        Let's now read the data
        >>> hpmap = HealpixFitsMap(input_filename=filename)
        >>> print(hpmap.nside)
        16

        If the data is already loaded, it won't reload it by default
        >>> hpmap.load_healpix_fits_map()
        External data already present in memory

        But you can force it
        >>> hpmap.load_healpix_fits_map(force=True)
        """
        if self.I is None or force:
            if self.do_pol:
                self.I, self.Q, self.U = hp.read_map(
                    self.input_filename, (0, 1, 2), verbose=self.verbose
                )
            else:
                self.I = hp.read_map(
                    self.input_filename, field=0, verbose=self.verbose
                )
            self.nside = hp.npix2nside(len(self.I))
        else:
            print("External data already present in memory")

    def load_healpix_fits_map_from_alms(self, force=False):
        """
        Load from disk into memory alms and make sky maps.

        Parameters
        ----------
        force : bool
            If true, force to load the maps in memory even if it is already
            loaded. Default is False.

        Examples
        ----------
        Let's generate fake data
        >>> np.random.seed(548397)
        >>> sky_maps = create_sky_map('s4cmb/data/test_data_set_lensedCls.dat')
        >>> alms = hp.map2alm(sky_maps)
        >>> filenames = ['myalms_to_test_tlm.fits', 'myalms_to_test_elm.fits',
        ...     'myalms_to_test_blm.fits']
        >>> for fn, alm in zip(filenames, alms):
        ...     hp.write_alm(fn, alm)

        Let's now read the data
        >>> hpmap = HealpixFitsMap(input_filename=filenames)
        >>> print(hpmap.nside)
        16

        If the data is already loaded, it won't reload it by default
        >>> hpmap.load_healpix_fits_map_from_alms()
        External data already present in memory

        But you can force it
        >>> hpmap.load_healpix_fits_map_from_alms(force=True)

        You can also generate 2 sets of maps with different resolution which
        is useful for dichroic detectors
        >>> hpmap = HealpixFitsMap(input_filename=filenames, fwhm_in=3.5,
        ...     fwhm_in2=1.8, nside_in=16,)
        >>> hasattr(hpmap, 'Q2')
        True
        """
        if self.I is None or force:
            if self.do_pol:
                tlm = hp.read_alm(self.input_filename[0])
                elm = hp.read_alm(self.input_filename[1])
                blm = hp.read_alm(self.input_filename[2])

                self.I, self.Q, self.U = hp.alm2map(
                    [tlm, elm, blm],
                    nside=self.nside_in,
                    pixwin=False,
                    fwhm=self.fwhm_in / 60.0 * np.pi / 180.0,
                    sigma=None,
                    pol=True,
                    inplace=False,
                    verbose=self.verbose,
                )
                if self.fwhm_in2 is not None:
                    self.I2, self.Q2, self.U2 = hp.alm2map(
                        [tlm, elm, blm],
                        nside=self.nside_in,
                        pixwin=False,
                        fwhm=self.fwhm_in2 / 60.0 * np.pi / 180.0,
                        sigma=None,
                        pol=True,
                        inplace=False,
                        verbose=self.verbose,
                    )
            else:
                tlm = hp.read_alm(self.input_filename[0])

                self.I = hp.alm2map(
                    tlm,
                    nside=self.nside_in,
                    pixwin=False,
                    fwhm=self.fwhm_in / 60.0 * np.pi / 180.0,
                    sigma=None,
                    pol=False,
                    inplace=False,
                    verbose=self.verbose,
                )
                if self.fwhm_in2 is not None:
                    self.I2 = hp.alm2map(
                        tlm,
                        nside=self.nside_in,
                        pixwin=False,
                        fwhm=self.fwhm_in2 / 60.0 * np.pi / 180.0,
                        sigma=None,
                        pol=False,
                        inplace=False,
                        verbose=self.verbose,
                    )

            self.nside = hp.npix2nside(len(self.I))
        else:
            print("External data already present in memory")

    def create_healpix_fits_map(self, force=False):
        """
        Create sky maps from cl file.
        Do nothing if data already presents in the memory.

        Parameters
        ----------
        force : bool
            If true, force to recreate the maps in memory even
            if it is already loaded. Default is False.

        Examples
        ----------
        Let's generate the map from a CAMB file
        >>> filename = 's4cmb/data/test_data_set_lensedCls.dat'
        >>> hpmap = HealpixFitsMap(input_filename=filename, fwhm_in=3.5,
        ...     nside_in=16, map_seed=489237)
        >>> print(hpmap.nside)
        16

        If the data is already loaded, it won't reload it by default
        >>> hpmap.create_healpix_fits_map()
        External data already present in memory

        But you can force it
        >>> hpmap.create_healpix_fits_map(force=True)

        You can also load 2 sets of maps with different resolution, which
        is useful for dichroic detectors
        >>> filename = 's4cmb/data/test_data_set_lensedCls.dat'
        >>> hpmap = HealpixFitsMap(input_filename=filename, fwhm_in=3.5,
        ...     fwhm_in2=1.8, nside_in=16, map_seed=489237)
        >>> hasattr(hpmap, 'I2')
        True
        """
        if self.I is None or force:
            if self.do_pol:
                self.I, self.Q, self.U = create_sky_map(
                    self.input_filename,
                    nside=self.nside_in,
                    FWHM=self.fwhm_in,
                    seed=self.map_seed,
                    lmax=self.lmax,
                )
                if self.fwhm_in2 is not None:
                    self.I2, self.Q2, self.U2 = create_sky_map(
                        self.input_filename,
                        nside=self.nside_in,
                        FWHM=self.fwhm_in2,
                        seed=self.map_seed,
                        lmax=self.lmax,
                    )
            else:
                self.I = create_sky_map(
                    self.input_filename,
                    nside=self.nside_in,
                    FWHM=self.fwhm_in,
                    seed=self.map_seed,
                    lmax=self.lmax,
                )
                if self.fwhm_in2 is not None:
                    self.I2 = create_sky_map(
                        self.input_filename,
                        nside=self.nside_in,
                        FWHM=self.fwhm_in2,
                        seed=self.map_seed,
                        lmax=self.lmax,
                    )
            self.nside = hp.npix2nside(len(self.I))
        else:
            print("External data already present in memory")

    def set_leakage_to_zero(self):
        """
        Remove either I, Q or U to remove possible leakages

        Examples
        ----------
        Test with no input intensity
        >>> write_dummy_map('myfits_to_test_.fits')
        >>> hpmap = HealpixFitsMap('myfits_to_test_.fits', no_ileak=True)
        >>> print(hpmap.I)
        [ 0.  0.  0. ...,  0.  0.  0.]

        Test with no input polarisation
        >>> write_dummy_map('myfits_to_test_.fits')
        >>> hpmap = HealpixFitsMap('myfits_to_test_.fits', no_quleak=True)
        >>> print(hpmap.Q, hpmap.U)
        [ 0.  0.  0. ...,  0.  0.  0.] [ 0.  0.  0. ...,  0.  0.  0.]

        If you have two sets of maps, it will remove leakages from the two sets
        >>> filename = 's4cmb/data/test_data_set_lensedCls.dat'
        >>> hpmap = HealpixFitsMap(input_filename=filename, fwhm_in=3.5,
        ...     fwhm_in2=1.8, nside_in=16, map_seed=489237,
        ...     no_ileak=True, no_quleak=True)
        >>> print(hpmap.I, hpmap.I2)
        [ 0.  0.  0. ...,  0.  0.  0.] [ 0.  0.  0. ...,  0.  0.  0.]
        """
        # Set temperature to zero to avoid I->QU leakage
        if self.no_ileak:
            self.I[:] = 0.0
            if self.I2 is not None:
                self.I2[:] = 0.0

        # Set polarisation to zero to avoid QU leakage
        if self.no_quleak:
            if self.Q is not None:
                self.Q[:] = 0.0
            if self.U is not None:
                self.U[:] = 0.0

            if self.Q2 is not None:
                self.Q2[:] = 0.0
            if self.U2 is not None:
                self.U2[:] = 0.0

    def compute_intensity_derivatives(self, fromalm=False):
        """
        Compute derivatives of the input temperature map (healpix).
        Not updated for dichroic for the moment.

        Parameters
        ----------
        fromalm : bool, optional
            If True, loads alm file from disk instead of fourier
            transform the input map. Automatically turns True if you input
            alm files. False otherwise.

        Examples
        ----------
        >>> filename = 's4cmb/data/test_data_set_lensedCls.dat'
        >>> hpmap = HealpixFitsMap(input_filename=filename, fwhm_in=3.5,
        ...     nside_in=16, compute_derivatives=True, map_seed=489237)
        >>> hasattr(hpmap, 'dIdp')
        True

        """
        if fromalm:
            alm = hp.read_alm(self.input_filename[0])
        else:
            alm = hp.map2alm(self.I, self.lmax)

        # lmax = hp.Alm.getlmax(alm.size)
        if "T1" in self.derivatives_type:
            junk, self.dIdt, self.dIdp = hp.alm2map_der1(
                alm, self.nside_in, self.lmax
            )
        else:
            # computes first and second derivative as derivatives of spin-1
            # transform of a scalar field with _1Elm=sqrt(l(l+1))Ilm _1Blm=0
            l = np.arange(self.lmax + 1)
            grad = np.sqrt(l * (l + 1))
            curl = np.zeros_like(alm)
            dervs = alm2map_spin_der1([hp.almxfl(alm, grad), curl], self.nside_in, 1)
            self.dIdt = dervs[0][0]
            self.dIdp = dervs[0][1]
            self.d2Id2t = dervs[1][0]
            self.d2Id2p = dervs[2][1]
            self.d2Idpdt = dervs[2][0]

    def compute_pol_derivatives(self, fromalm=False):
        """
        Compute derivatives of the input polarization components (healpix).
        Not updated for dichroic for the moment.

        Parameters
        ----------
        fromalm : bool, optional
            If True, loads alm file from disk instead of fourier
            transform the input map. Automatically turns True if you input
            alm files. False otherwise.

        Examples
        ----------
        >>> filename = 's4cmb/data/test_data_set_lensedCls.dat'
        >>> hpmap = HealpixFitsMap(input_filename=filename, fwhm_in=3.5,
        ...     nside_in=16, compute_derivatives=True, map_seed=489237)
        >>> hasattr(hpmap, 'dIdp')
        True

        """
        if fromalm:
            Elm = hp.read_alm(self.input_filename[1])
            Blm = hp.read_alm(self.input_filename[2])
        else:
            alm = hp.map2alm([self.I, self.Q, self.U], self.lmax)
            Elm = alm[1]
            Blm = alm[2]
        # lmax = hp.Alm.getlmax(Elm.size)
        if "P1" in self.derivatives_type:
            out = alm2map_spin_der1([Elm, Blm], self.nside_in, 2)
            self.dQdt = out[1][0]
            self.dUdt = out[1][1]
            self.dQdp = out[2][0]
            self.dUdp = out[2][1]
        else:
            warnings.warn("""
            Computation of second order polarization derivatives not
            implemented yet. Set to 0.
            """)

            out = alm2map_spin_der1([Elm, Blm], self.nside_in, 2)
            self.dQdt = out[1][0]
            self.dUdt = out[1][1]
            self.dQdp = out[2][0]
            self.dUdp = out[2][1]
            self.d2Qd2t = np.zeros_like(self.dQdt)
            self.d2Qd2p = np.zeros_like(self.dQdt)
            self.d2Qdpdt = np.zeros_like(self.dQdt)
            self.d2Ud2t = np.zeros_like(self.dQdt)
            self.d2Ud2p = np.zeros_like(self.dQdt)
            self.d2Udpdt = np.zeros_like(self.dQdt) class Hardware:
    """ Class to load all the hardware and models of the instrument in once """

    def __init__(
        self,
        ncrate=1,
        ndfmux_per_crate=1,
        nsquid_per_mux=1,
        npair_per_squid=4,
        fp_size=60.0,
        fwhm=3.5,
        beam_seed=58347,
        projected_fp_size=3.0,
        pm_name="5params",
        type_hwp="CRHWP",
        freq_hwp=2.0,
        angle_hwp=0.0,
        verbose=False,
    ):
        """
        This class creates the data used to model the instrument:
        * focal plane
        * pointing model parameters of the telescope
        * beam parameters of the bolometers
        * polarisation angle of the bolometers

        Parameters
        ----------
        ncrate : int
            Number of crate plate.
        ndfmux_per_crate : int
            Number of MUX board per crate.
        nsquid_per_mux : int
            Number of SQUID per MUX.
        npair_per_squid : int
            Number of pair of bolometers per SQUID.
        fp_size : float, optional
            The size of the focal plane in cm. Default is 60 cm.
        fwhm : float, optional
            Full Width Half Maximum of the beam (in arcmin).
            Default = 3.5 arcmin.
        beam_seed : int
            Seed used to generate angle of rotation of beam axes.
            Default is 58347.
        projected_fp_size : float, optional
            Size of the focal plane on the sky (in degree). This has to
            do with the size of the mirror and the beam size.
            Default = 3 degrees.
        pm_name : string, optional
            The pointing model to load. Currently, only the five-parameter
            pointing model (Mangum 2001) is implemented (pm_name = 5params).
        verbose : boolean, optional
            If True, print out a number of useful comments for verboseging.
        type_hwp : string, optional
            Choose between CRWHP or stepped.
        freq_hwp : float, optional
            If type_hwp=CRHWP, then freq_hwp sets the spin frequency of the
            HWP [Hz].
        angle_hwp : float, optional
            If type_hwp=CRHWP, angle_hwp corresponds to the starting position.
            If type_hwp=stepped, angle_hwp corresponds to the step size
            which is updated after each scan (CES). [deg]

        Examples
        ----------
        >>> instrument = Hardware()
        """
        self.focal_plane = FocalPlane(
            ncrate,
            ndfmux_per_crate,
            nsquid_per_mux,
            npair_per_squid,
            fp_size,
            verbose,
        )

        self.beam_model = BeamModel(
            self.focal_plane, fwhm, beam_seed, projected_fp_size, verbose
        )

        self.pointing_model = PointingModel(pm_name)

        self.half_wave_plate = HalfWavePlate(type_hwp, freq_hwp, angle_hwp)

    def make_dichroic(
        self, fwhm=1.8, beam_seed=58347, projected_fp_size=3.0, shift_angle=45
    ):
        """
        Add a layer of detectors on top of the existing ones
        to have dichroic detectors. The new detectors will sit on top of the
        existing ones (that is same number and same location on the focal
        plane). The only difference is the beam size for the moment (and the
        projected size of the fp on sky).

        Parameters
        ----------
        fwhm : float, optional
            Full Width Half Maximum of the beam (in arcmin).
            Default = 3.5 arcmin.
        beam_seed : int
            Seed used to generate angle of rotation of beam axes.
            Default is 58347.
        projected_fp_size : float, optional
            Size of the focal plane on the sky (in degree). This has to
            do with the size of the mirror and the beam size.
            Default = 3 degrees.
        shift_angle : int, optional
            Shift the polarisation angle with respect to the ones from the
            first frequency channel.

        Examples
        ----------
        >>> instrument = Hardware()
        >>> instrument.make_dichroic()
        >>> hasattr(instrument, 'focal_plane2')
        True
        """
        self.focal_plane2 = copy.copy(self.focal_plane)

        # Shift the polarisation angles (and clip it btw 0 and 360 deg)
        self.focal_plane2.bolo_polangle = (
            np.array(self.focal_plane2.bolo_polangle) + shift_angle
        ) % 360

        self.beam_model2 = BeamModel(
            self.focal_plane2, fwhm, beam_seed, projected_fp_size
        ) class ScanningStrategy:
    """ Class to handle the scanning strategy of the telescope """

    def __init__(
        self,
        nces=12,
        start_date="2013/1/1 00:00:00",
        telescope_longitude="-67:46.816",
        telescope_latitude="-22:56.396",
        telescope_elevation=5200.0,
        name_strategy="deep_patch",
        sampling_freq=30.0,
        sky_speed=0.4,
        ut1utc_fn="s4cmb/data/ut1utc.ephem",
        language="python",
        verbose=False,
    ):
        """
        A scanning strategy consists in defining the site of observation
        on earth for which we will make the observation, the region
        of the sky to observe, and the schedule of observations.

        Parameters
        ----------
        nces : int, optional
            Number of scans to generate.
        start_date : string, optional
            Starting date for observations. The format is: YYYY/M/D HH:MM:SS.
        telescope_longitude : string, optional
            Longitute (angle) of the telescope. String form: 0:00:00.0.
        telescope_latitude : string, optional
            Latitude (angle) of the telescope. String form: 0:00:00.0.
        telescope_elevation : float, optional
            Height above sea level (in meter).
        name_strategy : string, optional
            Name of a pre-defined scanning strategy to define the boundaries
            of the scan: elevation, azimuth, and time. Only available for the:
            moment name_strategy = deep_patch.
        sampling_freq : float, optional
            Sampling frequency of the bolometers in Hz.
        sky_speed : float, optional
            Azimuth speed of the telescope in deg/s.
        ut1utc_fn : string, optional
            File containing time correction to UTC (leap second).
            See http://tycho.usno.navy.mil/leapsec.html.
            This is not used here, but pass to the pointing module later on.
        language : string, optional
            Language used for core computations. For big experiments, the
            computational time can be big, and some part of the code can be
            speeded up by interfacing python with Fortran.
            Default is python (i.e. no interfacing and can be slow).
            Choose language=fortran otherwise. Note that for fortran codes you
            need first to compile it. See the setup.py or
            the provided Makefile.
        verbose : bool
            If True, print out several messages to ease the debug.
            Default is False.

        """
        self.nces = nces
        self.start_date = start_date
        self.name_strategy = name_strategy
        self.sampling_freq = sampling_freq
        self.sky_speed = sky_speed
        self.language = language

        self.ut1utc_fn = ut1utc_fn
        if not os.path.isfile(self.ut1utc_fn):
            # server in mantainance mode as of 202007. Message updated
            # url = 'http://tycho.usno.navy.mil/leapsec.html'
            url = 's4cmb.scanning_strategy.update_ut1utc'
            msg = 'The path {} does not point to a valid file! see ' + \
                's4cmb/data/ut1utc.ephem provided with the package. ' + \
                'For more information, see {}.'
            raise Exception(msg.format(self.ut1utc_fn, url))

        self.verbose = verbose

        self.telescope_location = self.define_telescope_location(
            telescope_longitude, telescope_latitude, telescope_elevation
        )

        self.define_boundary_of_scan()

    def define_telescope_location(
        self,
        telescope_longitude="-67:46.816",
        telescope_latitude="-22:56.396",
        telescope_elevation=5200.0,
    ):
        """
        Routine to define the site of observation on earth for which
        positions are to be computed. The location of the Polarbear telescope
        is entered as default.

        Parameters
        ----------
        telescope_longitude : str
            Longitute (angle) of the telescope. String form: 0:00:00.0.
        telescope_latitude : str
            Latitude (angle) of the telescope. String form: 0:00:00.0.
        telescope_elevation : float
            Height above sea level (in meter).

        Returns
        ----------
        location : Observer instance
            An `Observer` instance allows you to compute the positions of
            celestial bodies as seen from a particular latitude and longitude
            on the Earth's surface.

        Examples
        ----------
        >>> scan = ScanningStrategy()
        >>> telescope_location = scan.define_telescope_location()
        >>> telescope_location.elevation
        5200.0

        """
        location = ephem.Observer()
        location.long = telescope_longitude
        location.lat = telescope_latitude
        location.elevation = telescope_elevation

        return location

    def define_boundary_of_scan(self):
        """
        Given a pre-defined scanning strategy,
        define the boundaries of the scan: elevation, azimuth, and time.
        For a custom usage (advanced users), modify this routine.

        Examples
        ----------
        >>> scan = ScanningStrategy(name_strategy='deep_patch')
        >>> scan.elevation # doctest: +NORMALIZE_WHITESPACE
        [30.0, 45.5226, 47.7448, 49.967,
         52.1892, 54.4114, 56.6336, 58.8558,
         61.078, 63.3002, 65.5226, 35.2126]

        Only few scanning strategies are currently available:
        >>> scan = ScanningStrategy(name_strategy='shallow_patch')
        >>> scan = ScanningStrategy(name_strategy='small_aperture')
        >>> scan = ScanningStrategy(name_strategy='custom')
        >>> scan = ScanningStrategy(name_strategy='toto')
        ... # doctest: +NORMALIZE_WHITESPACE, +ELLIPSIS
        Traceback (most recent call last):
         ...
        ValueError: Only name_strategy = deep_patch or shallow_patch or
        custom are currently available. For another usage (advanced users),
        modify this routine.

        >>> scan.allowed_scanning_strategies
        ['deep_patch', 'shallow_patch', 'small_aperture', 'custom']
        """
        self.allowed_scanning_strategies = [
            "deep_patch",
            "shallow_patch",
            "small_aperture",
            "custom",
        ]

        if self.name_strategy == "deep_patch":
            self.elevation = [
                30.0,
                45.5226,
                47.7448,
                49.967,
                52.1892,
                54.4114,
                56.6336,
                58.8558,
                61.078,
                63.3002,
                65.5226,
                35.2126,
            ]

            self.az_min = [
                134.2263,
                162.3532,
                162.3532,
                162.3532,
                162.3532,
                162.3532,
                162.3532,
                162.3532,
                162.3532,
                162.3532,
                162.3532,
                204.7929,
            ]

            self.az_max = [
                154.2263,
                197.3532,
                197.3532,
                197.3532,
                197.3532,
                197.3532,
                197.3532,
                197.3532,
                197.3532,
                197.3532,
                197.3532,
                224.7929,
            ]

            self.begin_LST = [
                "17:07:54.84",
                "22:00:21.76",
                "22:00:21.76",
                "22:00:21.76",
                "22:00:21.76",
                "22:00:21.76",
                "22:00:21.76",
                "22:00:21.76",
                "22:00:21.76",
                "22:00:21.76",
                "22:00:21.76",
                "2:01:01.19",
            ]

            self.end_LST = [
                "22:00:21.76",
                "02:01:01.19",
                "02:01:01.19",
                "02:01:01.19",
                "02:01:01.19",
                "02:01:01.19",
                "02:01:01.19",
                "02:01:01.19",
                "02:01:01.19",
                "02:01:01.19",
                "02:01:01.19",
                "6:53:29.11",
            ]

            self.dec_min = None
            self.dec_max = None
            self.begin_RA = None
            self.end_RA = None
            self.orientation = None

            # Center of the patch in RA/Dec
            self.ra_mid = 0.0
            self.dec_mid = -57.5
        elif self.name_strategy == "small_aperture":
            self.elevation = [50.0] * 6

            self.az_min = [37.1, 282.8, 110.3, 197.0, 110.3, 197.0]

            self.az_max = [77.7, 322.9, 163.0, 250.2, 163.0, 250.2]

            self.begin_LST = [
                "08:30:00.00",
                "12:30:00.00",
                "16:30:00.00",
                "02:30:00.00",
                "16:30:00.00",
                "22:30:00.00",
            ]

            self.end_LST = [
                "12:29:59.32",
                "16:29:59.32",
                "02:29:59.30",
                "08:29:59.98",
                "22:29:59.99",
                "08:29:59.30",
            ]

            self.dec_min = None
            self.dec_max = None
            self.begin_RA = None
            self.end_RA = None
            self.orientation = None

            # Center of the patch in RA/Dec
            self.ra_mid = 0.0
            self.dec_mid = 0.0
        elif self.name_strategy == "shallow_patch":
            self.elevation = [30.0, 30.0, 30.0, 30.0, 30.0, 30.0]

            self.dec_min = [
                "-5:00:00",
                "-5:00:00",
                "-15:00:00",
                "-15:00:00",
                "-60:00:00",
                "-60:00:00",
            ]

            self.dec_max = [
                "20:00:00",
                "20:00:00",
                "20:00:00",
                "20:00:00",
                "-15:00:00",
                "-15:00:00",
            ]

            self.begin_RA = [
                "7:40:00",
                "7:40:00",
                "20:00:00",
                "20:00:00",
                "20:00:00",
                "20:00:00",
            ]

            self.end_RA = [
                "15:20:00",
                "15:20:00",
                "5:40:00",
                "5:40:00",
                "5:40:00",
                "5:40:00",
            ]

            self.orientation = ["west", "east", "west", "east", "west", "east"]

            self.az_min = None
            self.az_max = None
            self.begin_LST = None
            self.end_LST = None

            # Center of the patch in RA/Dec
            self.ra_mid = 0.0
            self.dec_mid = 0.0
        elif self.name_strategy == "custom":
            self.elevation = None
            self.dec_min = None
            self.dec_max = None
            self.begin_RA = None
            self.end_RA = None
            self.orientation = None
            self.az_min = None
            self.az_max = None
            self.begin_LST = None
            self.end_LST = None

            # Center of the patch in RA/Dec
            self.ra_mid = 0.0
            self.dec_mid = 0.0
        else:
            raise ValueError(
                """Only name_strategy = deep_patch or
                shallow_patch or custom are
                currently available. For another usage
                (advanced users), modify this routine."""
            )

    def run_one_scan(self, scan_file, scan_number):
        """
        Generate one observation (i.e. one CES) of the telescope.

        Parameters
        ----------
        scan_file : dictionary
            Empty dictionary which will contain the outputs of the scan.
        scan_number : int
            Index of the scan (between 0 and nces - 1).

        Returns
        ----------
        bool : bool
            Returns True if the scan has been generated, and False if the scan
            already exists on the disk.
        """
        # Check if we have too much/enough information to make a scan
        msg = "You cannot specify azimuth and declination!"
        assert (getattr(self, "az_min") and not getattr(self, "dec_min")) or (
            not getattr(self, "az_min") and getattr(self, "dec_min")
        ), msg
        assert (getattr(self, "az_max") and not getattr(self, "dec_max")) or (
            not getattr(self, "az_max") and getattr(self, "dec_max")
        ), msg

        msg = "You cannot specify LST and RA!"
        assert (getattr(self, "begin_LST") and not getattr(self, "begin_RA")) or (
            not getattr(self, "begin_LST") and getattr(self, "begin_RA")
        ), msg
        assert (getattr(self, "end_LST") and not getattr(self, "end_RA")) or (
            not getattr(self, "end_LST") and getattr(self, "end_RA")
        ), msg

        if getattr(self, "az_min") and getattr(self, "az_max"):
            msg = "You need to define timing bounds!"
            assert getattr(self, "begin_LST") is not None, msg
            assert getattr(self, "end_LST") is not None, msg
            azLST = True
            RADEC = False
            if self.verbose:
                print("Using azimuth and LST bounds")

        elif getattr(self, "dec_min") and getattr(self, "dec_max"):
            msg = "You need to define RA bounds!"
            assert getattr(self, "begin_RA") is not None, msg
            assert getattr(self, "end_RA") is not None, msg
            msg = "You need to define orientation of scan (east/west)!"
            assert getattr(self, "orientation") is not None, msg
            azLST = False
            RADEC = True
            if self.verbose:
                print("Using RA and Dec bounds")

        # Figure out the elevation to run the scan
        el = self.elevation[scan_number]

        # Define the sampling rate in Hz
        sampling_freq = self.sampling_freq

        #############################
        # Define geometry of the scan
        #############################
        if azLST:
            # Define geometry of the scan by figuring out the azimuth bounds
            az_mean = (self.az_min[scan_number] + self.az_max[scan_number]) * 0.5
            az_throw = (self.az_max[scan_number] - self.az_min[scan_number]) / np.cos(
                el / radToDeg
            )

        elif RADEC:
            # If given bounds in declination, make bounds in azimuth
            # note there is no sanity checking here!
            az_array = np.linspace(0.0, 180.0, endpoint=True, num=360)
            ra_array = np.zeros(az_array.shape)
            dec_array = np.zeros(az_array.shape)
            dec_min = ephem.degrees(self.dec_min[scan_number])
            dec_max = ephem.degrees(self.dec_max[scan_number])
            if self.orientation[scan_number] == "west":
                az_array += 180.0

            for i in range(0, az_array.shape[0]):
                ra_array[i], dec_array[i] = self.telescope_location.radec_of(
                    az_array[i] / radToDeg, el / radToDeg
                )

            az_allowed = np.asarray(
                [
                    az_array[i]
                    for i in range(0, az_array.shape[0])
                    if (dec_array[i] > dec_min and dec_array[i] < dec_max)
                ]
            )

            if az_allowed.shape[0] < 2:
                ms = "Invalid combination of declination bounds and elevation."
                print(ms)
                exit(1)

            az_max = np.max(az_allowed)
            az_min = np.min(az_allowed)
            az_mean = (az_min + az_max) * 0.5
            az_throw = az_max - az_min

        #############################
        # Define the timing bounds!
        #############################
        if azLST:
            LST_now = float(self.telescope_location.sidereal_time()) / (2 * np.pi)
            begin_LST = float(ephem.hours(self.begin_LST[scan_number])) / (2 * np.pi)
            end_LST = float(ephem.hours(self.end_LST[scan_number])) / (2 * np.pi)

            if begin_LST > end_LST:
                begin_LST -= 1.0

            # Reset the date to correspond to the sidereal time to start
            self.telescope_location.date -= (
                (LST_now - begin_LST) * sidDayToSec
            ) * ephem.second

            # Figure out how long to run the scan for
            num_pts = int((end_LST - begin_LST) * sidDayToSec * sampling_freq)

        if RADEC:
            self.telescope_location.horizon = el * ephem.degree

            # Define a fixed source at the ra, dec that we want to scan
            target_min_ra = ephem.FixedBody()
            target_max_ra = ephem.FixedBody()

            # Figure out where we are looking now
            ra_target, dec_target = self.telescope_location.radec_of(
                az_mean / radToDeg, el / radToDeg
            )

            # Instantiate the targets
            target_min_ra._dec = dec_target
            target_max_ra._dec = dec_target
            target_min_ra._ra = self.begin_RA[scan_number]
            target_max_ra._ra = self.end_RA[scan_number]

            # Compute initial RA
            target_min_ra.compute(self.telescope_location)
            target_max_ra.compute(self.telescope_location)
            if self.orientation[scan_number] == "east":
                self.telescope_location.date = self.telescope_location.next_rising(
                    target_min_ra
                )

                # Recompute coodinates in the light of change of date
                target_min_ra.compute(self.telescope_location)
                target_max_ra.compute(self.telescope_location)

                # Update number of time samples for the scan
                ra_max_rising = self.telescope_location.next_rising(target_max_ra)
                diff_ra = ra_max_rising - self.telescope_location.date
                num_pts = int(diff_ra / ephem.second * sampling_freq)

            if self.orientation[scan_number] == "west":
                self.telescope_location.date = self.telescope_location.next_setting(
                    target_min_ra
                )

                # Recompute coodinates in the light of change of date
                target_min_ra.compute(self.telescope_location)
                target_max_ra.compute(self.telescope_location)

                # Update number of time samples for the scan
                ra_max_setting = self.telescope_location.next_setting(target_max_ra)
                diff_ra = ra_max_setting - self.telescope_location.date
                num_pts = int(diff_ra / ephem.second * sampling_freq)

        # Run the scan!
        pb_az_dir = 1.0
        upper_az = az_mean + az_throw / 2.0
        lower_az = az_mean - az_throw / 2.0
        az_speed = self.sky_speed / np.cos(el / radToDeg)
        running_az = az_mean

        # Initialize arrays
        pb_az_array = np.zeros(num_pts)
        pb_mjd_array = np.zeros(num_pts)
        pb_ra_array = np.zeros(num_pts)
        pb_dec_array = np.zeros(num_pts)
        pb_el_array = np.ones(num_pts) * el

        # Subscan boundaries
        pb_subscans = []
        pb_subscans.append(0)

        # Loop over time samples
        # begin_lst = str(self.telescope_location.sidereal_time())
        # Pad scans 10 seconds on either side
        time_padding = 10.0 / 86400.0

        # Start of the scan
        pb_az_array[0] = running_az
        pb_mjd_array[0] = date_to_mjd(self.telescope_location.date)

        # Initialize the time
        scan_file["firstmjd"] = pb_mjd_array[0] - time_padding

        # Update before starting the loop
        running_az += az_speed * pb_az_dir / sampling_freq
        self.telescope_location.date += ephem.second / sampling_freq

        if self.language == "python":
            for t in range(0, num_pts):
                # Set the Azimuth and time
                pb_az_array[t] = running_az

                pb_ra_array[t], pb_dec_array[t] = self.telescope_location.radec_of(
                    pb_az_array[t] * np.pi / 180.0, pb_el_array[t] * np.pi / 180.0
                )

                # Case to change the direction of the scan
                if running_az > upper_az:
                    pb_az_dir = -1.0
                    pb_subscans.append(t)
                elif running_az < lower_az:
                    pb_az_dir = 1.0
                    pb_subscans.append(t)

                running_az += az_speed * pb_az_dir / sampling_freq

                # Increment the time by one second / sampling rate
                if t > 0:
                    pb_mjd_array[t] = (
                        pb_mjd_array[t - 1] + ephem.second / sampling_freq
                    )

                # Increment the time by one second / sampling rate
                self.telescope_location.date += ephem.second / sampling_freq

        elif self.language == "fortran":
            second = 1.0 / 24.0 / 3600.0
            scanning_strategy_f.run_one_scan_f(
                pb_az_array,
                pb_mjd_array,
                running_az,
                upper_az,
                lower_az,
                az_speed,
                pb_az_dir,
                second,
                sampling_freq,
                num_pts,
            )

        # Do not use that for precision - it truncates values
        self.telescope_location.date += num_pts * ephem.second / sampling_freq

        # Save in file
        scan_file["nces"] = self.nces
        scan_file["CES"] = scan_number
        scan_file["sample_rate"] = sampling_freq
        scan_file["sky_speed"] = self.sky_speed
        scan_file["lastmjd"] = pb_mjd_array[-1] + time_padding

        scan_file["azimuth"] = pb_az_array * np.pi / 180
        scan_file["elevation"] = pb_el_array * np.pi / 180
        scan_file["clock-utc"] = pb_mjd_array

        scan_file["RA"] = pb_ra_array
        scan_file["Dec"] = pb_dec_array

        scan_file["nts"] = len(pb_mjd_array)

        scan_file["subscans"] = zip(pb_subscans[:-1], np.array(pb_subscans[1:]) - 1)

        if self.verbose:
            print("+-----------------------------------+")
            print(
                " CES starts at %s and finishes at %s"
                % (
                    mjd_to_greg(scan_file["firstmjd"]),
                    mjd_to_greg(scan_file["lastmjd"]),
                )
            )
            print(
                " It lasts %.3f hours"
                % ((scan_file["lastmjd"] - scan_file["firstmjd"]) * 24)
            )
            print("+-----------------------------------+")

        # Add one day before the next CES (to avoid conflict of time)
        self.telescope_location.date += 24 * ephem.second * 3600

        # Add the scan into the instance
        # self._update('scan{}'.format(scan_number), scan_file)

        return True

    def run(self):
        """
        Generate all the observations (i.e. all CES) of the telescope.

        Examples
        ----------
        >>> scan = ScanningStrategy(sampling_freq=1., nces=2,
        ...     name_strategy='deep_patch')
        >>> scan.run()
        >>> print(scan.scan0['firstmjd'], scan.scan0['lastmjd'])
        56293.6202546 56293.8230093

        By default, the language used for the core computation is the Python.
        It can be quite slow for heavy configuration, and one can set up
        the language to fortran for speeding up the computation (x1000).
        Note that for using fortran codes you need first to compile it.
        See the setup.py or the provided Makefile.
        >>> scan = ScanningStrategy(sampling_freq=1., nces=2,
        ...     language='fortran', name_strategy='shallow_patch')
        >>> scan.run()
        >>> print(round(scan.scan0['firstmjd'], 2))
        56293.37

        Note that you can create your own scanning strategy. First choose
        the custom ones (set everything to None):
        >>> scan = ScanningStrategy(sampling_freq=1., nces=2,
        ...     language='fortran', name_strategy='custom')

        And then define your own parameters. Example:
        >>> scan.nces = 1
        >>> scan.elevation = [30.]
        >>> scan.az_min = [60.]
        >>> scan.az_max = [100.]
        >>> scan.begin_LST = ['17:00:00.00']
        >>> scan.end_LST = ['22:00:00.00']
        >>> scan.run()

        Note that you To create a scanning strategy within s4cmb,
        you need to specify either
        * Elevations + minimum and maximum azimuths (spatial bounds) +
            begining and end Local Sidereal Times (timing bounds)
        * Elevations + minimum and maximum declinations +
            begining and end Right Ascensions (spatial & timing bounds) +
            orientations (east/west)

        """
        # Initialise the date and loop over CESes
        self.telescope_location.date = self.start_date
        for CES_position in range(self.nces):
            # Initialise the starting date of observation
            # It will be updated then automatically
            setattr(self, "scan{}".format(CES_position), {})

            # Create the scan strategy
            self.run_one_scan(
                getattr(self, "scan{}".format(CES_position)), CES_position
            )

    def visualize_my_scan(
        self,
        nside,
        reso=6.9,
        xsize=900,
        rot=[0, -57.5],
        nfid_bolometer=6000,
        fp_size=180.0,
        boost=1.0,
        fullsky=False,
        flatsky=False,
        nest=False,
    ):
        """
        Simple map-making: project time ordered data into sky maps for
        visualisation. It works only in pure python (i.e. if you set
        language='python' when initialising the scanning_strategy class).

        Parameters
        ----------
        nside : int
            Resolution of the healpix map.
        reso : float
            Resolution of the projected map (pixel size) in arcmin.
        xsize : int
            Number of pixel per row for the projected map.
        rot : 1d array of float
            Coordinates of the center of
            the projected map [lon, lat] in degree
        nfid_bolometer : int
            Fiducial number of bolometers for visualisation.
            By default, the scans are done for one reference bolometer.
        fp_size : float
            Size of the focal plane on the sky, in arcmin.
        boost : int
            Boost factor to artificially increase the number of hits.
            It doesn't change the shape of the survey (just the amplitude).
        fullsky : boolean
            If True, visualising full sky. If False, visualising area defined by
            reso and my xsize.
        flatsky : boolean
            If True, visualising projection in a square array with the specified
            resolution. If False, visualising using healpy functions.
        nest : boolean
            If flatsky projection, using the defined nesting in projection.

        Outputs
        ----------
        nhit: Sky map with cumulative hit counts.
            If flatsky is True, a square numpy array of size xsize x xsize is
            returned. Else, a 1D healpy array is returned.

        Examples
        ----------
        >>> import matplotlib
        >>> matplotlib.use("Agg") # remove this if you want to display
        >>> scan = ScanningStrategy(sampling_freq=1., nces=1)
        >>> scan.run()
        >>> nhit = scan.visualize_my_scan(512)

        """
        import pylab as pl

        if self.language != "python":
            raise ValueError(
                """Visualisation is available only in pure
                python because we do not provide (yet)
                RA and Dec in fortran. Relaunch
                using language='python' in the class
                ScanningStrategy."""
            )

        npix = hp.pixelfunc.nside2npix(nside)
        nhit = np.zeros(npix)
        for scan_number in range(self.nces):
            scan = getattr(self, "scan{}".format(scan_number))

            num_pts = len(scan["clock-utc"])
            pix_global = hp.pixelfunc.ang2pix(
                nside, (np.pi / 2.0) - scan["Dec"], scan["RA"]
            )

            # Boresight pointing healpix maps
            nhit_loc = np.zeros(npix, dtype=np.int32)
            scanning_strategy_f.mapmaking(pix_global, nhit_loc, npix, num_pts)

            # Fake large focal plane with many bolometers for visualisation.
            nhit_loc = convolve_focalplane(nhit_loc, nfid_bolometer, fp_size, boost)

            nhit += nhit_loc

        if self.verbose:
            print(
                "Stats: nhits = {}/{} (fsky={}%), max hit = {}".format(
                    len(nhit[nhit > 0]),
                    len(nhit),
                    round(len(nhit[nhit > 0]) / len(nhit) * 100, 2),
                    int(np.max(nhit)),
                )
            )

        title = """
        nbolos = {},
        fp size = {} arcmin,
        nhit boost = {}
        """.format(nfid_bolometer, fp_size, boost)
        if flatsky:
            if fullsky:
                # projecting full sky with the given resolution onto a squared array.
                A_sky_deg2 = 360.0 ** 2 / np.pi
                N_pix = (
                    60.0 / reso * np.sqrt(A_sky_deg2)
                )  # number of pixels on a side
                N_pix = np.ceil(N_pix)  # rounding up, integer number required
                #########
                # for full sky:
                lonra = [-180, 180]
                latra = [-90, 90]
                #########
                hpcp = hp.projector.CartesianProj(
                    xsize=N_pix, ysize=N_pix, lonra=lonra, latra=latra
                )
                flat_hits = np.flipud(
                    hpcp.projmap(
                        nhit,
                        lambda x, y, z:
                            hp.pixelfunc.vec2pix(nside, x, y, z, nest=nest)
                    )
                )  # flipping array just for aesthetics
                # plot
                pl.imshow(flat_hits, cmap=pl.cm.viridis)
                pl.colorbar()
                title_ = """
                npix = {},
                nbolos = {},
                fp size = {} arcmin,
                nhit boost = {}
                """.format(int(N_pix), nfid_bolometer, fp_size, boost)
                pl.title(title_)
                nhit = flat_hits  # will return a flat hits map
            else:
                nhit[nhit == 0] = hp.UNSEEN
                nhit = hp.gnomview(
                    nhit,
                    rot=rot,
                    reso=reso,
                    xsize=xsize,
                    cmap=pl.cm.viridis,
                    title=title,
                    return_projected_map=True,
                )  # will return a flat hits map
                hp.graticule(verbose=self.verbose)
        else:
            if fullsky:
                nhit[nhit == 0] = hp.UNSEEN
                hp.mollview(
                    nhit,
                    rot=rot,
                    cmap=pl.cm.viridis,
                    title=title
                )
                hp.graticule(verbose=self.verbose)
            else:
                nhit[nhit == 0] = hp.UNSEEN
                hp.gnomview(
                    nhit,
                    rot=rot,
                    reso=reso,
                    xsize=xsize,
                    cmap=pl.cm.viridis,
                    title=title,
                )
                hp.graticule(verbose=self.verbose)

        pl.show()

        return nhit

    def _update(self, name, value):
        """
        Wrapper around setattr function.
        Set a named attribute on an object.

        Parameters
        ----------
        name : string
            The name of the new attribute
        value : *
            The value of the attribute.
        """
        setattr(self, name, value) class TimeOrderedDataPairDiff:
    """ Class to handle Time-Ordered Data (TOD) """

    def __init__(
        self,
        hardware,
        scanning_strategy,
        HealpixFitsMap,
        CESnumber,
        projection="healpix",
        nside_out=None,
        pixel_size=None,
        width=140.0,
        cut_pixels_outside=True,
        array_noise_level=None,
        array_noise_seed=487587,
        array_noise_level2=None,
        array_noise_seed2=56736,
        nclouds=None,
        corrlength=None,
        alpha=None,
        f0=None,
        amp_atm=None,
        mapping_perpair=False,
        mode="standard",
        store_pointing_matrix_input=False,
        verbose=False,
        perturb_el=False,
        perturb_az=False,
        seed_pointing=0.0,
        mu_pointing=0.0,
        sigma_pointing=13.0,
        perturb_pol_angs=False,
        seed_pa=0.0,
        pa_mu=0.0,
        pa_sig=0.0,
        pa_diff=False,
    ):
        """
        C'est parti!

        Parameters
        ----------
        hardware : Hardware instance
            Instance of Hardware containing instrument parameters and models.
        scanning_strategy : ScanningStrategy instance
            Instance of ScanningStrategy containing scan parameters.
        HealpixFitsMap : HealpixFitsMap instance
            Instance of HealpixFitsMap containing input sky parameters.
        CESnumber : int
            Number of the scan to simulate. Must be between 0 and
            scanning_strategy.nces - 1.
        projection : string, optional
            Type of projection for the output map. Currently available:
            healpix, flat. Here is a warning: Because of projection artifact,
            if you choose flat projection, then we will scan the sky *as if
            it was centered in [0., 0.]*. Therefore, one cannot for the moment
            compared directly healpix and flat runs.
        nside_out : int, optional
            The resolution for the output maps if projection=healpix.
            Default is nside of the input map.
        pixel_size : float, optional
            The pixel size for the output maps if projection=flat.
            In arcmin. Default is resolution of the input map.
        width : float, optional
            Width for the output map in degree.
        cut_pixels_outside : bool, optional
            The mapping time -> map is done relatively to the sky patch
            defined by (width, ra_src, dec_src). So if for some reason your
            scanning strategy goes outside the defined sky patch, the routine
            will assign -1 to the pixel index
            (or crash if cut_pixels_outside is False).
            Long story short: make sure that (width, ra_src, dec_src) returns a
            sky patch bigger than what has been defined in the
            scanning strategy, or you will have truncated output sky maps.
            Why this is so? To save memory by not having full sky maps stored
            in the memory and in the same time be able to coadd maps from
            different scans easily. If you have a better idea, let me know!
        array_noise_level : float, optional
            Noise level for the whole array in [u]K.sqrt(s). If not None, it
            will inject on-the-fly noise in time-domain while scanning
            the input map (map2tod). WARNING: units has to be same as the
            input map! Note also that it corresponds to the polarisation level.
        array_noise_seed : int, optional
            Seed used to generate random numbers to simulate noise.
            From this single seed, we generate a list of seeds
            for all detectors. Has an effect only if array_noise_level is
            provided.
        nclouds : int, optional
            Number of clouds, that is number of correlated regions
            in the focal plane.
        corrlength : int, optional
            Correlation length in time over which atmosphere signal
            is supposed to be constant. Few minutes is typical.
            Units are seconds.
        alpha : int or float, optional
            Value of the 1/f slope for correlated noise simulation.
        f0 : float, optional
            Knee frequency in Hz for correlated noise simulation.
        amp_atm : float, optional
            Typical value of atmospheric fluctuations in [u]K^2.s.
            WARNING: units has to be same as the input map!
        mapping_perpair : bool, optional
            If True, assume that you want to process pairs of bolometers
            one-by-one, that is pairs are uncorrelated. Default is False (and
            should be False unless you know what you are doing).
        mode : string, optional
            Choose between `standard` (1 frequency band) and `dichroic`
            (2 frequency bands). If `dichroic` is chosen, make sure your
            hardware can handle it (see instrument.py) and HealpixFitsMap
            should contain the inputs maps at different frequency.
        store_pointing_matrix_input : bool, optional
            Store pointing matrix used to scan the input map
        perturb_el : boolean, optional
            If True, perturb one of a few pointing instances.
            If False, sets error arrays to 0 so that the additional
            pointing instances are the same as the original ones.
        perturb_az : boolean, optional
            If True, perturb one of a few pointing instances.
            If False, sets error arrays to 0 so that the additional
            pointing instances are the same as the original ones.
        seed_pointing : integer, optional
            Seed for boresight pointing perturbation. Doesn't affect anything
            if perturb_el == False or perturb_az == False.
        mu_pointing : integer, optional
            Mean for boresight error distribution in arcseconds.
            Doesn't affect anything perturb_el == False or perturb_az == False.
        sigma_pointing : integer, optional
            Width for boresight error distribution in arcseconds.
            Doesn't affect anything perturb_el == False or perturb_az == False.
        perturb_pol_angs : bool, optional
            If True, perturbing polarization angles.
        seed_pa : int, optional
            seed for perturbing polarization angle
            (ch is added to it for randomization)
        pa_mu : float, optional
            mean of pa syst perturbation.
        pa_sig : float, optional
            width of pa syst perturbation.
        pa_diff : bool, optional
            if True, doing differential pa perturbation
        """
        # Initialise args
        self.verbose = verbose
        self.hardware = hardware
        self.scanning_strategy = scanning_strategy
        self.HealpixFitsMap = HealpixFitsMap
        self.mapping_perpair = mapping_perpair

        # pol ang syst:
        self.perturb_pol_angs = perturb_pol_angs
        self.seed_pa = seed_pa
        self.pa_mu = pa_mu
        self.pa_sig = pa_sig
        self.pa_diff = pa_diff

        self.perturb_el = perturb_el
        self.perturb_az = perturb_az

        # Check if you can run dichroic detectors
        self.mode = mode
        if self.mode == "dichroic" and (
            not hasattr(self.HealpixFitsMap, "I2")
        ):
            raise IOError("You need two sets of maps for dichroic detectors!")
        if self.mode == "dichroic" and (
            not hasattr(self.hardware, "focal_plane2")
        ):
            raise IOError("You need two sets of det for dichroic detectors!")
        if self.mode == "dichroic" and (
            not hasattr(self.hardware, "beam_model2")
        ):
            raise IOError("You need two sets of det for dichroic detectors!")

        self.width = width
        self.cut_pixels_outside = cut_pixels_outside
        self.store_pointing_matrix_input = store_pointing_matrix_input

        # Check the projection
        self.projection = projection
        assert self.projection in ["healpix", "flat"], ValueError(
            """Projection <{}> for
            the output map not understood!
            Choose among ['healpix', 'flat']."""
        ).format(self.projection)

        # Check if the number of CES is consistent with the scanning strategy
        self.CESnumber = CESnumber
        assert self.CESnumber < self.scanning_strategy.nces, ValueError(
            "The scan index must be between 0 and {}.".format(
                self.scanning_strategy.nces - 1
            )
        )

        # Initialise internal parameters
        self.scan = getattr(
            self.scanning_strategy, "scan{}".format(self.CESnumber)
        )
        self.nsamples = self.scan["nts"]
        self.npair = self.hardware.focal_plane.npair
        self.pair_list = np.reshape(
            self.hardware.focal_plane.bolo_index_in_fp, (self.npair, 2)
        )

        # boresight pointing perturbation errors
        if self.perturb_az or self.perturb_el:
            mu_pointing = mu_pointing / 3600 * np.pi / 180  # arcsec to radians
            # arcsec to degrees and also dividing by sqrt(2) so that
            # the overal error on position is np.sqrt(\delta_az^2+\delta_el^2)
            sigma_pointing = (
                sigma_pointing / (3600 * np.sqrt(2)) * np.pi / 180
            )
        if self.perturb_az:
            # perturb azimuth
            seed_pointing_az = seed_pointing
            state_for_pointing_errors_az = np.random.RandomState(
                seed_pointing_az
            )
            self.err_azimuth = state_for_pointing_errors_az.normal(
                mu_pointing, sigma_pointing, self.scan["nts"]
            )
        else:
            self.err_azimuth = np.zeros(self.scan["nts"])

        if self.perturb_el:
            # perturb elevation
            seed_pointing_el = (
                seed_pointing + 1234567890
            )
            # different seed for each perturbation;
            # tried to avoid having another seed input
            state_for_pointing_errors_el = np.random.RandomState(
                seed_pointing_el
            )
            self.err_elevation = state_for_pointing_errors_el.normal(
                mu_pointing, sigma_pointing, self.scan["nts"]
            )
        else:
            self.err_elevation = np.zeros(self.scan["nts"])

        # Pre-compute boresight pointing objects
        self.get_boresightpointing()

        # Polarisation angles: intrinsic and HWP angles
        self.get_angles()

        # Position of bolometers in the focal plane
        # TODO move that elsewhere...
        self.ypos = self.hardware.beam_model.ypos
        self.xpos = self.hardware.beam_model.xpos
        self.xpos = self.xpos / np.cos(self.ypos)

        # Initialise pointing matrix, that is the matrix to go from time
        # to map domain, for all pairs of detectors.
        if not self.mapping_perpair:
            self.point_matrix = np.zeros(
                (self.npair, self.nsamples), dtype=np.int32
            )

        else:
            self.point_matrix = np.zeros((1, self.nsamples), dtype=np.int32)

        # If set, stores the pointing matrix used to scan the map
        if self.store_pointing_matrix_input:
            self.point_matrix_input = np.zeros(
                self.point_matrix.shape, dtype=np.int32
            )

        # Initialise the mask for timestreams
        self.wafermask_pixel = self.get_timestream_masks()

        # Boundaries for subscans (t_beg, t_end)
        self.subscans = self.scan["subscans"]

        # Get observed pixels in the input map
        if nside_out is None:
            self.nside_out = self.HealpixFitsMap.nside
        else:
            self.nside_out = nside_out
        if pixel_size is None:
            self.pixel_size = hp.nside2resol(self.HealpixFitsMap.nside)
        else:
            self.pixel_size = pixel_size * am2rad

        self.obspix, self.npixsky = self.get_obspix(
            self.width,
            self.scanning_strategy.ra_mid,
            self.scanning_strategy.dec_mid,
        )

        # Get timestream weights
        self.sum_weight, self.diff_weight = self.get_weights()

        # Get detector gains
        self.set_detector_gains()

        # Prepare noise simulator if needed
        self.array_noise_level = array_noise_level
        self.array_noise_seed = array_noise_seed
        self.nclouds = nclouds
        self.corrlength = corrlength
        self.alpha = alpha
        self.f0 = f0
        self.amp_atm = amp_atm
        if self.array_noise_level is not None and self.alpha is None:
            self.noise_generator = WhiteNoiseGenerator(
                array_noise_level=self.array_noise_level,
                ndetectors=2 * self.npair,
                ntimesamples=self.nsamples,
                array_noise_seed=self.array_noise_seed,
                sampling_freq=self.scanning_strategy.sampling_freq,
            )
        elif self.array_noise_level is not None and self.alpha is not None:
            self.noise_generator = CorrNoiseGenerator(
                array_noise_level=self.array_noise_level,
                ndetectors=2 * self.npair,
                ntimesamples=self.nsamples,
                array_noise_seed=self.array_noise_seed,
                sampling_freq=self.scanning_strategy.sampling_freq,
                nclouds=self.nclouds,
                f0=self.f0,
                amp_atm=self.amp_atm,
                corrlength=self.corrlength,
                alpha=self.alpha,
            )

        else:
            self.noise_generator = None

        self.array_noise_level2 = array_noise_level2
        self.array_noise_seed2 = array_noise_seed2
        if self.array_noise_level2 is not None and self.alpha is None:
            self.noise_generator2 = WhiteNoiseGenerator(
                array_noise_level=self.array_noise_level2,
                ndetectors=2 * self.npair,
                ntimesamples=self.nsamples,
                array_noise_seed=self.array_noise_seed2,
                sampling_freq=self.scanning_strategy.sampling_freq,
            )

        elif self.array_noise_level2 is not None and self.alpha is not None:
            self.noise_generator2 = CorrNoiseGenerator(
                array_noise_level=self.array_noise_level2,
                ndetectors=2 * self.npair,
                ntimesamples=self.nsamples,
                array_noise_seed=self.array_noise_seed2,
                sampling_freq=self.scanning_strategy.sampling_freq,
                nclouds=self.nclouds,
                f0=self.f0,
                amp_atm=self.amp_atm,
                corrlength=self.corrlength,
                alpha=self.alpha,
            )

        else:
            self.noise_generator2 = None

    def get_angles(self):
        """
        Retrieve polarisation angles: intrinsic (focal plane) and HWP angles,
        and initialise total polarisation angle.
        """
        self.hwpangle = self.hardware.half_wave_plate.compute_HWP_angles(
            sample_rate=self.scan["sample_rate"], size=self.nsamples
        )

        self.intrinsic_polangle = self.hardware.focal_plane.bolo_polangle
        self.intrinsic_polangle2 = None
        if self.mode == "dichroic":
            self.intrinsic_polangle2 = self.hardware.focal_plane2.bolo_polangle

        # Will contain the total polarisation angles for all bolometers
        # That is PA + intrinsic + 2 * HWP
        if not self.mapping_perpair:
            self.pol_angs = np.zeros((self.npair, self.nsamples))
            self.pol_angs2 = None
            if self.mode == "dichroic":
                self.pol_angs2 = np.zeros((self.npair, self.nsamples))
        else:
            self.pol_angs = np.zeros((1, self.nsamples))
            self.pol_angs2 = None
            if self.mode == "dichroic":
                self.pol_angs2 = np.zeros((1, self.nsamples))

    def get_timestream_masks(self):
        """
        Define the masks for all the timestreams.
        1 if the time sample should be included, 0 otherwise.
        Set to ones for the moment.
        """
        if not self.mapping_perpair:
            return np.ones((self.npair, self.nsamples), dtype=int)
        else:
            return np.ones((1, self.nsamples), dtype=int)

    def get_obspix(self, width, ra_src, dec_src):
        """
        Return the index of observed pixels within a given patch
        defined by (`ra_src`, `dec_src`) and `width`.
        This will be the sky patch that will be returned.

        Parameters
        ----------
        width : float
            Width of the patch in degree.
        ra_src : float
            RA of the center of the patch in degree.
        dec_src : float
            Dec of the center of the patch in degree.

        Returns
        ----------
        obspix : 1d array
            The indices of the observed pixels in the input map. Same for
            all output projection as the only input projection is healpix.
        nskypix : int
            Number of sky pixels in our patch given the boundaries (specified
            via `width`). For healpix projection, this is the length of the
            obspix. For flat projection, this is the number of square pixels
            within the patch.

        Examples
        ----------
        Healpix projection
        >>> inst, scan, sky_in = load_fake_instrument()
        >>> tod = TimeOrderedDataPairDiff(inst, scan, sky_in, CESnumber=0)
        >>> obspix, npix = tod.get_obspix(10., 0., 0.)
        >>> print(obspix)
        [1376 1439 1440 1504 1567 1568 1632 1695]

        Flat projection
        >>> inst, scan, sky_in = load_fake_instrument()
        >>> tod = TimeOrderedDataPairDiff(inst, scan, sky_in, CESnumber=0,
        ...     projection='flat')
        >>> obspix, npix = tod.get_obspix(10., 0., 0.)
        >>> print(obspix) # the same as healpix
        [1376 1439 1440 1504 1567 1568 1632 1695]
        >>> print(npix, len(obspix))
        16 8
        """
        # Change to radian
        ra_src = ra_src * d2r
        dec_src = dec_src * d2r
        # TODO implement the first line correctly...
        try:
            xmin, xmax, ymin, ymax = np.array(width) * d2r
        except TypeError:
            xmin = xmax = ymin = ymax = d2r * width / 2.0

        # If map bound crosses zero make coordinates
        # bounds monotonic in [-pi,pi]
        ra_min = ra_src - xmin
        if (ra_src + xmax) >= 2 * np.pi:
            ra_max = (ra_src + xmax) % (2 * np.pi)
            ra_min = ra_min if ra_min <= np.pi else ra_min - 2 * np.pi
        else:
            ra_min = ra_src - xmin
            ra_max = ra_src + xmax

        dec_min = max([(dec_src - ymin), -np.pi / 2])
        dec_max = min([dec_src + ymax, np.pi / 2])

        self.xmin = ra_min
        self.xmax = ra_max
        self.ymin = dec_min
        self.ymax = dec_max

        obspix = input_sky.get_obspix(
            ra_min, ra_max, dec_min, dec_max, self.nside_out
        )

        if self.projection == "flat":
            npixsky = (
                int(
                    round(
                         (self.xmax - self.xmin + self.pixel_size) / self.pixel_size
                    )
                ) ** 2
            )
        elif self.projection == "healpix":
            npixsky = len(obspix)

        return obspix, npixsky

    def get_weights(self):
        """
        Return the noise weights of the sum and difference timestreams
        (in 1/noise units).
        For the moment, there is one number per pair for the whole scan.
        Typically, this can be the (mean) PSD of the timestream.

        Default for the moment is 1 (i.e. no weights).

        Returns
        ----------
        sum_weight : 1d array
            Weights for the sum of timestreams (size: npair)
        diff_weight : 1d array
            Weights for the difference of timestreams (size: npair)
        """
        if not self.mapping_perpair:
            return np.ones((2, self.npair), dtype=int)
        else:
            return np.ones((2, 1), dtype=int)

    def set_detector_gains(self, new_gains=None, new_gains2=None):
        """
        Set the gains of the detectors (unitless).
        Default is 1., that is perfectly calibrated.

        Parameters
        ----------
        new_gains : 1d array
            Array containing the gain value for all detectors
            (1 number per detector) for the 1st frequency channel.
        new_gains2 : 1d array
            Array containing the gain value for all detectors
            (1 number per detector) for the 2nd frequency channel.

        Examples
        ----------
        >>> inst, scan, sky_in = load_fake_instrument()
        >>> tod = TimeOrderedDataPairDiff(inst, scan, sky_in, CESnumber=1)
        >>> print(tod.gain[0])
        1.0

        Change the value of gains
        >>> new_gains = np.ones(2 * tod.npair) * 2.
        >>> tod.set_detector_gains(new_gains=new_gains)
        >>> print(tod.gain[0])
        2.0

        For dichroic
        >>> inst, scan, sky_in = load_fake_instrument(fwhm_in2=1.8)
        >>> tod = TimeOrderedDataPairDiff(inst, scan, sky_in,
        ...     mode='dichroic', CESnumber=1)
        >>> assert tod.gain2 is not None

        Change the value of gains for both frequency channels
        >>> new_gains = np.ones(2 * tod.npair) * 2.
        >>> new_gains2 = np.ones(2 * tod.npair) * 4.
        >>> tod.set_detector_gains(new_gains=new_gains, new_gains2=new_gains2)
        """
        if new_gains is not None:
            assert len(new_gains) == 2 * self.npair, ValueError(
                "You have to provide {} new gain values!".format(
                    2 * self.npair
                )
            )
            self.gain = new_gains
        else:
            self.gain = np.ones(2 * self.npair)

        self.gain2 = None
        if self.mode == "dichroic":
            if new_gains2 is not None:
                assert len(new_gains2) == 2 * self.npair, ValueError(
                    "You have to provide {} new gain values!".format(
                        2 * self.npair
                    )
                )
                self.gain2 = new_gains2
            else:
                self.gain2 = np.ones(2 * self.npair)

    def set_detector_gains_perpair(self, new_gains=None, new_gains2=None):
        """
        Set the gains of all 2 pair detectors for each timestep (unitless).
        This is particularly useful to introduce drifts for example.
        Default is 1., that is perfectly calibrated.

        Parameters
        ----------
        new_gains : 2d array of size (2, nsamples)
            Array containing the gain value for all detectors
            (nsamples number per detector) for the 1st frequency channel.
        new_gains2 : 2d array of size (2, nsamples)
            Array containing the gain value for all detectors
            (nsamples number per detector) for the 2nd frequency channel.

        Examples
        ----------
        >>> inst, scan, sky_in = load_fake_instrument()
        >>> tod = TimeOrderedDataPairDiff(inst, scan, sky_in, CESnumber=1)
        >>> print(tod.gain[0])
        1.0

        Change the value of gains every other sample.
        >>> new_gains = np.ones((2, tod.nsamples))
        >>> new_gains[:, ::2] = 2.
        >>> tod.set_detector_gains_perpair(new_gains=new_gains)
        >>> print(tod.gain[0][0:4])
        [ 2.  1.  2.  1.]

        Dichroic
        >>> inst, scan, sky_in = load_fake_instrument(fwhm_in2=1.8)
        >>> tod = TimeOrderedDataPairDiff(inst, scan, sky_in,
        ...     mode='dichroic', CESnumber=1)

        Change the value of gains every other sample.
        >>> new_gains = np.ones((2, tod.nsamples))
        >>> new_gains[:, ::2] = 2.
        >>> new_gains2 = np.ones((2, tod.nsamples))
        >>> new_gains2[:, ::2] = 4.
        >>> tod.set_detector_gains_perpair(
        ...     new_gains=new_gains, new_gains2=new_gains2)
        """
        if new_gains is not None:
            msg = "You have to provide ({}, {}) new gain values!"
            assert new_gains.shape == (2, self.nsamples), ValueError(
                msg.format(2, self.nsamples)
            )
            self.gain = new_gains
        else:
            self.gain = np.ones((2, self.nsamples))

        self.gain2 = None
        if self.mode == "dichroic":
            if new_gains2 is not None:
                msg = "You have to provide ({}, {}) new gain values!"
                assert new_gains2.shape == (2, self.nsamples), ValueError(
                    msg.format(2, self.nsamples)
                )
                self.gain2 = new_gains2
            else:
                self.gain2 = np.ones((2, self.nsamples))

    def get_boresightpointing(self):
        """
        Initialise the boresight pointing for all the focal plane bolometers.
        The actual pointing (RA/Dec/Parallactic angle) is computed on-the-fly
        when we load the data.

        Note:
        For healpix projection, our (ra_src, dec_src) = (0, 0) and we
        rotate the input map while for flat we true center of the patch.
        This is to avoid projection artifact by operating a rotation
        of the coordinates to (0, 0) in flat projection (scan around equator).
        """
        lat = float(self.scanning_strategy.telescope_location.lat)\
            * 180.0 / np.pi

        if self.projection == "healpix":
            ra_src = 0.0
            dec_src = 0.0
        elif self.projection == "flat":
            ra_src = self.scanning_strategy.ra_mid * np.pi / 180.0
            dec_src = self.scanning_strategy.dec_mid * np.pi / 180.0

            # Perform a rotation of the input to put the point
            # (ra_src, dec_src) at (0, 0).
            # r = hp.Rotator(rot=[ra_src, self.scanning_strategy.dec_mid])
            # theta, phi = hp.pix2ang(self.HealpixFitsMap.nside,
            #                         range(12 * self.HealpixFitsMap.nside**2))
            # t, p = r(theta, phi, inv=True)
            # pix = hp.ang2pix(self.HealpixFitsMap.nside, t, p)
            #
            # # Apply the rotation to our maps
            # self.HealpixFitsMap.I = self.HealpixFitsMap.I[pix]
            # self.HealpixFitsMap.Q = self.HealpixFitsMap.Q[pix]
            # self.HealpixFitsMap.U = self.HealpixFitsMap.U[pix]

        self.pointing = Pointing(
            az_enc=self.scan["azimuth"],  # degrees, size of nts
            el_enc=self.scan["elevation"],  # degrees, size of nts
            time=self.scan["clock-utc"],
            value_params=self.hardware.pointing_model.value_params,
            allowed_params=self.hardware.pointing_model.allowed_params,
            ut1utc_fn=self.scanning_strategy.ut1utc_fn,
            lat=lat,
            ra_src=ra_src,
            dec_src=dec_src,
        )

        # If boresight pointing systematics are present define a new pointing
        # class with a perturbed pointing to be used when creating TODs with
        # map2tod method
        if self.perturb_el or self.perturb_az:
            self.pointing_perturbed = Pointing(
                az_enc=(
                    self.scan["azimuth"] + self.err_azimuth
                ),  # degrees, size of nts
                el_enc=(
                    self.scan["elevation"] + self.err_elevation
                ),  # degrees, size of nts
                time=self.scan["clock-utc"],
                value_params=self.hardware.pointing_model.value_params,
                allowed_params=self.hardware.pointing_model.allowed_params,
                ut1utc_fn=self.scanning_strategy.ut1utc_fn,
                lat=lat,
                ra_src=ra_src,
                dec_src=dec_src,
            )
        else:
            self.pointing_perturbed = None

    def compute_simpolangle(self, ch, parallactic_angle, polangle_err=False):
        """
        Compute the full polarisation angles used to generate timestreams.
        The polarisation angle contains intrinsic polarisation angle (from
        focal plane design), parallactic angle (from pointing), and the angle
        from the half-wave plate.

        Parameters
        ----------
        ch : int
            Channel index in the focal plane.
        parallactic_angle : 1d array
            All parallactic angles for detector ch.
        polangle_err : bool, optional
            If True, inject systematic effect.
            TODO: remove that in the systematic module.

        Returns
        ----------
        pol_ang : 1d array
            Vector containing the values of the polarisation angle for the
            whole scan.

        Examples
        ----------
        >>> inst, scan, sky_in = load_fake_instrument()
        >>> tod = TimeOrderedDataPairDiff(inst, scan, sky_in, CESnumber=0)
        >>> angles, angles2 = tod.compute_simpolangle(ch=0,
        ...     parallactic_angle=np.array([np.pi] * tod.nsamples))
        >>> print(angles[:4])
        [ 0.          0.31415927  0.62831853  0.9424778 ]
        >>> assert angles2 is None

        If you have dichroic detectors, the routine will return angles for
        the two frequency channels
        >>> inst, scan, sky_in = load_fake_instrument(fwhm_in2=1.8)
        >>> tod = TimeOrderedDataPairDiff(inst, scan, sky_in,
        ...     mode='dichroic', CESnumber=0)
        >>> angles, angles2 = tod.compute_simpolangle(ch=0,
        ...     parallactic_angle=np.array([np.pi] * tod.nsamples))
        >>> assert angles2 is not None
        """

        if polangle_err:
            state_for_polang = np.random.RandomState(self.seed_pa)

            if self.pa_sig == 0:
                dpolang = (
                    np.ones(len(self.intrinsic_polangle)) * self.pa_mu
                )  # in deg
            else:
                dpolang = state_for_polang.normal(
                    self.pa_mu, self.pa_sig, len(self.intrinsic_polangle)
                )  # in deg

            if self.pa_diff:
                # differential pol ang symmetrically
                dpolang = dpolang
            else:
                dpolang[1::2] = dpolang[::2]

            # inject perturbations
            intrinsic_polangle = np.array(self.intrinsic_polangle) + dpolang
        else:
            intrinsic_polangle = self.intrinsic_polangle

        # calc total pa (with or without perturbations)

        ang_pix = (90.0 - intrinsic_polangle[ch]) * d2r

        # Demodulation or pair diff use different convention
        # for the definition of the angle.
        if not hasattr(self, "dm"):
            pol_ang = parallactic_angle + ang_pix + 2.0 * self.hwpangle
        else:
            pol_ang = parallactic_angle - ang_pix - 2.0 * self.hwpangle

        pol_ang2 = None

        if self.mode == "dichroic":
            # TODO: add support for different polarization angle perturbations
            # for each frequency of the dichroic. Assume equal perturbations
            # for the moment.
            if polangle_err:
                # inject perturbations
                intrinsic_polangle2 = np.array(self.intrinsic_polangle2) + \
                    dpolang
            else:
                intrinsic_polangle2 = self.intrinsic_polangle2
            ang_pix2 = (90.0 - intrinsic_polangle2[ch]) * d2r

            # Demodulation or pair diff use different convention
            # for the definition of the angle.
            if not hasattr(self, "dm"):
                pol_ang2 = parallactic_angle + ang_pix2 + 2.0 * self.hwpangle
            else:
                pol_ang2 = parallactic_angle - ang_pix2 - 2.0 * self.hwpangle

        return pol_ang, pol_ang2

    def return_parallactic_angle(self, ch, frequency_channel=1):
        """
        Return the parallactic angles (orientations of the pixel on sky)
        for one specific channel (full timestream).

        Parameters
        ----------
        ch : int
            Index of the bolometer (within the focal plane).
        frequency_channel : int, optional
            When using dichroic detectors, you need to specify the channel.

        Returns
        ----------
        pa : 1d array
            Parallactic angles (timestream) for detector ch. [radian]

        Examples
        ----------
        Simple TES bolometers
        >>> inst, scan, sky_in = load_fake_instrument()
        >>> tod = TimeOrderedDataPairDiff(inst, scan, sky_in, CESnumber=1)
        >>> pa = tod.return_parallactic_angle(0)

        Dichroic detectors
        >>> inst, scan, sky_in = load_fake_instrument(fwhm_in2=1.8)
        >>> tod = TimeOrderedDataPairDiff(inst, scan, sky_in,
        ...     mode='dichroic', CESnumber=1)
        >>> pa1 = tod.return_parallactic_angle(0, frequency_channel=1)
        >>> pa2 = tod.return_parallactic_angle(0, frequency_channel=2)
        """
        if frequency_channel == 1:
            pa = self.pol_angs
        elif frequency_channel == 2:
            pa = self.pol_angs2

        if self.mapping_perpair is True:
            if frequency_channel == 1:
                ang_pix = (90.0 - self.intrinsic_polangle[ch]) * d2r
            elif frequency_channel == 2:
                ang_pix = (90.0 - self.intrinsic_polangle2[ch]) * d2r

            if not hasattr(self, "dm"):
                return pa - ang_pix - 2 * self.hwpangle
            else:
                return pa + ang_pix
        else:
            if frequency_channel == 1:
                ang_pix = (90.0 - self.intrinsic_polangle[2 * ch]) * d2r
            elif frequency_channel == 2:
                ang_pix = (90.0 - self.intrinsic_polangle2[2 * ch]) * d2r
            if not hasattr(self, "dm"):
                return pa[ch] - ang_pix - 2 * self.hwpangle
            else:
                return pa[ch] + ang_pix

    def get_pixel_indices(self, ra, dec):
        # Retrieve corresponding pixels on the sky, and their index locally.
        if self.projection == "flat":
            # ??
            xmin = -self.width / 2.0 * np.pi / 180.0
            ymin = -self.width / 2.0 * np.pi / 180.0
            index_global, index_local = build_pointing_matrix(
                ra,
                dec,
                nside_in=self.HealpixFitsMap.nside,
                nside_out=self.nside_out,
                xmin=xmin,
                ymin=ymin,
                pixel_size=self.pixel_size,
                npix_per_row=int(np.sqrt(self.npixsky)),
                projection=self.projection,
                cut_pixels_outside=self.cut_pixels_outside,
            )
        elif self.projection == "healpix":
            index_global, index_local = build_pointing_matrix(
                ra,
                dec,
                nside_in=self.HealpixFitsMap.nside,
                nside_out=self.nside_out,
                obspix=self.obspix,
                ext_map_gal=self.HealpixFitsMap.ext_map_gal,
                projection=self.projection,
                cut_pixels_outside=self.cut_pixels_outside,
            )
        return index_global, index_local

    def map2tod(self, ch, demod_ts=False):
        """
        Scan the input sky maps to generate timestream for channel ch.
        /!\ this is currently the bottleneck in computation. Need to speed
        up this routine!

        Parameters
        ----------
        ch : int
            Channel index in the focal plane.
        demod_ts : bool
            If set returns the timestream of the channel as a demodulated
            timestream from single detector observation if a CRWHP is set.
            Useful to avoid artifacts due to CRHWP demodulation. See the
            description of the demodulate_timestreams method of the
            Demodulation class for more details.
        Returns
        ----------
        ts : 1d array or array of 1d array
            The timestream for detector ch. If `self.HealpixFitsMap.do_pol` is
            True it returns intensity+polarisation, otherwise just intensity.
            If dichroic, ts = np.array([ts_freq1, ts_freq2]). If demod_ts is
            set returns a 3d array or an array of 3d array.

        Examples
        ----------
        Simple TES bolometers
        >>> inst, scan, sky_in = load_fake_instrument()
        >>> tod = TimeOrderedDataPairDiff(inst, scan, sky_in, CESnumber=1)
        >>> d = tod.map2tod(0)
        >>> print(d.shape)
        (115200,)

        Dichroic detectors
        >>> inst, scan, sky_in = load_fake_instrument(fwhm_in2=1.8)
        >>> tod = TimeOrderedDataPairDiff(inst, scan, sky_in,
        ...     mode='dichroic', CESnumber=1)
        >>> d = tod.map2tod(0)
        >>> print(d.shape)
        (2, 115200)

        Dichroic detectors with different noise levels
        >>> inst, scan, sky_in = load_fake_instrument(fwhm_in2=1.8)
        >>> tod = TimeOrderedDataPairDiff(inst, scan, sky_in,
        ...     array_noise_level=2.5, array_noise_seed=487587,
        ...     array_noise_level2=25., array_noise_seed2=56736,
        ...     mode='dichroic', CESnumber=1)
        >>> d = tod.map2tod(0)
        """
        # Use bolometer beam offsets.
        azd, eld = self.xpos[ch], self.ypos[ch]

        # Compute pointing for detector ch
        ra, dec, pa = self.pointing.offset_detector(azd, eld)

        # Compute pointing matrix for baseline pointing
        index_global, index_local = self.get_pixel_indices(ra, dec)

        # Perturbed boresight pointing values if requested
        if self.pointing_perturbed is not None:
            ra_perturbed, dec_perturbed, pa_in = \
                self.pointing_perturbed.offset_detector(azd, eld)

            # When pointing is perturbed the pointing matrix and pa used to
            # scan the sky are different from those used in the sky
            # reconstruction map2tod.
            # Dedicated variables are defined.
            index_global_in, index_local_in = \
                self.get_pixel_indices(ra_perturbed, dec_perturbed)

        else:
            # Otherwise the pointing matrix and pa used to scan the sky will be
            # the same used in the sky map reconstruction.
            pa_in = pa
            index_global_in = index_global
            index_local_in = index_local

        if ((self.projection == 'healpix') & (index_local is None)):
            # Using a pointer not to increase memory usage
            index_local = index_global
            # Redefines pointing value when pointing is perturbed
            if index_local_in is None:
                index_local_in = index_global_in

        # For flat projection, one needs to flip the sign of U
        # w.r.t to the full-sky basis (IAU angle convention)
        if self.projection == "flat":
            sign = -1.0
        elif self.projection == "healpix":
            sign = 1.0
        self.Usign = sign

        # Store pointing matrix and parallactic angle for tod2map operations
        if ch % 2 == 0:
            if (self.xpos[ch] != self.xpos[ch + 1]) or (
                self.ypos[ch] != self.ypos[ch + 1]
            ):
                # Compute pointing matrix for the pair center if beam
                # offsets for top and bottom detector do not coincide.
                azd = 0.5 * (self.xpos[ch] + self.xpos[ch + 1])
                eld = 0.5 * (self.ypos[ch] + self.ypos[ch + 1])
                ra, dec, pa_pair = self.pointing.offset_detector(azd, eld)
                index_global_pair, index_local = self.get_pixel_indices(
                    ra, dec
                )
                if index_local is None:
                    # Using a pointer not to increase memory usage
                    index_local = index_global_pair
                if self.store_pointing_matrix_input:
                    if self.pointing_perturbed is not None:
                        ra, dec, pa_pair_in = \
                            self.pointing_perturbed.offset_detector(azd, eld)
                        index_global_pair_in, index_local_pair_in = \
                            self.get_pixel_indices(ra, dec)
                        if index_local_pair_in is None:
                            index_local_pair_in = index_global_pair_in
                    else:
                        index_global_pair_in = index_global_pair
                        index_local_pair_in = index_local

            else:
                # For consistency when storing input_map scanning after
                # introducing differential beam systematics.
                index_global_pair = index_global

            # Store list of pixels to be mapped only for top bolometers or pair
            # center in presence of differential pointing
            if ch % 2 == 0 and not self.mapping_perpair:
                self.point_matrix[int(ch / 2)] = index_local
                if self.store_pointing_matrix_input:
                    index_global_pair_in[index_local_pair_in == -1] = -1
                    self.point_matrix_input[int(ch / 2)] = index_global_pair_in
            elif ch % 2 == 0 and self.mapping_perpair:
                self.point_matrix[0] = index_local
                if self.store_pointing_matrix_input:
                    index_global_pair_in[index_local_pair_in == -1] = -1
                    self.point_matrix_input[0] = index_global_pair_in

        # Default gain for a detector is 1.,
        # but you can change it using set_detector_gains or
        # set_detector_gains_perpair.
        if len(self.gain) == self.npair * 2:
            norm = self.gain[ch]
        elif len(self.gain) == 2:
            # print("Gains per pair")
            if ch % 2 == 0:
                norm = self.gain[0]
            else:
                norm = self.gain[1]

        # Noise simulation
        if self.noise_generator is not None:
            noise = self.noise_generator.simulate_noise_one_detector(ch)
        else:
            noise = 0.0

        if self.noise_generator2 is not None:
            noise2 = self.noise_generator2.simulate_noise_one_detector(ch)
        else:
            noise2 = 0.0

        if self.HealpixFitsMap.do_pol:
            # Compute polarization angle per detector used to scan the input
            # sky (accounts thus pointing errors (boresight and differential)
            # and pol. angle perturbations).
            # WARNING: pol_ang2 is None if mode == 'standard'
            pol_ang_in, pol_ang2_in = self.compute_simpolangle(
                ch, pa_in, polangle_err=self.perturb_pol_angs)

            cos2pol_ang_in = np.cos(2 * pol_ang_in)
            sin2pol_ang_in = np.sin(2 * pol_ang_in)

            # Use pa of pair center for tod2map if differential pointing
            # is used (and thus the pa_pair is defined). Throws an exception if
            # this is not the case and uses the same angle for the
            # the top and bottom detectors.
            # In this case, we use the pointing objects specifically defined
            # for the sky map reconstruction (i.e that do not include boresight
            # and differential pointing errors).
            try:
                pol_ang_pair, pol_ang2_pair = self.compute_simpolangle(
                    ch, pa_pair, polangle_err=False)
            except UnboundLocalError:
                # WARNING: pol_ang2 is None if mode == 'standard'
                pol_ang, pol_ang2 = self.compute_simpolangle(
                    ch, pa, polangle_err=False)
                pol_ang_pair = pol_ang
                pol_ang2_pair = pol_ang2

            # Defines values of pol. angle to be used in sky map
            # reconstruction.
            # For demodulation, HWP angles are not included at the level
            # of the pointing matrix (convention).
            if hasattr(self, "dm"):
                pol_ang_out = pol_ang_pair + 2.0 * self.hwpangle
            else:
                pol_ang_out = pol_ang_pair

            if (hasattr(self, "dm") and demod_ts):
                # Store list of polangle only for tod2map operations.
                # Stores data only for top bolometers without
                # additional hwp rotation in the pointing matrix required
                # when performing demodulation.
                if ch % 2 == 0 and not self.mapping_perpair:
                    self.pol_angs[int(ch / 2)] = pol_ang_pair
                elif ch % 2 == 0 and self.mapping_perpair:
                    self.pol_angs[0] = pol_ang_pair
                nt = self.nsamples
                # defines perfectly demodulated timestreams
                ts1 = np.zeros((3, nt))
                ts1[0] = self.HealpixFitsMap.I[index_global_in]
                ts1[1] = (cos2pol_ang_in * self.HealpixFitsMap.Q[index_global_in])\
                    + (sign * sin2pol_ang_in * self.HealpixFitsMap.U[index_global_in])
                ts1[2] = (sin2pol_ang_in * self.HealpixFitsMap.Q[index_global_in])\
                    - (sign * cos2pol_ang_in * self.HealpixFitsMap.U[index_global_in])
            else:
                # Store list of polangle only for top bolometers
                if ch % 2 == 0 and not self.mapping_perpair:
                    self.pol_angs[int(ch / 2)] = pol_ang_out
                elif ch % 2 == 0 and self.mapping_perpair:
                    self.pol_angs[0] = pol_ang_out
                ts1 = self.HealpixFitsMap.I[index_global_in] \
                    + self.HealpixFitsMap.Q[index_global_in] \
                    * cos2pol_ang_in \
                    + sign \
                    * self.HealpixFitsMap.U[index_global_in] \
                    * sin2pol_ang_in \
                    + noise
            ts1 = ts1 * norm

            if self.mode == "standard":
                return ts1

            elif self.mode == "dichroic":

                cos2pol_ang2_in = np.cos(2 * pol_ang2_in)
                sin2pol_ang2_in = np.sin(2 * pol_ang2_in)

                # For demodulation, HWP angles are not included at the level
                # of the pointing matrix (convention).
                if hasattr(self, "dm"):
                    pol_ang_out2 = pol_ang2_pair + 2.0 * self.hwpangle
                else:
                    pol_ang_out2 = pol_ang2_pair

                if (hasattr(self, "dm") and demod_ts):
                    # Store list of polangle only for tod2map operations.
                    # Stores data only for top bolometers without
                    # additional hwp rotation in the pointing matrix required
                    # when performing demodulation.
                    if ch % 2 == 0 and not self.mapping_perpair:
                        self.pol_angs2[int(ch / 2)] = pol_ang2_pair
                    elif ch % 2 == 0 and self.mapping_perpair:
                        self.pol_angs2[0] = pol_ang2_pair
                    nt = self.nsamples
                    # defines perfectly demodulated timestreams
                    ts2 = np.zeros((3, nt))
                    ts2[0] = self.HealpixFitsMap.I2[index_global_in]
                    ts2[1] = cos2pol_ang2_in \
                        * self.HealpixFitsMap.Q2[index_global_in]\
                        + sign * sin2pol_ang2_in\
                        * self.HealpixFitsMap.U2[index_global_in]
                    ts2[2] = sin2pol_ang2_in \
                        * self.HealpixFitsMap.Q2[index_global_in]\
                        - sign * cos2pol_ang2_in\
                        * self.HealpixFitsMap.U2[index_global_in]
                else:
                    # Store list polangle only for top bolometers
                    if ch % 2 == 0 and not self.mapping_perpair:
                        self.pol_angs2[int(ch / 2)] = pol_ang_out2
                    elif ch % 2 == 0 and self.mapping_perpair:
                        self.pol_angs2[0] = pol_ang_out2

                    ts2 = self.HealpixFitsMap.I2[index_global_in]\
                        + self.HealpixFitsMap.Q2[index_global_in]\
                        * cos2pol_ang2_in\
                        + sign\
                        * self.HealpixFitsMap.U2[index_global_in]\
                        * sin2pol_ang2_in\
                        + noise2
                    ts2 = ts2 * norm
                return np.array([ts1, ts2])

        else:
            ts1 = norm * (
                self.HealpixFitsMap.I[index_global_in] + noise
            )
            if self.mode == "standard":
                return ts1
            elif self.mode == "dichroic":
                ts2 = norm * (
                    self.HealpixFitsMap.I2[index_global_in] + noise2
                )
                return np.array([ts1, ts2])

    def tod2map(
        self, waferts, output_maps, gdeprojection=False, frequency_channel=1
    ):
        """
        Project time-ordered data into sky maps for the whole array.
        Maps are updated on-the-fly. Massive speed-up thanks to the
        interface with fortran. Memory consuming though...

        Parameters
        ----------
        waferts : ndarray
            Array of timestreams. Size (ndetectors, ntimesamples).
        output_maps : OutputSkyMap instance
            Instance of OutputSkyMap which contains the sky maps. The
            coaddition of data is done on-the-fly directly.
        gdeprojection : bool, optional
            If True, perform a deprojection of a constant contribution in the
            polarisation timestream: d^{-} = G + Qcos(2*theta) + U*sin(2theta).
            Work only for pair difference.
        frequency_channel : int, optional
            If you are processing dichroic pixels, you need to specify the
            index of the frequency channel (1 or 2). Default is 1.

        Examples
        ----------
        HEALPIX: Test the routines MAP -> TOD -> MAP (mapping all detectors
        at once).
        >>> inst, scan, sky_in = load_fake_instrument()
        >>> tod = TimeOrderedDataPairDiff(inst, scan, sky_in,
        ...     CESnumber=0, projection='healpix', mapping_perpair=False)
        >>> d = np.array([tod.map2tod(det) for det in range(2 * tod.npair)])
        >>> m = OutputSkyMap(projection=tod.projection,
        ...     nside=tod.nside_out, obspix=tod.obspix)
        >>> tod.tod2map(d, m)

        HEALPIX: Test the routines MAP -> TOD -> MAP with mapping_perpair to
        save memory (only 2 TOD at a time are loaded)
        >>> inst, scan, sky_in = load_fake_instrument()
        >>> tod = TimeOrderedDataPairDiff(inst, scan, sky_in,
        ...     CESnumber=0, projection='healpix', mapping_perpair=True)
        >>> m = OutputSkyMap(projection=tod.projection,
        ...     nside=tod.nside_out, obspix=tod.obspix)
        >>> for pair in tod.pair_list:
        ...   d = np.array([tod.map2tod(det) for det in pair])
        ...   tod.tod2map(d, m)

        Check intensity map
        >>> sky_out = np.zeros(12 * tod.nside_out**2)
        >>> sky_out[tod.obspix] = m.get_I()
        >>> mask = sky_out != 0.0
        >>> assert np.allclose(sky_out[mask], sky_in.I[mask])

        Check polarisation maps
        >>> sky_out = np.zeros((2, 12 * tod.nside_out**2))
        >>> sky_out[0][tod.obspix], sky_out[1][tod.obspix] = m.get_QU()
        >>> mask = (sky_out[0] != 0.0) * (sky_out[1] != 0.0)
        >>> assert np.allclose(sky_out[0][mask], sky_in.Q[mask])
        >>> assert np.allclose(sky_out[1][mask], sky_in.U[mask])

        FLAT: Test the routines MAP -> TOD -> MAP.
        >>> inst, scan, sky_in = load_fake_instrument()
        >>> tod = TimeOrderedDataPairDiff(inst, scan, sky_in,
        ...     CESnumber=0, projection='flat')
        >>> d = np.array([tod.map2tod(det) for det in range(2 * tod.npair)])
        >>> m = OutputSkyMap(projection=tod.projection,
        ...     npixsky=tod.npixsky, pixel_size=tod.pixel_size)
        >>> tod.tod2map(d, m)

        HEALPIX + Dichroic
        >>> inst, scan, sky_in = load_fake_instrument(fwhm_in2=1.8)
        >>> tod = TimeOrderedDataPairDiff(inst, scan, sky_in, mode='dichroic',
        ...     CESnumber=0, projection='healpix', mapping_perpair=True)
        >>> m1 = OutputSkyMap(projection=tod.projection,
        ...     nside=tod.nside_out, obspix=tod.obspix)
        >>> m2 = OutputSkyMap(projection=tod.projection,
        ...     nside=tod.nside_out, obspix=tod.obspix)
        >>> for pair in tod.pair_list:
        ...   d = np.array([tod.map2tod(det) for det in pair])
        ...   d1 = d[:, 0] # first frequency channel
        ...   d2 = d[:, 1] # second frequency channel
        ...   tod.tod2map(d1, m1)
        ...   tod.tod2map(d2, m2)

        HEALPIX + deprojection
        >>> inst, scan, sky_in = load_fake_instrument()
        >>> tod = TimeOrderedDataPairDiff(inst, scan, sky_in,
        ...     CESnumber=0, projection='healpix', mapping_perpair=True)
        >>> m = OutputSkyMapIGQU(projection=tod.projection,
        ...     nside=tod.nside_out, obspix=tod.obspix)
        >>> for pair in tod.pair_list:
        ...   d = np.array([tod.map2tod(det) for det in pair])
        ...   tod.tod2map(d, m, gdeprojection=True)

        """
        if frequency_channel == 1:
            pol_angs = self.pol_angs
        elif frequency_channel == 2:
            pol_angs = self.pol_angs2

        nbolofp = waferts.shape[0]
        npixfp = nbolofp / 2
        nt = int(waferts.shape[-1])

        # Check sizes
        msg = """
            Most likely you set mapping_perpair wrongly when
            initialising your TOD.
            The way you loop over focal plane pixel to do the mapmaking
            depends on the value of mapping_perpair parameter.
            If mapping_perpair=False, one first load all the pixel and
            then you need to perform the mapmaking with all the pixels
            at once. If mapping_perpair=True, one should load
            pair-by-pair and the mapmaking is done pair-by-pair.
            See so_MC_crosstalk.py vs so_MC_gain_drift.py to see both
            approaches (s4cmb-resources/Part2), and example in doctest above.
        """
        assert npixfp == self.point_matrix.shape[0], msg
        assert nt == self.point_matrix.shape[1], msg

        assert npixfp == pol_angs.shape[0], msg
        assert nt == pol_angs.shape[1], msg

        assert npixfp == self.diff_weight.shape[0], msg
        assert npixfp == self.sum_weight.shape[0], msg
        point_matrix = self.point_matrix.flatten()
        pol_angs = pol_angs.flatten()
        waferts = waferts.flatten()
        diff_weight = self.diff_weight.flatten()
        sum_weight = self.sum_weight.flatten()
        wafermask_pixel = self.wafermask_pixel.flatten()

        if hasattr(self, "dm") and (gdeprojection is False):
            tod_f.tod2map_hwp_f(
                output_maps.d0,
                output_maps.d4r,
                output_maps.d4i,
                output_maps.w0,
                output_maps.w4,
                output_maps.nhit,
                point_matrix,
                pol_angs,
                waferts,
                diff_weight,
                sum_weight,
                nt,
                wafermask_pixel,
                npixfp,
                self.npixsky,
            )
        elif (hasattr(output_maps, "dm") is False) and gdeprojection:
            tod_f.tod2map_pair_gdeprojection_f(
                output_maps.d,
                output_maps.w,
                output_maps.dm,
                output_maps.dc,
                output_maps.ds,
                output_maps.wm,
                output_maps.cc,
                output_maps.cs,
                output_maps.ss,
                output_maps.c,
                output_maps.s,
                output_maps.nhit,
                point_matrix,
                pol_angs,
                waferts,
                diff_weight,
                sum_weight,
                nt,
                wafermask_pixel,
                npixfp,
                self.npixsky,
            )
        else:
            tod_f.tod2map_pair_f(
                output_maps.d,
                output_maps.w,
                output_maps.dc,
                output_maps.ds,
                output_maps.cc,
                output_maps.cs,
                output_maps.ss,
                output_maps.nhit,
                point_matrix,
                pol_angs,
                waferts,
                diff_weight,
                sum_weight,
                nt,
                wafermask_pixel,
                npixfp,
                self.npixsky,
            )
        # Garbage collector guard
        wafermask_pixel class OutputSkyMap:
    """ Class to handle sky maps generated by tod2map """

    def __init__(
        self,
        projection,
        obspix=None,
        npixsky=None,
        nside=None,
        pixel_size=None,
        demodulation=False,
    ):
        """
        Initialise all maps: weights, projected TOD, and Stokes parameters.

        Parameters
        ----------
        projection : string
            Type of projection among [healpix, flat].
        obspix : 1d array, optional
            List of indices of observed pixels if projection=healpix. No effect
            if projection=flat.
        npixsky : int, optional
            The number of observed sky pixels in projection=flat.
            npixsky is by default len(obspix) if projection=healpix.
        nside : int, optional
            The resolution for the output map if projection=healpix. No effect
            if projection=flat.
        pixel_size : float, optional
            The size of pixels in arcmin if projection=flat. No effect
            if projection=healpix.
        demodulation : bool, False
            If True, assume you are working with demodulated timestreams.
            Default is False.
        """
        self.nside = nside
        self.projection = projection
        self.obspix = obspix
        self.npixsky = npixsky
        self.nside = nside
        self.pixel_size = pixel_size
        self.demodulation = demodulation

        if self.projection == "healpix":
            assert self.obspix is not None, ValueError(
                """You need to provide the list (obspix)
                of observed pixel if projection=healpix!"""
            )
            assert self.nside is not None, ValueError(
                """You need to provide the resolution (nside)
                of the map if projection=healpix!"""
            )
            self.npixsky = len(self.obspix)
            self.pixel_size = hp.nside2resol(nside, arcmin=True)

        elif self.projection == "flat":
            assert self.npixsky is not None, ValueError(
                """You need to provide the number of
                observed pixels (npixsky) if projection=flat."""
            )
            assert self.pixel_size is not None, ValueError(
                """You need to provide the size of
                pixels (pixel_size) in arcmin if projection=flat.
                """
            )

        self.initialise_sky_maps()

        if not self.demodulation:
            self.to_coadd = "d dc ds w cc cs ss nhit"
        else:
            self.to_coadd = "d0 d4r d4i w0 w4 nhit"

    def initialise_sky_maps(self):
        """
        Create empty sky maps. This includes:
        * d : projected weighted sum of timestreams
        * dc : projected noise weighted difference of timestreams
            multiplied by cosine
        * ds : projected noise weighted difference of timestreams
            multiplied by sine
        * nhit : projected hit counts.
        * w : projected (inverse) noise weights and hits.
        * cc : projected noise weighted cosine**2
        * ss : projected noise weighted sine**2
        * cs : projected noise weighted cosine * sine.

        Work with pair differenced data only.
        See initialise_sky_maps_demod for demodulated data.

        """
        if self.demodulation:
            self.initialise_sky_maps_demod()
        else:
            # To accumulate A^T N^-1 d
            self.d = np.zeros(self.npixsky)
            self.dc = np.zeros(self.npixsky)
            self.ds = np.zeros(self.npixsky)

            # To accumulate A^T N^-1 A
            self.w = np.zeros(self.npixsky)
            self.cc = np.zeros(self.npixsky)
            self.cs = np.zeros(self.npixsky)
            self.ss = np.zeros(self.npixsky)

            self.nhit = np.zeros(self.npixsky, dtype=np.int32)

    def set_idet(self):
        """
        Return the inverse determinant of the weight matrix for each sky pixel.
        Bad pixels have their inverse determinant set to zero.

        Returns
        ----------
        idet : 1D arrray
            Array containing the inverse determinant of the weight matrix for
            each sky pixel.
        """
        testcc = self.cc * self.ss - self.cs * self.cs
        idet = np.zeros(testcc.shape)
        inonzero = testcc != 0.0
        idet[inonzero] = 1.0 / testcc[inonzero]

        thresh = np.finfo(np.float32).eps
        try:
            izero = np.abs(testcc) < thresh
        except FloatingPointError:
            izero = inan = np.isnan(testcc)
            izero[~inan] = np.abs(testcc[~inan]) < thresh

        idet[izero] = 0.0
        self.idet = idet

    def get_I(self):
        """
        Solve for the intensity map I from projected sum timestream map d
        and weights w: w * I = d. Work with pair differenced data only.
        See get_I_demod for demodulated data.

        Returns
        ----------
        I : 1d array
            Intensity map. Note that only the observed pixels defined in
            obspix are returned (and not the full sky map).

        Examples
        ----------
        Pair difference data
        >>> m1 = OutputSkyMap(projection='healpix',
        ...     nside=16, obspix=np.array([0, 1, 2, 3]), demodulation=False)
        >>> I = m1.get_I()
        """
        if self.demodulation:
            return self.get_I_demod()

        hit = self.w > 0
        I = np.zeros_like(self.d)
        I[hit] = self.d[hit] / self.w[hit]
        return I

    def get_QU(self, force=True):
        """
        Solve for the polarisation maps from projected difference timestream
        maps and weights:

        [cc cs]   [Q]   [dc]
        [cs ss] * [U] = [ds]

        Work with pair differenced data only. See get_QU_demod for
        demodulated data.

        Returns
        ----------
        Q : 1d array
            Stokes Q map. Note that only the observed pixels defined in
            obspix are returned (and not the full sky map).
        U : 1d array
            Stokes U map. Note that only the observed pixels defined in
            obspix are returned (and not the full sky map).

        Examples
        ----------
        Pair difference data
        >>> m1 = OutputSkyMap(projection='healpix',
        ...     nside=16, obspix=np.array([0, 1, 2, 3]), demodulation=False)
        >>> Q, U = m1.get_QU()
        """
        if self.demodulation:
            return self.get_QU_demod()

        if (not hasattr(self, "idet")) or force:
            self.set_idet()

        Q = self.idet * (self.ss * self.dc - self.cs * self.ds)
        U = self.idet * (-self.cs * self.dc + self.cc * self.ds)

        return Q, U

    def get_IQU(self):
        """
        Solve for the temperature and polarisation maps from
        projected sum and difference timestream maps and weights:

        [w 0  0 ]   [I]   [d ]
        [0 cc cs]   [Q]   [dc]
        [0 cs ss] * [U] = [ds]

        or demodulated data:

        [w0  0  0]   [I]   [d0 ]
        [ 0 w4  0]   [Q]   [d4r]
        [ 0  0 w4] * [U] = [d4i]

        Examples
        ----------
        Pair difference data
        >>> m1 = OutputSkyMap(projection='healpix',
        ...     nside=16, obspix=np.array([0, 1, 2, 3]), demodulation=False)
        >>> I, Q, U = m1.get_IQU()

        Demodulation
        >>> m1 = OutputSkyMap(projection='healpix',
        ...     nside=16, obspix=np.array([0, 1, 2, 3]), demodulation=True)
        >>> I, Q, U = m1.get_IQU()

        Returns
        ----------
        I : 1d array
            Intensity map. Note that only the observed pixels defined in
            obspix are returned (and not the full sky map).
        Q : 1d array
            Stokes Q map. Note that only the observed pixels defined in
            obspix are returned (and not the full sky map).
        U : 1d array
            Stokes U map. Note that only the observed pixels defined in
            obspix are returned (and not the full sky map).
        """
        I = self.get_I()
        Q, U = self.get_QU()
        return I, Q, U

    def coadd(self, other, to_coadd=None):
        """
        Add other\'s vectors into our vectors.

        Note:
        You do not need this routine most of the case as
        tod2map operates a coaddition internally if you pass the
        same OutputSkyMap instance.

        Parameters
        ----------
        other : OutputSkyMap instance
            Instance of OutputSkyMap to be coadded with this one.
        to_coadd : string, optional
            String with names of vectors to coadd separated by a space.
            Names must be attributes of other and self.

        Examples
        ---------
        Coadd two maps together.
        >>> m1 = OutputSkyMap(projection='healpix',
        ...     nside=16, obspix=np.array([0, 1, 2, 3]))
        >>> m1.nhit = np.ones(4)
        >>> m2 = OutputSkyMap(projection='healpix',
        ...     nside=16, obspix=np.array([0, 1, 2, 3]))
        >>> m2.nhit = np.ones(4)
        >>> m1.coadd(m2)
        >>> print(m1.nhit)
        [ 2.  2.  2.  2.]

        Same idea for demodulated data.
        >>> m1 = OutputSkyMap(projection='healpix',
        ...     nside=16, obspix=np.array([0, 1, 2, 3]), demodulation=True)
        >>> m1.nhit = np.ones(4)
        >>> m2 = OutputSkyMap(projection='healpix',
        ...     nside=16, obspix=np.array([0, 1, 2, 3]), demodulation=True)
        >>> m2.nhit = np.ones(4)
        >>> m1.coadd(m2)
        >>> print(m1.nhit)
        [ 2.  2.  2.  2.]

        Same idea in flat sky.
        >>> m1 = OutputSkyMap(projection='flat',
        ...     npixsky=4, pixel_size=2.)
        >>> m1.nhit = np.ones(4)
        >>> m2 = OutputSkyMap(projection='flat',
        ...     npixsky=4, pixel_size=2.)
        >>> m2.nhit = np.ones(4)
        >>> m1.coadd(m2)
        >>> print(m1.nhit)
        [ 2.  2.  2.  2.]
        """
        assert np.all(self.obspix == other.obspix), ValueError(
            "To add maps together, they must have the same obspix!"
        )

        if to_coadd is None:
            to_coadd = self.to_coadd

        to_coadd_split = to_coadd.split(" ")
        for k in to_coadd_split:
            a = getattr(self, k)
            b = getattr(other, k)
            a += b

    def coadd_MPI(self, other, MPI, to_coadd=None):
        """
        Crappy way of coadding vectors through different processors.

        Parameters
        ----------
        other : OutputSkyMap instance
            Instance of OutputSkyMap to be coadded with this one.
        MPI : module
            Module for communication. It has been tested through mpi4py only
            for the moment.
        to_coadd : string, optional
            String with names of vectors to coadd separated by a space.
            Names must be attributes of other and self.

        Examples
        ---------
        Coadd maps from different processors together.
        >>> from mpi4py import MPI
        >>> m = OutputSkyMap(projection='healpix',
        ...     nside=16, obspix=np.array([0, 1, 2, 3]))
        >>> # do whatever you want with the maps
        >>> m.coadd_MPI(m, MPI)
        """
        if to_coadd is None:
            to_coadd = self.to_coadd

        to_coadd_split = to_coadd.split(" ")
        for k in to_coadd_split:
            setattr(
                self,
                k,
                MPI.COMM_WORLD.allreduce(getattr(other, k), op=MPI.SUM),
            )

    def pickle_me(
        self, fn, shrink_maps=True, crop_maps=False, epsilon=0.0, verbose=False
    ):
        """
        Save data into pickle file.
        Work with pair differenced data only. See pickle_me_demod for
        demodulated data.

        Parameters
        ----------
        fn: string
            The name of the file where data will be stored.
        shrink_maps : bool, optional
            If true, remove unecessary zeros before storing the maps.
        crop_me : False or int
            If not False, crop the output maps as map[:crop_me][:crop_me].
        epsilon : float, optional
            Threshold for selecting the pixels in polarisation.
            0 <= epsilon < 1/4. The higher the more selective.

        """
        if self.demodulation:
            self.pickle_me_demod(
                fn,
                shrink_maps=shrink_maps,
                crop_maps=crop_maps,
                epsilon=epsilon,
                verbose=verbose,
            )
        else:
            try:
                I, Q, U = self.get_IQU()
                wP = qu_weight_mineig(
                    self.cc, self.cs, self.ss, epsilon=epsilon, verbose=verbose
                )

                data = {
                    "I": I,
                    "Q": Q,
                    "U": U,
                    "wI": self.w,
                    "wP": wP,
                    "nhit": self.nhit,
                    "projection": self.projection,
                    "nside": self.nside,
                    "pixel_size": self.pixel_size,
                    "obspix": self.obspix,
                }
            except Exception as e:
                print("Exception error: ", e)
                I, G, Q, U = self.get_IQU()  # if using IGQU class
                wP = qu_weight_mineig(
                    self.cc, self.cs, self.ss, epsilon=epsilon, verbose=verbose
                )

                data = {
                    "I": I,
                    "G": G,
                    "Q": Q,
                    "U": U,
                    "wI": self.w,
                    "wP": wP,
                    "nhit": self.nhit,
                    "projection": self.projection,
                    "nside": self.nside,
                    "pixel_size": self.pixel_size,
                    "obspix": self.obspix,
                }

            if shrink_maps and self.projection == "flat":
                data = shrink_me(data, based_on="wP")
            elif crop_maps is not False and self.projection == "flat":
                data = crop_me(data, based_on="wP", npix_per_row=crop_maps)

            with open(fn, "wb") as f:
                pickle.dump(data, f, protocol=2)

    def initialise_sky_maps_demod(self):
        """
        Create empty sky maps. This includes:
        * d : projected weighted sum of timestreams
        * dc : projected noise weighted difference of timestreams
            multiplied by cosine
        * ds : projected noise weighted difference of timestreams
            multiplied by sine
        * nhit : projected hit counts.
        * w : projected (inverse) noise weights and hits.
        * cc : projected noise weighted cosine**2
        * ss : projected noise weighted sine**2
        * cs : projected noise weighted cosine * sine.

        Work with demodulated data only. See initialise_sky_maps for
        pair differenced.

        """
        # To accumulate A^T N^-1 d
        self.d0 = np.zeros(self.npixsky)
        self.d4r = np.zeros(self.npixsky)
        self.d4i = np.zeros(self.npixsky)

        # To accumulate A^T N^-1 A
        self.w0 = np.zeros(self.npixsky)
        self.w4 = np.zeros(self.npixsky)

        self.nhit = np.zeros(self.npixsky, dtype=np.int32)

    def get_I_demod(self):
        """
        Solve for the intensity map I from projected sum timestream map d
        and weights w0: w0 * I = d0. Work with demodulated data only.

        Returns
        ----------
        I : 1d array
            Intensity map. Note that only the observed pixels defined in
            obspix are returned (and not the full sky map).

        Examples
        ----------
        Demodulated data
        >>> m1 = OutputSkyMap(projection='healpix',
        ...     nside=16, obspix=np.array([0, 1, 2, 3]), demodulation=True)
        >>> I = m1.get_I()
        """
        hit = self.w0 > 0
        I = np.zeros_like(self.d0)
        I[hit] = self.d0[hit] / self.w0[hit]
        return I

    def get_QU_demod(self):
        """
        Solve for the polarisation maps from projected difference timestream
        maps and weights:

        [w4  0]   [Q]   [d4r]
        [0  w4] * [U] = [d4i]

        Work with demodulated data only.

        Returns
        ----------
        Q : 1d array
            Stokes Q map. Note that only the observed pixels defined in
            obspix are returned (and not the full sky map).
        U : 1d array
            Stokes U map. Note that only the observed pixels defined in
            obspix are returned (and not the full sky map).

        Examples
        ----------
        Demodulated data
        >>> m1 = OutputSkyMap(projection='healpix',
        ...     nside=16, obspix=np.array([0, 1, 2, 3]), demodulation=True)
        >>> Q, U = m1.get_QU()
        """
        hit = self.w4 > 0

        Q = np.zeros_like(self.d4r)
        U = np.zeros_like(self.d4i)

        Q[hit] = self.d4r[hit] / self.w4[hit]
        U[hit] = self.d4i[hit] / self.w4[hit]

        return Q, U

    def pickle_me_demod(
        self, fn, shrink_maps=True, crop_maps=False, epsilon=0.0, verbose=False
    ):
        """
        Save data into pickle file.
        Work with demodulated data only.

        Parameters
        ----------
        fn: string
            The name of the file where data will be stored.
        shrink_maps : bool, optional
            If true, remove unecessary zeros before storing the maps.
        crop_me : False or int
            If not False, crop the output maps as map[:crop_me][:crop_me].
        epsilon : float, optional
            Threshold for selecting the pixels in polarisation.
            0 <= epsilon < 1/4. The higher the more selective.

        """
        I, Q, U = self.get_IQU()

        data = {
            "I": I,
            "Q": Q,
            "U": U,
            "wI": self.w0,
            "wP": self.w4,
            "nhit": self.nhit,
            "projection": self.projection,
            "nside": self.nside,
            "pixel_size": self.pixel_size,
            "obspix": self.obspix,
        }

        if shrink_maps and self.projection == "flat":
            data = shrink_me(data, based_on="wP")
        elif crop_maps is not False and self.projection == "flat":
            data = crop_me(data, based_on="wP", npix_per_row=crop_maps)

        with open(fn, "wb") as f:
            pickle.dump(data, f, protocol=2) def inject_crosstalk_inside_SQUID(
    bolo_data,
    squid_ids,
    bolo_ids,
    mu=-3.0,
    sigma=1.0,
    radius=1,
    beta=2,
    seed=5438765,
    new_array=None,
    language="python",
):
    """
    Introduce leakage between neighboring bolometers within a SQUID.
    You have to provide the list of bolometers, to which SQUID they
    belong, and the index of bolometers within each SQUID.

    Parameters
    ----------
    bolo_data : list of arrays
        Contains the time-ordered data of bolometers.
    squid_ids : list of strings
        Contains the SQUID id for each bolometer.
        Should have the same shape[0] as bolo_data.
    bolo_ids : list of strings
        Contains the positions of bolometers within their SQUID.
        Should have the same shape[0] as bolo_data.
    mu : float, optional
        Mean of the Gaussian used to generate the level of leakage,
        in percent. E.g. mu=1.0 means the leakage coefficients will be
        centred around 1%.
    sigma : float, optional
        Width of the Gaussian used to generate the level of leakage,
        in percent. E.g. sigma=1.0 means the leakage coefficients will be
        of order X \pm 1%.
    radius : int, optional
        Controls the number of bolometers talking within a SQUID.
        radius=n means bolometers [N-n, N-n+1, ..., N, ..., N+n-1, N+n]
        will be crosstalking. `radius` must be at most equal to the number of
        bolometers within a SQUID.
    beta : int, optional
        Exponent controling the attenuation. The idea is that crosstalk depends
        on the separation length. Bolometer $b_i$ crosstalking with bolometers
        $b_{i \pm d}$ ($d \leq FF$) with leakage coefficients
        $\alpha_{i \pm d}$ will see his timestream modified as:
            $b_i = b_i + \dfrac{\alpha_{i \pm d}}{d^\beta}  b_{i \pm d}$.
        Assuming some fixed frequency spacing, hardware consideration
        suggests $\beta = 2$.
    seed : int, optional
        Control the random seed used to generate leakage coefficients.
    new_array : None or bolo_data-like array, optional
        If not None, return a new array of timestreams with the modifications.
        Modify bolo_data directly otherwise. Default is None.
    language : string, optional
        Language to perform computations to be chosen in ['python', 'fortran'].
        Default is python. [fortran is very slow...]

    Example
    ----------
    >>> from s4cmb.tod import load_fake_instrument
    >>> from s4cmb.tod import TimeOrderedDataPairDiff
    >>> inst, scan, sky_in = load_fake_instrument()
    >>> tod = TimeOrderedDataPairDiff(inst, scan, sky_in, CESnumber=0)
    >>> d = np.array([tod.map2tod(det) for det in range(2 * tod.npair)])
    >>> print(round(d[0][0], 3))
    40.95

    Inject crosstalk between neighbour bolometers (radius=1)
    >>> squid_ids = inst.focal_plane.get_indices('Sq')
    >>> bolo_ids = inst.focal_plane.bolo_index_in_squid
    >>> inject_crosstalk_inside_SQUID(d, squid_ids, bolo_ids, radius=1)
    >>> print(round(d[0][0], 3))
    39.561

    One can also keep the original timestreams, and return a new array
    containing modified timestreams.
    >>> d = np.array([tod.map2tod(det) for det in range(2 * tod.npair)])
    >>> d_new = np.zeros_like(d)
    >>> inject_crosstalk_inside_SQUID(d, squid_ids, bolo_ids, radius=1,
    ...     new_array=d_new, language='python')
    >>> print(round(d[0][0], 3), round(d_new[0][0], 3))
    ... #doctest: +NORMALIZE_WHITESPACE
    40.95 39.561

    For large number of bolometers per SQUID, you would prefer fortran
    to python to perform the loops. Choose python otherwise.
    """
    # Make mu and sigma unitless (user provides percentage)
    mu = mu / 100.0
    sigma = sigma / 100.0

    combs = {}
    for bolo in range(len(bolo_data)):
        sq = squid_ids[bolo]
        if sq not in combs:
            combs[sq] = []
        combs[sq].append((bolo_ids[bolo], bolo))

    tsout = 0.0 + bolo_data
    # How much to leak from one bolometer to its neighboring channels
    state = np.random.RandomState(seed)
    cross_amp = state.normal(mu, sigma, len(bolo_data))

    if language == "python":
        for sq in combs:
            for ch, i in combs[sq]:
                for ch2, i2 in combs[sq]:
                    separation_length = abs(ch - ch2)
                    if separation_length > 0 and separation_length <= radius:
                        tsout[i] += (
                            cross_amp[i2] / separation_length ** beta * tsout[i2]
                        )

    elif language == "fortran":
        # F2PY convention
        tsout = np.array(tsout, order="F")
        for sq in combs:
            local_indices = np.array(combs[sq]).flatten()[::2]
            global_indices = np.array(combs[sq]).flatten()[1::2]
            systematics_f.inject_crosstalk_inside_squid_f(
                tsout,
                local_indices,
                global_indices,
                radius,
                cross_amp,
                beta,
                len(local_indices),
                len(bolo_data),
                len(bolo_data[0]),
            )

    if new_array is not None:
        new_array[:] = tsout
    else:
        bolo_data[:] = tsout def import_string_as_module(fn_full):
    """
    Import module from its name given as a string.

    Parameters
    ----------
    fn_full: string
        Python filename containing parameters that you
        want to import.

    Returns
    ----------
    params: module
        Module containing your parameters.

    Examples
    ----------
    >>> fn_full = 'examples/inifiles/simple_parameters.py'
    >>> params = import_string_as_module(fn_full)
    >>> 'do_pol' in dir(params)
    True

    """
    # Import parameters from the user parameter file
    fn_short = os.path.basename(fn_full).split(".py")[0]
    sys.path.insert(
        0, os.path.realpath(os.path.join(os.getcwd(), os.path.dirname(fn_full)))
    )
    params = importlib.import_module(fn_short)
    return params def create_batch(batch_file, name_out, params_s4cmb, params_xpure):
    """
    Create submission file for the software xpure

    Parameters
    ----------
    batch_file : string
        Filename.
    params_s4cmb : NormaliseS4cmbParser instance
        Object with s4cmb parameter values.
    params_xpure : NormaliseXpureParser instance
        Object with xpure parameter values.
    """
    host = os.environ["NERSC_HOST"]
    if host == "edison":
        params_xpure.nproc_per_node = 24
    elif host == "cori":
        params_xpure.nproc_per_node = 32

    with open(batch_file, "w") as f:
        print("#!/bin/bash -l", file=f)
        print("#SBATCH -p {}".format(params_xpure.queue), file=f)
        print("#SBATCH -N {}".format(params_xpure.node), file=f)
        print("#SBATCH -t {}".format(params_xpure.time), file=f)
        print("#SBATCH -J {}".format(name_out), file=f)
        if host == "cori":
            print("#SBATCH -C haswell", file=f)

        print(" ", file=f)
        print("#### TO BE CHANGED BY USER ####", file=f)
        print("# Name for files and folders", file=f)
        print("name={}".format(name_out), file=f)

        print(" ", file=f)
        print("# Radius for apodization", file=f)
        print("radius={}".format(params_xpure.radius_apodization), file=f)

        print(" ", file=f)
        print("# Parameters for XPURE", file=f)
        print("LMAX_USER={}".format(params_xpure.lmax_user), file=f)
        print("NSIDE_USER={}".format(params_s4cmb.nside_out), file=f)

        print(" ", file=f)
        print("# Number of simulation (put 1 if not simulation)", file=f)
        print("NSIMU_USER=1", file=f)

        print(" ", file=f)
        print(
            "# Number of masks (put one if several maps but one common masks)",
            file=f,
        )
        print("NMASK_USER=1", file=f)

        print(" ", file=f)
        print("# Number of maps (if > 1, will do the cross correlations)", file=f)
        print("NMAPS_USER=1", file=f)

        print(" ", file=f)
        print("# Beam and bins", file=f)
        print("BEAMFILE={}".format(params_xpure.beam_file), file=f)
        print("BINTAB={}".format(params_xpure.bin_file), file=f)

        print(" ", file=f)
        print("# XPURE running mode (0=xpol, 1=xpure, 2=hybrid)", file=f)
        print("MODE_XPURE={}".format(params_xpure.xpure_mode), file=f)

        print(" ", file=f)
        txt = "# FULL (0=myapodizemask, create_mll and XPURE) or FAST "\
            + "(1=only XPURE) or SEMI-FAST (=2 create_mll and XPURE)"
        print(
            txt,
            file=f,
        )
        print("FAST={}".format(params_xpure.fast), file=f)

        print(" ", file=f)
        print("#### END TO BE CHANGED BY USER ####", file=f)

        print(" ", file=f)
        print("#ENVIRONEMENT", file=f)
        print("SCRATCH2=/global/cscratch1/sd/peloton", file=f)
        if host == "edison":
            print("BINDIR=${HOME}/mapmaker/xpure/xpure/trunk/build/edison", file=f)
        elif host == "cori":
            print("BINDIR=${HOME}/mapmaker/xpure/xpure/trunk/build/cori", file=f)
        print("THEORYDIR=${SCRATCH2}/s4cmb/additional_files", file=f)
        print("OUTPUTMLL=${SCRATCH2}/s4cmb/xpure/cls/${name}", file=f)
        print("OUTPUTCL=${SCRATCH2}/s4cmb/xpure/cls/${name}", file=f)
        print("mkdir -p ${OUTPUTMLL}", file=f)
        print("mkdir -p ${OUTPUTCL}", file=f)
        print("MASKDIR=${SCRATCH2}/s4cmb/xpure/masks/${name}", file=f)
        print("MAPDIR=${SCRATCH2}/s4cmb/xpure/maps/${name}", file=f)

        print(" ", file=f)
        print("##########################", file=f)
        print("# MASK (COMMON MASK)", file=f)
        print("##########################", file=f)
        print("cd $MASKDIR", file=f)

        print(" ", file=f)
        print('binary_name=$(basename $(find -name "Iw*${name}_norm.fits"))', file=f)
        print('weight_name=$(basename $(find -name "Iw*${name}.fits"))', file=f)
        print("BINARY_MASK_I1=${MASKDIR}/${binary_name}", file=f)
        print("WEIGHT_I1=${MASKDIR}/${weight_name}", file=f)
        print("APODIZED_MASK_I1=${MASKDIR}/ap30_${weight_name}", file=f)

        print(" ", file=f)
        print('binary_name=$(basename $(find -name "Pw*${name}_norm.fits"))', file=f)
        print('weight_name=$(basename $(find -name "Pw*${name}.fits"))', file=f)
        print("BINARY_MASK_P1=${MASKDIR}/${binary_name}", file=f)
        print("WEIGHT_P1=${MASKDIR}/${weight_name}", file=f)
        print("APODIZED_MASK_P1=${MASKDIR}/ap30_${weight_name}", file=f)

        print(" ", file=f)
        print("## Intensity", file=f)
        print("output_spin0_I1=${MASKDIR}/spin0_I_${name}.fits", file=f)
        print("output_spin1_I1=${MASKDIR}/spin1_I_${name}.fits", file=f)
        print("output_spin2_I1=${MASKDIR}/spin2_I_${name}.fits", file=f)

        print(" ", file=f)
        print("## Polarization", file=f)
        print("output_spin0_P1=${MASKDIR}/spin0_P_${name}.fits", file=f)
        print("output_spin1_P1=${MASKDIR}/spin1_P_${name}.fits", file=f)
        print("output_spin2_P1=${MASKDIR}/spin2_P_${name}.fits", file=f)

        print(" ", file=f)
        print("##########################", file=f)
        print("# MAPS", file=f)
        print("##########################", file=f)
        print("cd ${MAPDIR}", file=f)

        print(" ", file=f)
        print("for (( i=1; i<=${NMAPS_USER}; i++)); do", file=f)
        print('	MAPS[${i}]=$(basename $(find -name "IQU*${name}.fits"))', file=f)
        print("done", file=f)

        print(" ", file=f)
        print("# CHECK", file=f)
        print("for (( i=1; i<=${NMAPS_USER}; i++)); do", file=f)
        print("	echo ${MAPS[${i}]}", file=f)
        print("done", file=f)

        print(" ", file=f)
        print("cd $PBS_O_WORKDIR", file=f)

        print(" ", file=f)
        print(
            "################################",
            file=f,
        )
        print("# From binary mask to binary apodized mask", file=f)
        print(
            "################################",
            file=f,
        )
        print('if [ "${FAST}" -eq "0" ]', file=f)
        print("	then", file=f)
        txt = "	time srun -N 1 -n {} ${{BINDIR}}/myapodizemask "\
            + "${{BINARY_MASK_I1}} ${{APODIZED_MASK_I1}} -minpix 1 -inside 1 "\
            + "-radius ${{radius}} & time srun -N 1 -n {} "\
            + "${{BINDIR}}/myapodizemask ${{BINARY_MASK_P1}} "\
            + "${{APODIZED_MASK_P1}} -minpix 1 -inside 1 -radius ${{radius}}"
        print(
            txt.format(
                params_xpure.nproc_apo, params_xpure.nproc_apo
            ),
            file=f,
        )
        print("	wait", file=f)
        print("else", file=f)
        print('	echo "Go fast - skip myapodizemask"', file=f)
        print("fi", file=f)

        print(" ", file=f)
        print(
            "################################",
            file=f,
        )
        print("# Computation of spin window function", file=f)
        print(
            "################################",
            file=f,
        )

        print(" ", file=f)
        print("cat > param_all_I11${name}.par << EOF", file=f)
        print("#healpix parameters", file=f)
        print("##########", file=f)
        print("nside = ${NSIDE_USER}", file=f)
        print("lmax = ${LMAX_USER} #should not exceed 2*Nside", file=f)

        print(" ", file=f)
        print("#input file parameters", file=f)
        print("############", file=f)
        print("maskBinary = ${BINARY_MASK_I1}", file=f)
        print("window_spin0 = ${APODIZED_MASK_I1}", file=f)
        print("inverseNoise = ${WEIGHT_I1}", file=f)

        print(" ", file=f)
        print("#output file parameters", file=f)
        print("############", file=f)
        print("output_spin0 = ${output_spin0_I1}", file=f)
        print("output_spin1 = ${output_spin1_I1}", file=f)
        print("output_spin2 = ${output_spin2_I1}", file=f)
        print("EOF", file=f)

        print(" ", file=f)
        print("cat > param_all_P11${name}.par << EOF", file=f)
        print("#healpix parameters", file=f)
        print("##########", file=f)
        print("nside = ${NSIDE_USER}", file=f)
        print("lmax = ${LMAX_USER} #should not exceed 2*Nside", file=f)

        print(" ", file=f)
        print("#input file parameters", file=f)
        print("############", file=f)
        print("maskBinary = ${BINARY_MASK_P1}", file=f)
        print("window_spin0 = ${APODIZED_MASK_P1}", file=f)
        print("inverseNoise = ${WEIGHT_P1}", file=f)

        print(" ", file=f)
        print("#output file parameters", file=f)
        print("############", file=f)
        print("output_spin0 = ${output_spin0_P1}", file=f)
        print("output_spin1 = ${output_spin1_P1}", file=f)
        print("output_spin2 = ${output_spin2_P1}", file=f)
        print("EOF", file=f)

        print(" ", file=f)
        print('if [ "${FAST}" -eq "0" ]', file=f)
        print("	then", file=f)
        txt = "	time srun -N {} -n {} ${{BINDIR}}/scalar2spin "\
            + "param_all_I11${{name}}.par >& output_scalar2spinI11${{name}} "\
            + "& time srun -N {} -n {} "\
            + "${{BINDIR}}/scalar2spin param_all_P11${{name}}.par "\
            + ">& output_scalar2spinP11${{name}}"
        print(
            txt.format(
                int(params_xpure.nproc_scalar_to_spin // params_xpure.nproc_per_node),
                params_xpure.nproc_scalar_to_spin,
                int(params_xpure.nproc_scalar_to_spin // params_xpure.nproc_per_node),
                params_xpure.nproc_scalar_to_spin,
            ),
            file=f,
        )
        print("	wait", file=f)
        print("else", file=f)
        print('	echo "Go fast - skip scalar2spin"', file=f)
        print("fi", file=f)

        print(" ", file=f)
        print(
            "################################",
            file=f,
        )
        print("# X2pure", file=f)
        print(
            "################################",
            file=f,
        )
        print("cd ${OUTPUTMLL}", file=f)

        print(" ", file=f)
        print(
            "################################",
            file=f,
        )
        print("# Mode-mode mixing matrix", file=f)
        print(
            "################################",
            file=f,
        )
        print("cat > createMll.par << EOF", file=f)

        print(" ", file=f)
        print("##### MODE #######", file=f)
        print("# 0 : Standard formalism", file=f)
        print("# 1 : Pure formalism", file=f)
        print("# 2 : Hybrid formalism", file=f)
        print("##############", file=f)
        print("mode = ${MODE_XPURE}", file=f)

        print(" ", file=f)
        print("###### SETUP #####", file=f)
        print("nside = ${NSIDE_USER}", file=f)
        print("lmax = ${LMAX_USER}", file=f)
        print("nmask = ${NMASK_USER}", file=f)

        print(" ", file=f)
        print("EOF", file=f)

        print(" ", file=f)
        print("for (( i=1; i<=${NMASK_USER}; i++)); do", file=f)
        print("	ind=$(($i - 1))", file=f)
        print("	cat >> createMll.par << EOF", file=f)

        print(" ", file=f)
        print("maskfile${i}_T  = ${output_spin0_I1}", file=f)

        print(" ", file=f)
        print("maskfile${i}_E_spin0 = ${output_spin0_P1}", file=f)
        print("maskfile${i}_E_spin1 = ${output_spin1_P1}", file=f)
        print("maskfile${i}_E_spin2 = ${output_spin2_P1}", file=f)

        print(" ", file=f)
        print("maskfile${i}_B_spin0 = ${output_spin0_P1}", file=f)
        print("maskfile${i}_B_spin1 = ${output_spin1_P1}", file=f)
        print("maskfile${i}_B_spin2 = ${output_spin2_P1}", file=f)

        print(" ", file=f)
        print("mllfile_TT_TT_${i} = ${OUTPUTMLL}/mll_TT_TT_BinMask${i}.fits", file=f)

        print(" ", file=f)
        print("mllfile_EE_EE_${i} = ${OUTPUTMLL}/mll_spinEE_EE_pcg${i}.fits", file=f)
        print("mllfile_EE_BB_${i} = ${OUTPUTMLL}/mll_spinEE_BB_pcg${i}.fits", file=f)
        print("mllfile_EE_EB_${i} = ${OUTPUTMLL}/mll_spinEE_EB_pcg${i}.fits", file=f)
        print("mllfile_BB_BB_${i} = ${OUTPUTMLL}/mll_spinBB_BB_pcg${i}.fits", file=f)
        print("mllfile_BB_EE_${i} = ${OUTPUTMLL}/mll_spinBB_EE_pcg${i}.fits", file=f)
        print("mllfile_BB_EB_${i} = ${OUTPUTMLL}/mll_spinBB_EB_pcg${i}.fits", file=f)

        print(" ", file=f)
        print("mllfile_TE_TE_${i} = ${OUTPUTMLL}/mll_spinTE_TE_pcg${i}.fits", file=f)
        print("mllfile_TE_TB_${i} = ${OUTPUTMLL}/mll_spinTE_TB_pcg${i}.fits", file=f)
        print("mllfile_TB_TE_${i} = ${OUTPUTMLL}/mll_spinTB_TE_pcg${i}.fits", file=f)
        print("mllfile_TB_TB_${i} = ${OUTPUTMLL}/mll_spinTB_TB_pcg${i}.fits", file=f)

        print(" ", file=f)
        print("mllfile_EB_EB_${i} = ${OUTPUTMLL}/mll_spinEB_EB_pcg${i}.fits", file=f)
        print("mllfile_EB_EE_${i} = ${OUTPUTMLL}/mll_spinEB_EE_pcg${i}.fits", file=f)
        print("mllfile_EB_BB_${i} = ${OUTPUTMLL}/mll_spinEB_BB_pcg${i}.fits", file=f)

        print(" ", file=f)
        print("EOF", file=f)

        print(" ", file=f)
        print("done", file=f)

        print(" ", file=f)
        print('if [ "${FAST}" -eq "0" ]', file=f)
        print("        then", file=f)
        txt = "	time srun -N {} -n {} "\
            + "${{BINDIR}}/x2pure_create_mll createMll.par"
        print(
            txt.format(
                int(params_xpure.nproc_mll // params_xpure.nproc_per_node),
                params_xpure.nproc_mll,
            ),
            file=f,
        )

        print("	rm -f createMll.par", file=f)
        print('elif [ "${FAST}" -eq "2" ]', file=f)
        print("        then", file=f)
        txt = "        time srun -N {} -n {} "\
            + "${{BINDIR}}/x2pure_create_mll createMll.par"
        print(
            txt.format(
                int(params_xpure.nproc_mll // params_xpure.nproc_per_node),
                params_xpure.nproc_mll,
            ),
            file=f,
        )
        print("        rm -f createMll.par", file=f)
        print("else", file=f)
        print('	echo "Go fast - skip x2pure_create_mll"', file=f)
        print("fi", file=f)

        print(" ", file=f)
        print(
            "#################################",
            file=f,
        )
        print("# X2pure", file=f)
        print(
            "################################",
            file=f,
        )
        print("for (( n=0; n<${NSIMU_USER}; n++ )); do", file=f)

        print(" ", file=f)
        print(
            '    echo "************************ simu $n ************************"',
            file=f,
        )

        print(" ", file=f)
        print("    num=$n", file=f)

        print(" ", file=f)
        print("    #CREATE PARAMETER FILE", file=f)
        print("    cat > xpure.par << _EOF_", file=f)

        print(" ", file=f)
        print("mode = ${MODE_XPURE}", file=f)

        print(" ", file=f)
        print("nside = ${NSIDE_USER}", file=f)
        print("nmaps = ${NMAPS_USER}", file=f)
        print("nmasks = ${NMASK_USER}", file=f)

        print(" ", file=f)
        print("_EOF_", file=f)

        print(" ", file=f)
        print("    for ((j=1; j<=${NMAPS_USER}; j++)); do", file=f)

        print(" ", file=f)
        print("	cat >> xpure.par << _EOF_", file=f)

        print(" ", file=f)
        print("bellfile${j} = ${THEORYDIR}/${BEAMFILE}", file=f)
        print("mapfile${j} = ${MAPDIR}/${MAPS[${j}]}", file=f)

        print(" ", file=f)
        print("_EOF_", file=f)

        print(" ", file=f)
        print("    done", file=f)

        print(" ", file=f)
        print("    cat >> xpure.par << _EOF_", file=f)

        print(" ", file=f)
        print("lmaxSim = ${LMAX_USER}", file=f)

        print(" ", file=f)
        print("_EOF_", file=f)

        print(" ", file=f)
        print("    for(( i=1; i<=${NMASK_USER}; i++)); do", file=f)
        print("	ind=$(($i - 1))", file=f)
        print("	cat >> xpure.par << EOF", file=f)

        print(" ", file=f)
        print("maskfile${i}_T  = ${output_spin0_I1}", file=f)

        print(" ", file=f)
        print("maskfile${i}_E_spin0 = ${output_spin0_P1}", file=f)
        print("maskfile${i}_E_spin1 = ${output_spin1_P1}", file=f)
        print("maskfile${i}_E_spin2 = ${output_spin2_P1}", file=f)

        print(" ", file=f)
        print("maskfile${i}_B_spin0 = ${output_spin0_P1}", file=f)
        print("maskfile${i}_B_spin1 = ${output_spin1_P1}", file=f)
        print("maskfile${i}_B_spin2 = ${output_spin2_P1}", file=f)

        print(" ", file=f)
        print("mllfile_TT_TT_${i} = ${OUTPUTMLL}/mll_TT_TT_BinMask${i}.fits", file=f)

        print(" ", file=f)
        print("mllfile_EE_EE_${i} = ${OUTPUTMLL}/mll_spinEE_EE_pcg${i}.fits", file=f)
        print("mllfile_EE_BB_${i} = ${OUTPUTMLL}/mll_spinEE_BB_pcg${i}.fits", file=f)
        print("mllfile_EE_EB_${i} = ${OUTPUTMLL}/mll_spinEE_EB_pcg${i}.fits", file=f)
        print("mllfile_BB_BB_${i} = ${OUTPUTMLL}/mll_spinBB_BB_pcg${i}.fits", file=f)
        print("mllfile_BB_EE_${i} = ${OUTPUTMLL}/mll_spinBB_EE_pcg${i}.fits", file=f)
        print("mllfile_BB_EB_${i} = ${OUTPUTMLL}/mll_spinBB_EB_pcg${i}.fits", file=f)

        print(" ", file=f)
        print("mllfile_TE_TE_${i} = ${OUTPUTMLL}/mll_spinTE_TE_pcg${i}.fits", file=f)
        print("mllfile_TE_TB_${i} = ${OUTPUTMLL}/mll_spinTE_TB_pcg${i}.fits", file=f)
        print("mllfile_TB_TE_${i} = ${OUTPUTMLL}/mll_spinTB_TE_pcg${i}.fits", file=f)
        print("mllfile_TB_TB_${i} = ${OUTPUTMLL}/mll_spinTB_TB_pcg${i}.fits", file=f)

        print(" ", file=f)
        print("mllfile_EB_EB_${i} = ${OUTPUTMLL}/mll_spinEB_EB_pcg${i}.fits", file=f)
        print("mllfile_EB_EE_${i} = ${OUTPUTMLL}/mll_spinEB_EE_pcg${i}.fits", file=f)
        print("mllfile_EB_BB_${i} = ${OUTPUTMLL}/mll_spinEB_BB_pcg${i}.fits", file=f)

        print(" ", file=f)
        print("EOF", file=f)
        print("    done", file=f)

        print(" ", file=f)
        print("    cat >> xpure.par << _EOF_", file=f)
        print("noise_biasT_1 = 0.", file=f)
        print("noise_biasT_2_2 = 0.", file=f)
        print("noise_biasT_1_2 = 0.", file=f)

        print(" ", file=f)
        print("noise_biasP_1 = 0.", file=f)
        print("noise_biasP_2_2 = 0.", file=f)
        print("noise_biasP_1_2 = 0.", file=f)

        print(" ", file=f)
        print("bintab = ${THEORYDIR}/${BINTAB}", file=f)
        print("#mask_list = ${INPUT}/bin2mask_43bins.fits", file=f)
        print("_EOF_", file=f)

        print(" ", file=f)
        print("cat >> xpure.par << _EOF_", file=f)
        print("pseudofile = ${OUTPUTCL}/pseudopure_pcg_${name}_$num", file=f)
        print("cellfile = ${OUTPUTCL}/cellpure_pbear_${name}_$num", file=f)
        print("lmax = ${LMAX_USER}", file=f)
        print("_EOF_", file=f)
        print("    #RUN", file=f)
        print(
            "    time srun -N {} -n {} ${{BINDIR}}/x2pure xpure.par".format(
                int(params_xpure.nproc_xpure // params_xpure.nproc_per_node),
                params_xpure.nproc_xpure,
            ),
            file=f,
        )

        print(" ", file=f)
        print("done", file=f) def write_maps_a_la_xpure(OutputSkyMap, name_out, output_path):
    """
    Grab sky data from OutputSkyMap and write them into files readable by
    the software xpure.

    Parameters
    ----------
    OutputSkyMap : OutputSkyMap instance
        Instance of OutputSkyMap containing map data.
    name_out : string
        File name (.fits) where to store the data.
    output_path : string
        Folder where to put the data.

    """
    safe_mkdir(os.path.join(output_path, name_out))

    fits_I = np.zeros((12 * OutputSkyMap.nside * OutputSkyMap.nside))
    fits_Q = np.zeros((12 * OutputSkyMap.nside * OutputSkyMap.nside))
    fits_U = np.zeros((12 * OutputSkyMap.nside * OutputSkyMap.nside))

    I, Q, U = OutputSkyMap.get_IQU()
    fits_I[OutputSkyMap.obspix] = I
    fits_Q[OutputSkyMap.obspix] = Q
    fits_U[OutputSkyMap.obspix] = U

    full_path = os.path.join(output_path, name_out, "IQU_{}.fits".format(name_out))

    write_healpix_cmbmap(
        full_path,
        data=[fits_I, fits_Q, fits_U],
        fits_IDL=False,
        coord=None,
        colnames=["I", "Q", "U"],
        partial=False,
        nest=False,
    )

    del fits_I, fits_Q, fits_U def write_weights_a_la_xpure(OutputSkyMap, name_out, output_path, epsilon, HWP=False):
    """
    Grab weight and mask from OutputSkyMap and write them into files
    readable by the software xpure.

    Parameters
    ----------
    OutputSkyMap : OutputSkyMap instance
        Instance of OutputSkyMap containing map data.
    name_out : string
        File name (.fits) where to store the data.
    output_path : string
        Folder where to put the data.
    epsilon : float
        Threshold for selecting the pixels. 0 <= epsilon < 1/4.
        The higher the more selective.
    HWP : bool
        If True, use demodulation syntax for the weights (w0, w4).
        Default is False (pair difference syntax: w, cc, ss, cs)

    """
    safe_mkdir(os.path.join(output_path, name_out))

    # Intensity
    weight = np.zeros((12 * OutputSkyMap.nside * OutputSkyMap.nside))

    if not HWP:
        weight[OutputSkyMap.obspix] = OutputSkyMap.w
    else:
        weight[OutputSkyMap.obspix] = OutputSkyMap.w0

    full_path = os.path.join(output_path, name_out, "Iw_{}.fits".format(name_out))
    write_healpix_cmbmap(
        full_path,
        data=weight,
        fits_IDL=False,
        coord=None,
        colnames=["I"],
        partial=False,
        nest=False,
    )

    binary = np.zeros((12 * OutputSkyMap.nside * OutputSkyMap.nside))
    mask = np.where((weight > 0))[0]
    binary[mask] = 1.0
    full_path = os.path.join(
        output_path, name_out, "Iw_{}_norm.fits".format(name_out)
    )
    write_healpix_cmbmap(
        full_path,
        data=binary,
        fits_IDL=False,
        coord=None,
        colnames=["I"],
        partial=False,
        nest=False,
    )

    # Polarisation
    weight = np.zeros((12 * OutputSkyMap.nside * OutputSkyMap.nside))
    if not HWP:
        weight[OutputSkyMap.obspix] = qu_weight_mineig(
            OutputSkyMap.cc, OutputSkyMap.cs, OutputSkyMap.ss, epsilon
        )
    else:
        weight[OutputSkyMap.obspix] = OutputSkyMap.w4

    full_path = os.path.join(output_path, name_out, "Pw_{}.fits".format(name_out))
    write_healpix_cmbmap(
        full_path,
        data=weight,
        fits_IDL=False,
        coord=None,
        colnames=["P"],
        partial=False,
        nest=False,
    )

    binary = np.zeros((12 * OutputSkyMap.nside * OutputSkyMap.nside))
    mask = np.where((weight > 0))[0]
    binary[mask] = 1.0
    full_path = os.path.join(
        output_path, name_out, "Pw_{}_norm.fits".format(name_out)
    )
    write_healpix_cmbmap(
        full_path,
        data=binary,
        fits_IDL=False,
        coord=None,
        colnames=["P"],
        partial=False,
        nest=False,
    )
from mpi4py import MPI
from s4cmb.input_sky import HealpixFitsMap
from s4cmb.instrument import Hardware
from s4cmb.scanning_strategy import ScanningStrategy
from s4cmb.tod import TimeOrderedDataPairDiff
from s4cmb.tod import OutputSkyMap
from s4cmb.systematics import inject_crosstalk_inside_SQUID
from s4cmb.config_s4cmb import import_string_as_module
from s4cmb.xpure import create_batch
from s4cmb.xpure import write_maps_a_la_xpure
from s4cmb.xpure import write_weights_a_la_xpure
import commands
import os
import healpy as hp
import numpy as np
import argparse
                            no_quleak=params.no_quleak)

    ## Initialise our instrument
    inst = Hardware(ncrate=params.ncrate,
                    ndfmux_per_crate=params.ndfmux_per_crate,
                    nsquid_per_mux=params.nsquid_per_mux,
                    npair_per_squid=params.npair_per_squid,
                    fp_size=params.fp_size,
                    fwhm=params.fwhm,
                    beam_seed=params.beam_seed,
                    projected_fp_size=params.projected_fp_size,
                    pm_name=params.pm_name,
                    type_hwp=params.type_hwp,
                    freq_hwp=params.freq_hwp,
                    angle_hwp=params.angle_hwp,
                    verbose=params.verbose)

    ## Initialize our scanning strategy
    scan = ScanningStrategy(nces=params.nces,
                            start_date=params.start_date,
                            telescope_longitude=params.telescope_longitude,
                            telescope_latitude=params.telescope_latitude,
                            telescope_elevation=params.telescope_elevation,
                            name_strategy=params.name_strategy,
                            sampling_freq=params.sampling_freq,
                            sky_speed=params.sky_speed,
                            ut1utc_fn=params.ut1utc_fn,
                            language=params.language)
    scan.run()

    ## Let's now generate our TOD from our input sky, instrument,
    ## and scanning strategy.
    if params.verbose:
        print("Proc [{}] doing scans".format(rank), range(
            rank, scan.nces, size))

    ## Get SQUID and bolo ID
    squid_ids = inst.focal_plane.get_indices('Sq')
    bolo_ids = inst.focal_plane.bolo_index_in_squid

    ## Noise seeds
    state_for_noise = np.random.RandomState(params.array_noise_seed)
    seeds_for_noise = state_for_noise.randint(0, 1e6, scan.nces)
    for pos_CES, CESnumber in enumerate(range(rank, scan.nces, size)):
        if params.verbose:
            print("Proc [{}] with seeds ".format(rank),
                  seeds_for_noise[CESnumber], seeds_for_noise)
        tod = TimeOrderedDataPairDiff(
            inst, scan, sky_in,
            CESnumber=CESnumber,
            projection=params.projection,
            nside_out=params.nside_out,
            pixel_size=params.pixel_size,
            width=params.width,
            array_noise_level=params.array_noise_level,
            array_noise_seed=seeds_for_noise[CESnumber],
            mapping_perpair=False)

        ## Initialise map containers for each processor
        if pos_CES == 0:
            sky_out_tot = OutputSkyMap(projection=tod.projection,
                                       nside=tod.nside_out,
                                       obspix=tod.obspix,
                                       npixsky=tod.npixsky,
                                       pixel_size=tod.pixel_size)

        ## Scan input map to get TODs
        d = np.array([
            tod.map2tod(det) for det in range(inst.focal_plane.nbolometer)])

        ## Inject crosstalk
        inject_crosstalk_inside_SQUID(d,
                                      squid_ids,
                                      bolo_ids,
                                      radius=args.radius,
                                      mu=args.mu,
                                      sigma=args.sigma,
                                      beta=args.beta,
                                      seed=args.seed)

        ## Project TOD to maps
        tod.tod2map(d, sky_out_tot)

    MPI.COMM_WORLD.barrier()

    ## Coaddition over all processors.
    ## Note that all processors will then have the coadded data.
    ## If you want informations at the level of each CES (or group of),
    ## use instead:
    ## final_map = OutputSkyMap(nside=nside_out, obspix=tod.obspix)
    ## final_map.coadd_MPI(sky_out_tot, MPI=MPI)
    sky_out_tot.coadd_MPI(sky_out_tot, MPI=MPI)

    if rank == 0:
        if params.projection == 'flat':
            name_out = '{}_{}_{}'.format(params.tag,
                                         params.name_instrument,
                                         params.name_strategy)
            sky_out_tot.pickle_me(
                '{}/sim{:03d}_{}.pkl'.format(
                    args.folder_out, args.sim_number, name_out),
                shrink_maps=False, crop_maps=2**12,
                epsilon=0., verbose=False)

        elif params.projection == 'healpix':
            ## Save data on disk into fits file for later use in xpure
            name_out = 'sim{:03d}_{}_{}_{}'.format(
                args.sim_number,
                params.tag,
                params.name_instrument,
                params.name_strategy)

            write_maps_a_la_xpure(sky_out_tot, name_out=name_out,
                                  output_path='xpure/maps')
            write_weights_a_la_xpure(sky_out_tot, name_out=name_out,
                                     output_path='xpure/masks',
                                     epsilon=0.08, HWP=False)

            if args.inifile_xpure is not None:
                ## Import parameters from the user parameter file

                params_xpure = import_string_as_module(args.inifile_xpure)