#!/usr/bin/python
# Copyright (c) 2016-2021 Julien Peloton, Giulio Fabbian.
#
# This file is part of s4cmb
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
"""
Script to simulate time-ordered data generated by a CMB experiment
scanning the sky.

Author: Julien Peloton, peloton@lal.in2p3.fr
        Giulio Fabbian, g.fabbian@sussex.ac.uk
"""
from __future__ import division, absolute_import, print_function

import sys
import os

import numpy as np
import healpy as hp

# Python3 does not have the cPickle module
try:
    import cPickle as pickle
except ImportError:
    import pickle

from numpy.fft import fft, fftfreq, fftshift

from scipy.signal import firwin
from scipy import fftpack

from s4cmb.detector_pointing import Pointing
from s4cmb.detector_pointing import radec2thetaphi
from s4cmb import input_sky
from s4cmb.tod_f import tod_f
from s4cmb.xpure import qu_weight_mineig

d2r = np.pi / 180.0
am2rad = np.pi / 180.0 / 60.0


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,
            )
        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,
                nclouds=self.nclouds,
                f0=self.f0,
                amp_atm=self.amp_atm,
                corrlength=self.corrlength,
                alpha=self.alpha,
                sampling_freq=self.scanning_strategy.sampling_freq,
            )
        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,
            )
        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,
                nclouds=self.nclouds,
                f0=self.f0,
                amp_atm=self.amp_atm,
                corrlength=self.corrlength,
                alpha=self.alpha,
                sampling_freq=self.scanning_strategy.sampling_freq,
            )
        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,
        )

        # boresight pointing systematics
        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
        # TODO: add polarization angle perturbations for the second frequency.
        if self.mode == "dichroic":
            ang_pix2 = (90.0 - self.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):
        """
        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.

        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]).

        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
        index_global, index_local = self.get_pixel_indices(ra, dec)

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

            # Perturbed boresight pointing values
            index_global_perturbed, index_local_perturbed = self.get_pixel_indices(
                ra_perturbed, dec_perturbed
            )

        else:
            ra_perturbed = ra
            dec_perturbed = dec
            # pa_perturbed = pa
            index_global_perturbed = index_global
            # index_local_perturbed = index_local

        if ((self.projection == 'healpix') & (index_local is None)):
            # Using a pointer not to increase memory usage
            index_local = index_global

            # Perturbed boresight pointing values
            # index_local_perturbed = index_global_perturbed

        # 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
            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 differntial 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[index_local == -1] = -1
                    self.point_matrix_input[int(ch / 2)] = index_global_pair
            elif ch % 2 == 0 and self.mapping_perpair:
                self.point_matrix[0] = index_local
                if self.store_pointing_matrix_input:
                    index_global_pair[index_local == -1] = -1
                    self.point_matrix_input[0] = index_global_pair

        # 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:
            # pol_ang2 is None if mode == 'standard'
            pol_ang, pol_ang2 = self.compute_simpolangle(
                ch, pa, polangle_err=False
            )

            if self.perturb_pol_angs:
                do_polangle_err = True
            else:
                do_polangle_err = False

            pol_ang_perturbed, pol_ang2_perturbed = self.compute_simpolangle(
                ch, pa, polangle_err=do_polangle_err
            )

            # Use pa of pair center for tod2map if differential pointing
            # is used. Otherwise the polarization angle of a pair is the
            # same as the top and bottom detectors.
            # TODO: I do not understand this try/except.
            # TODO: remove the bare exception!!!
            try:
                pol_ang_pair, pol_ang2_pair = self.compute_simpolangle(
                    ch, pa_pair, polangle_err=False
                )
            except UnboundLocalError:
                pol_ang_pair = pol_ang
                pol_ang2_pair = pol_ang2

            # 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

            # 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]\
                + self.HealpixFitsMap.Q[index_global] * np.cos(2 * pol_ang)\
                + sign * self.HealpixFitsMap.U[index_global] * np.sin(2 * pol_ang)\
                + noise
            ts1 = ts1 * norm

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

            elif self.mode == "dichroic":
                # 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

                # 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_perturbed]\
                    + self.HealpixFitsMap.Q2[index_global_perturbed]\
                    * np.cos(2 * pol_ang2_perturbed)\
                    + sign\
                    * self.HealpixFitsMap.U2[index_global_perturbed]\
                    * np.sin(2 * pol_ang2_perturbed)\
                    + noise2
                ts2 = ts2 * norm
                return np.array([ts1, ts2])

        else:
            ts1 = norm * (
                self.HealpixFitsMap.I[index_global_perturbed] + noise
            )
            if self.mode == "standard":
                return ts1
            elif self.mode == "dichroic":
                ts2 = norm * (
                    self.HealpixFitsMap.I2[index_global_perturbed] + 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:
            print("gdep!")
            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 TimeOrderedDataDemod(TimeOrderedDataPairDiff):
    """ Class to """

    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,
        mapping_perpair=False,
        mode="standard",
        verbose=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.
        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.
        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.

        Examples
        ----------
        HEALPIX: Test the routines MAP -> TOD -> MAP.
        >>> inst, scan, sky_in = load_fake_instrument()
        >>> tod = TimeOrderedDataDemod(inst, scan, sky_in,
        ...     CESnumber=0, projection='healpix', mapping_perpair=True)
        >>> d = np.array([tod.map2tod(det) for det in range(2)])
        >>> d = tod.demodulate_timestreams(d)
        >>> m = OutputSkyMap(projection=tod.projection,
        ...     nside=tod.nside_out, obspix=tod.obspix, demodulation=True)
        >>> tod.tod2map(d, m)

        Same with dichroic detectors
        >>> inst, scan, sky_in = load_fake_instrument(fwhm_in2=1.8)
        >>> tod = TimeOrderedDataDemod(inst, scan, sky_in, mode='dichroic',
        ...     CESnumber=0, projection='healpix', mapping_perpair=True)
        >>> d = np.array([tod.map2tod(det) for det in range(2)])
        >>> d1 = d[:, 0]
        >>> d2 = d[:, 1]
        >>> d1 = tod.demodulate_timestreams(d1)
        >>> d2 = tod.demodulate_timestreams(d2)
        >>> m1 = OutputSkyMap(projection=tod.projection,
        ...     nside=tod.nside_out, obspix=tod.obspix, demodulation=True)
        >>> m2 = OutputSkyMap(projection=tod.projection,
        ...     nside=tod.nside_out, obspix=tod.obspix, demodulation=True)
        >>> tod.tod2map(d1, m1, frequency_channel=1)
        >>> tod.tod2map(d2, m2, frequency_channel=2)
        """
        TimeOrderedDataPairDiff.__init__(
            self,
            hardware,
            scanning_strategy,
            HealpixFitsMap,
            CESnumber,
            projection=projection,
            nside_out=nside_out,
            pixel_size=pixel_size,
            width=width,
            cut_pixels_outside=cut_pixels_outside,
            array_noise_level=array_noise_level,
            array_noise_seed=array_noise_seed,
            array_noise_level2=array_noise_level2,
            array_noise_seed2=array_noise_seed2,
            mapping_perpair=mapping_perpair,
            mode=mode,
            verbose=verbose,
        )

        # Prepare the demodulation of timestreams
        self.dm = Demodulation(
            self.hardware.half_wave_plate.freq_hwp,
            self.scanning_strategy.sampling_freq,
            self.hwpangle,
            self.verbose,
        )

        # Prepare the filters use for the demodulation
        self.dm.prepfilter([4], [1.9])
        self.dm.prepfftedfilter(nt=self.nsamples)

    def demodulate_timestreams(self, ts):
        """
        Perform the demodulation of timestreams, that is split full timestream
        according to their HWP frequency dependency:
        d --> d0 (intensity) and d4r + j * d4i (polarisation).

        Careful, memory requirement explodes here! You go from an array of size
        (ndet, nbolometer) to (ndet, 3, nbolometer), that is you need 3x more
        memory.

        Parameters
        ----------
        ts : array of size (ndet, nbolometer)
            Array of timestreams.

        Returns
        ----------
        newts : array of size (ndet, 3, nbolometer)
        """
        outshape = (ts.shape[0], 3, ts.shape[1])
        newts = np.zeros(outshape)

        self.dm.b = ts
        # dm.b.copy()
        self.dm.br = ts

        # Do temperature
        self.dm.demod(0)
        newts[:, 0, :] = self.dm.bm.real

        # Do 4f component (effectively polarisation)
        self.dm.demod(4)
        newts[:, 1, :] = self.dm.bm.real
        newts[:, 2, :] = self.dm.bm.imag

        return newts


class Demodulation:
    """ Class to handle demodulation of timestreams """

    def __init__(self, hwp_freq, sampling_freq, hwp_angles, verbose=False):
        """
        Demodulation of timestreams means here that we estimate I, Q and U
        from a single detector timstream. That is possible realising that I
        is the 0f component while Q and U are 4f components wrt the HWP angle.

        The demodulation comports two parts. The first one is to prepare the
        low pass filter and band pass filter that will be used to split out
        I and Q&U. The second part is to rotate fields (effectively applying
        a complex exp) and convolving them with the filters.

        Parameters
        ----------
        hwp_freq : float
            HWP spin frequency in Hz.
        sampling_freq : float
            Detector sampling frequency in Hz.
        hwp_angles : 1D array
            Array containing the hwp angles sampled at the detector
            sampling frequency. In radians.

        For examples, see TimeOrderedDataDemod.

        """
        self.hwp_freq = hwp_freq
        self.sampling_freq = sampling_freq
        self.hwp_angles = hwp_angles
        self.verbose = verbose

        # In Hz
        med = np.median(np.diff(self.hwp_angles))
        self.speed = med * self.sampling_freq / (2 * np.pi)

        # Nyquist frequency at half the sampling
        self.nyq = self.sampling_freq / 2.0

    def prepfilter(self, modes, bands=None, numtaps=None, relative=True):
        """
        Prepare the filters used for the demodulation (real space).

        Parameters
        ----------
        modes : 1D array of int
            This is where you specify which multiplier of hwp_freq you care
            about. For example if you want the 4f component of your
            timestream, modes = [4]. Intensity (0) is never specified.
            For the moment, only 4f is implemented.
        bands : 1D array of int
            Band width around your frequency of interest [Hz]. If not
            specified, it is equal to the speed of the HWP.
        numtaps : int, optional
            Length of the filter (number of coefficients, i.e. the filter
            order + 1).  `numtaps` must be even if a passband includes the
            Nyquist frequency.
        relative : bool
            If True, bands and HWP frequency are relative to the HWP speed.
            Default is True.

        """
        if numtaps is None:
            if self.nyq < 100.0 / 6:
                numtaps = 255
            else:
                numtaps = 1023

        self.modes = np.array([modes]).flatten()

        if bands is None:
            self.bands = np.ones(self.modes.shape) * self.speed
        else:
            self.bands = np.array([bands]).flatten()
            if self.bands.size == 1:
                self.bands = np.tile(self.bands, self.modes.size)

            if relative:
                self.bands *= self.speed
                self.hwp_freq *= self.speed

        bad = self.modes * self.speed + self.bands > self.nyq

        if bad.any():
            print("You need freq = modes * speed + bands < nyquist")
            print("MODES", self.modes)
            print("SPEED", self.speed)
            print("BANDS", self.bands)
            print("FREQ", self.modes * self.speed + self.bands)
            print("NYQUIST", self.nyq)
            print("mode", self.modes[bad], "are over nyq.")
            print(
                "You need either to increase the sampling frequency of",
                "your detectors, or decrease the spin frequency of the HWP.",
            )
            self.modes = self.modes[~bad]
            self.bands = self.bands[~bad]

        self.lpf0 = firwin(numtaps, self.hwp_freq, nyq=self.nyq)

        # Construct the filters
        self.lpfs = []
        self.bpfs = []
        for mode, band in zip(self.modes, self.bands):
            self.lpfs.append(firwin(numtaps, band, nyq=self.nyq))
            bpfreal = firwin(
                numtaps,
                [self.speed * mode - band, self.speed * mode + band],
                nyq=self.nyq,
                pass_zero=False,
            )

            fbpf = fftpack.fft(bpfreal)
            fbpf[: int((numtaps + 1) / 2)] = 0.0
            self.bpfs.append(fftpack.ifft(fbpf))

        self.bpfs = np.array(self.bpfs)
        self.lpfs = np.array(self.lpfs)
        self.numtaps = numtaps

    def prepfftedfilter(self, nt):
        """
        Prepare the filters used for the demodulation (Fourier space).

        Parameters
        ----------
        nt : int
            Length of the timestreams to filter (number of time samples).
        """
        n = self.lpf0.size

        fftsize = int(2 ** np.ceil(np.log2(nt + 3 * n - 1)))
        self.flpf0 = fftpack.fft(self.lpf0, fftsize)

        self.fbpfs = []
        for bpf in self.bpfs:
            self.fbpfs.append(fftpack.fft(bpf, fftsize))
        self.fbpfs = np.asarray(self.fbpfs)

        self.flpfs = []
        for lpf in self.lpfs:
            self.flpfs.append(fftpack.fft(lpf, fftsize))
        self.flpfs = np.asarray(self.flpfs)

    def demod(self, mode, bpf=True, lpf=True):
        """
        Perform the demodulation of timestreams.

        TODO: self.bm and self.br are weirdly set. Get rid of self,
        put br in args are ts, and output bm.

        Parameters
        ----------
        mode : int
            mode=0 is intensity, mode=4 is for polarisation (complex).
        bpf : bool, optional
            If True, perform a Band Pass Filtering. Default is True.
        lpf : bool, optional
            If True, perform a Low Pass Filtering. Default is True.
        """
        # Temperature or not understood!
        if mode not in self.modes:
            if mode == 0:
                # if (self.br == 0.).all():
                #     self.bm = np.zeros(self.br.shape, dtype=self.br.dtype)
                #     return self.bm

                if hasattr(self, "flpf0"):
                    flpf0 = self.flpf0
                else:
                    flpf0 = None

                self.bm = convolvefilter(self.br, self.lpf0, flpf0)
                return self.bm
            else:
                print("Filters for", mode, "is not prepared!")
                return None

        # Polarisation
        else:
            i = np.where(self.modes == mode)[0][0]

        if hasattr(self, "bm"):
            del self.bm

        # if (self.br == 0.).all():
        #     self.bm = np.zeros(self.br.shape, dtype=self.br.dtype) + 0.j
        #     return self.bm

        if bpf:
            bpf = self.bpfs[i]
            try:
                fbpf = self.fbpfs[i]
            except (AttributeError, IndexError):
                fbpf = None
        else:
            bpf = None
        if lpf:
            lpf = self.lpfs[i]
            try:
                flpf = self.flpfs[i]
            except (AttributeError, IndexError):
                flpf = None
        else:
            lpf = None

        self.bm = demod(
            self.br,
            np.exp(1.0j * mode * self.hwp_angles),
            bpf=bpf,
            lpf=lpf,
            fbpf=fbpf,
            flpf=flpf,
        )

        return self.bm


def demod(x, e, bpf=None, lpf=None, fbpf=None, flpf=None):
    """
    Actual demodulation scheme for a timestream vector x.
    First multiply x by 2*e, and then apply Low Pass Filter and/or Band Pass
    Filter to 2*e*x to isolate the desired part of the signal.

    Parameters
    ----------
    x : 1D array
        Your timestream.
    e : 1D array
        Vector containing the value of a complex exponential: exp(i * theta)
    bpf : 1D array, optional
        Coefficients of a Band Pass Filter to apply to the data.
        See scipy.signal.firwin
    lpf : 1D array, optional
        Coefficients of a Low Pass Filter to apply to the data.
        See scipy.signal.firwin
    fbpf : float, optional
        FFTed coefficients of the Band Pass Filter to apply to the data
    flpf : float, optional
        FFTed coefficients of the Low Pass Filter to apply to the data

    Returns
    ----------
    u : 1D array
        Demodulated part of x.

    For examples, see TimeOrderedDataDemod.

    """
    u = np.ones(list(x.shape[:-1]) + [1], x.dtype) * 2.0 * e
    if bpf is not None:
        u *= convolvefilter(x, bpf, fbpf)
    else:
        u *= x
    if lpf is not None:
        u = convolvefilter(u, lpf, flpf)
    return u


def convolvefilter(x, f, ff=None, isreal=False):
    """
    Apply filter f to a timestream x (convolution)

    Parameters
    ----------
    x : 1D array
        Your timestream vector.
    f : 1D array
        Coefficients of a finite impulse response filter.
        See scipy.signal.firwin
    ff : 1D array, optional
        FFT of the coefficients of a finite impulse response filter.
        See scipy.fftpack.fft
    ireal : bool, optional
        If True, return only the real part of the filtered x.

    Returns
    ----------
    u : 1D array
        Filtered x.

    For examples, see TimeOrderedDataDemod.
    """
    # since f might not be numpy.ndarray
    assert np.ndim(f) == 1

    # Get shapes
    init_shape = x.shape
    n = f.size
    nt = init_shape[-1]

    # How many FFTs to perform
    fftsize = int(2 ** np.ceil(np.log2(nt + 3 * n - 1)))
    fftslice = (slice(None), slice((3 * n - 1) // 2, (3 * n - 1) // 2 + nt))

    x2 = x.reshape(-1, nt)
    u = np.zeros((init_shape[0], fftsize), (np.ones(1, x.dtype) + 0.0j).dtype)
    (u[:, :n].T)[:] = x2[:, 0]
    u[:, n: nt + n] = x2
    (u[:, nt + n: nt + n * 2].T)[:] = x2[:, -1]

    if ff is None:
        ff = fftpack.fft(f, fftsize)

    for u1 in u:
        fftpack.fft(u1, fftsize, overwrite_x=True)
        u1 *= ff
        fftpack.ifft(u1, fftsize, overwrite_x=True)

    # Condition to return the real part
    cond1 = x.dtype == float or x.dtype == np.float32
    cond2 = x.dtype == int or x.dtype == np.int32
    cond3 = f.dtype == float or f.dtype == np.float32
    cond4 = f.dtype == int or f.dtype == np.int32
    if isreal or (cond1 or cond2) and (cond3 or cond4):
        u = u.real

    return u[fftslice].reshape(init_shape)


class WhiteNoiseGenerator:
    """ Class to handle white noise """

    def __init__(
        self, array_noise_level, ndetectors, ntimesamples, array_noise_seed
    ):
        """
        This class is used to simulate time-domain noise.
        Usually, it is used in combination with map2tod to insert noise
        on-the-fly while scanning an input CMB map.

        Parameters
        ----------
        array_noise_level : float
            Noise level for a detector in [u]K.sqrt(s). WARNING: units has
            to be same as the input map!
        ndetectors : int
            Total number of detectors in the focal plane.
        ntimesamples : float
            Number of time samples per timestream (length of the observation).
        array_noise_seed : int
            Seed used to generate random numbers. From this single seed,
            we generate a list of seeds for all detectors.

        """
        self.array_noise_level = array_noise_level
        self.ndetectors = ndetectors
        self.ntimesamples = ntimesamples

        # Noise level for one detector
        self.detector_noise_level = self.array_noise_level * np.sqrt(
            self.ndetectors
        )

        self.array_noise_seed = array_noise_seed
        state = np.random.RandomState(self.array_noise_seed)
        self.noise_seeds = state.randint(0, 1e6, size=self.ndetectors)

    def simulate_noise_one_detector(self, ch):
        """
        Simulate white noise on-the-fly for one detector.

        Parameters
        ----------
        ch : int
            Index of the detector in the array.

        Returns
        ----------
        vec : 1d array
            Vector of noise of size ntimesamples.
            The level of noise is given by detector_noise_level in uK.sqrt(s).

        Examples
        ----------
        >>> wn = WhiteNoiseGenerator(3000., 2, 4, array_noise_seed=493875)
        >>> ts = wn.simulate_noise_one_detector(0)
        >>> print(ts) #doctest: +NORMALIZE_WHITESPACE
        [ -2185.65609023   5137.21044598  -5407.22292574  11020.59471471]
        """
        state = np.random.RandomState(self.noise_seeds[ch])
        vec = state.normal(size=self.ntimesamples)

        return self.detector_noise_level * vec


class CorrNoiseGenerator(WhiteNoiseGenerator):
    """ """

    def __init__(
        self,
        array_noise_level,
        ndetectors,
        ntimesamples,
        array_noise_seed,
        nclouds=10,
        f0=0.1,
        alpha=-4,
        amp_atm=1e2,
        corrlength=300,
        sampling_freq=8,
    ):
        """
        This class is used to simulate time-domain correlated noise.
        Usually, it is used in combination with map2tod to insert noise
        on-the-fly while scanning an input CMB map.

        We start first by simulating white noise at the level of detector
        timestream, and then we add the correlated part:

        $$\begin{equation}
        n = n^{w} + n^{corr}.
        \end{equation}$$

        The correlated part is based on an empirical model in
        frequency domain (power spectrum density of the timestream):

        $$\begin{equation}
        PSD = A \Big[ 1 + \Big( \dfrac{f}{f_0} \Big)^{\alpha} \Big].
        \end{equation}$$

        where $A$ is a typical amplitude for the fluctuation,
        $f_0$ is the knee frequency at which correlated noise stops
        being significant, and $\alpha$ is the value of the slope
        (typically negative). The correlated part of the timestream is
        then obtained by taking the inverse Fourier transform of this PSD
        (its square root to be precise).

        To take into account the fact that atmosphere is evolving over time,
        we assume a correlation length in time (typically 5min) over which
        the atmosphere signal is constant. Each detector will see the same
        signal, but with a different amplitude. We assume atmosphere is
        unpolarised, so bolometers from the same pair will see exactly the
        same signal. We compute the PSD over this period of time, and
        generate 1/f noise. At the end of this period, we re-generate new
        phases and we draw a new realisation of atmosphere, and so on.

        Parameters
        ----------
        array_noise_level : float
            White noise level for a detector in [u]K.sqrt(s).
            WARNING: units has to be same as the input map!
        ndetectors : int
            Total number of detectors in the focal plane.
        ntimesamples : float
            Number of time samples per timestream (length of the observation).
        array_noise_seed : int
            Seed used to generate random numbers. From this single seed,
            we generate a list of seeds for all detectors.
        nclouds : int, optional
            Number of clouds, that is number of correlated regions
            in the focal plane.
        f0 : float, optional
            Knee frequency in Hz.
        alpha : int or float, optional
            Value of the 1/f slope.
        amp_atm : float, optional
            Typical value of atmospheric fluctuations in [u]K^2.s.
            WARNING: units has to be same as the input map!
        corrlength : int, optional
            Correlation length in time over which atmosphere signal
            is supposed to be constant. Few minutes is typical.
            Units are seconds.
        sampling_freq : float, optional
            Sampling frequency of the detectors in Hz.

        """
        WhiteNoiseGenerator.__init__(
            self, array_noise_level, ndetectors, ntimesamples, array_noise_seed
        )
        self.nclouds = nclouds
        self.alpha = alpha
        self.sampling_freq = sampling_freq
        self.corrlength = int(corrlength * self.sampling_freq)
        self.f0 = f0
        self.amp_atm = amp_atm

        # Bolometers in a pair get the same seed for correlated noise
        self.pixel_noise_seeds = np.repeat(self.noise_seeds[::2], 2)

    def simulate_noise_one_detector(self, ch):
        """
        Simulate correlated noise on-the-fly for one detector.

        Parameters
        ----------
        ch : int
            Index of the detector in the array.

        Returns
        ----------
        vec : 1d array
            Vector of noise of size ntimesamples.
            The level of noise is given by detector_noise_level in uK.sqrt(s).

        Examples
        ----------
        >>> cn = CorrNoiseGenerator(3000., 2, 17000,
        ...     array_noise_seed=493875, nclouds=1, f0=0.5, amp_atm=1.,
        ...     corrlength=300, alpha=-4, sampling_freq=8.)
        >>> ts = cn.simulate_noise_one_detector(0)
        >>> print(ts) #doctest: +NORMALIZE_WHITESPACE
        [ -7536.5882971    -224.58319073 -10795.19644268 ...,
          -5528.66256308  -3161.93996673  -5174.84161989]
        """
        # White noise part
        state = np.random.RandomState(self.noise_seeds[ch])
        vec = state.normal(size=self.ntimesamples)
        wnoise = self.detector_noise_level * vec

        # Correlated part
        state = np.random.RandomState(self.pixel_noise_seeds[ch])
        # amps = 2 * (-0.5 + state.uniform(size=1))
        amps = state.uniform(size=1)

        corrdet = int(self.ndetectors / self.nclouds)
        state = np.random.RandomState(self.array_noise_seed + ch // corrdet)
        phases = 2 * np.pi * state.rand(self.ntimesamples)

        ts_corr = np.zeros(self.ntimesamples)
        for i in range(0, self.ntimesamples, self.corrlength):
            # Check that you have enough samples
            if self.ntimesamples - i < self.corrlength:
                step = self.ntimesamples - i
            else:
                step = self.corrlength

            # Get the PSD and the frequency range
            fs = fftfreq(step, 1.0 / self.sampling_freq)
            psd = np.zeros_like(fs)

            # Avoid zero frequency
            psd[1:] = self.amp_atm * (1 + (fs[1:] / self.f0) ** self.alpha)

            # Get the TOD from the PSD
            ts_corr[i: i + step] = corr_ts(
                PSD=psd, N=step, amp=amps, phase=phases[i: i + step]
            )

        # remove PSD normalisation and add white noise!
        return ts_corr / np.sqrt(self.sampling_freq) + wnoise


def corr_ts(PSD, N, amp, phase):
    """
    Generate a timestream based on its Power Spectrum Density.

    Parameters
    ----------
    PSD : 1d array
        Power Spectrum Density.
    N : int
        Length of the output timestream.
    amp : float
        Scaling factor (between 0 and 1).
    phase : 1d array
        Phase for the IFFT. Must have the same length as the PSD.

    Returns
    ----------
    ts : 1d array
        Timestream based the PSD. Length N, and units sqrt(PSD).
        WARNING: if your PSD is in uk^2.s, you need to normalise
        your timestream by the sampling rate of the detectors:
        ts <- ts/sqrt(sampling_freq).

    Examples
    ----------
    Assuming sampling at 2 Hz.
    >>> sampling = 2.
    >>> state = np.random.RandomState(0)
    >>> PSD = state.rand(10)
    >>> phase = 2 * np.pi * state.rand(10)
    >>> ts = corr_ts(PSD, N=10, amp=1., phase=phase)
    >>> print(ts / np.sqrt(sampling)) #doctest: +NORMALIZE_WHITESPACE
    [ 1.49579944 -2.75993169  0.19454693  0.96278942  0.20053976
      -0.49232085 0.07355544  0.68893301  0.46171998 -0.82563144]
    """
    Nf = len(PSD)
    PSD[0] = 0.0

    A = amp * N * np.sqrt(PSD)

    FFT = A * np.exp(1j * phase[:Nf])

    ts = np.fft.ifft(FFT, n=N)

    return np.real(ts)


def psdts(ts, sample_rate, NFFT=4096):
    """
    Compute the power spectrum (or power spectral density) of a timestream ts.

    Parameters
    ----------
    ts : 1d array
        Detector timestream.
    sample_rate : float
        Sample rate of the detector in Hz.
    NFFT : int, optional
        The number of points.

    Returns
    ----------
    fs : 1d array
        Frequency bin centers in Hz.
    psd : 1d array
        Power spectral density of the timestream.

    Examples
    ----------
    >>> state = np.random.RandomState(0)
    >>> ts = state.rand(100)
    >>> fs, psd = psdts(ts, sample_rate=1)
    >>> print(round(fs[-1], 2), round(psd[-1], 2))
    0.49 0.5
    """
    # Remove the mean
    ts -= np.mean(ts)

    fs, asd = compute_asd(ts, sample_rate=sample_rate, NFFT=NFFT)
    psd = asd ** 2

    return fs, psd


def compute_asd(x, sample_rate, NFFT=2048, is_complex=False):
    """
    Compute the amplitude spectral density (PSD = ASD**2).

    Parameters
    ----------
    x : 1d array
        Quantity for which you want the ASD.
        Typically a detector timestream. Mean of x should have been removed
        prior to passing x.
    sample_rate : float
        Sample rate of x in Hz.
    NFFT : int, optional
        The number of points.
    is_complex : bool, optional

    Returns
    ----------
    fs : 1d array
        Frequency bin centers in Hz.
    asd : 1d array
        Amplitude spectral density of x.
    """
    # Sanity check in the case of no chunking
    maxnfft = x.shape[0] - (x.shape[0] % 2)
    if NFFT > maxnfft or NFFT < 0:
        NFFT = maxnfft

    period = 1.0 / sample_rate

    window = np.blackman(NFFT)
    window_norm = 1.0 / np.average(window ** 2)

    n = len(x)
    # Truncates and rounds down
    nchunks = int((2 * n) / NFFT)
    if is_complex:
        PSD = np.zeros(NFFT)
    else:
        PSD = np.zeros(int(NFFT / 2))

    for i in range(0, nchunks - 1):
        chunk = window * x[i * int(NFFT / 2): (i + 2) * int(NFFT / 2)]

        fch = fft(chunk)
        if not is_complex:
            fch = fch[0: int(NFFT / 2)]

        cf = fch * np.conj(fch)

        if not is_complex:
            cf[1: (int(NFFT / 2) - 1)] *= 2
            cf = cf.real
        PSD += cf

    fs = fftfreq(NFFT, period)
    if not is_complex:
        fs = fs[0: int(NFFT / 2)]

    if nchunks != 1:
        PSD /= nchunks - 1

    # Numpy normalizes inverse transform
    PSD /= NFFT

    # Power per root Hz
    PSD *= 1.0 / sample_rate

    # Remove power reduction from window
    PSD *= window_norm

    if is_complex:
        fs = fftshift(fs)
        PSD = fftshift(PSD)

    # Convert to amplitude/rtHz
    return fs, PSD ** 0.5


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)


class OutputSkyMapIGQU(OutputSkyMap):
    """ Class to handle sky maps generated by tod2map + G deprojection """

    def __init__(
        self,
        projection,
        obspix=None,
        npixsky=None,
        nside=None,
        pixel_size=None,
    ):
        """
        Initialise all maps: weights, projected TOD, and Stokes parameters.
        Suitable if you want to perform a deprojection of a constant
        contribution in the polarisation timestream:
            d^{-} = G + Qcos(2*theta) + U*sin(2theta).
        Work only for pair difference.

        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.
        """
        OutputSkyMap.__init__(
            self,
            projection=projection,
            obspix=obspix,
            npixsky=npixsky,
            nside=nside,
            pixel_size=pixel_size,
            demodulation=False,
        )

        self.to_coadd = self.to_coadd + " dm wm c s"

        # Map and weight of deprojected spurious signal
        self.dm = np.zeros(self.npixsky)
        self.wm = np.zeros(self.npixsky)

        # cosine and sine (weights for cross term GxQ and GxU)
        self.c = np.zeros(self.npixsky)
        self.s = np.zeros(self.npixsky)

    def buildV(self, ipix):
        """Build vector of projected polarisation TOD for a given sky pixel

        It includes G template for deprojection (dm).

        Parameters
        ----------
        ipix : int
            Pixel index (relative to 0 in the observed pixels)

        Returns
        ----------
        vec : 1d array of length 3.
            vector of projected polarisation TOD for a given sky pixel

        Examples
        ----------
        >>> m1 = OutputSkyMapIGQU(projection='healpix',
        ...     nside=16, obspix=np.array(range(12*16**2)))
        >>> V = m1.buildV(1000)
        """
        return np.array([self.dm[ipix], self.dc[ipix], self.ds[ipix]])

    def buildP(self, ipix):
        """Build pixel weight matrix for polarisation.

        The [3x3] matrix contains the weights for G, Q and U components of
        the mapmaking:

        [wm  c  s]   [G]   [dm]
        [c  cc cs]   [Q]   [dc]
        [s  cs ss] * [U] = [ds]

        Parameters
        ----------
        ipix : int
            Pixel index (relative to 0 in the observed pixels)

        Returns
        ----------
        M : 3x3 array
            Pixel weight matrix for polarisation

        Examples
        ----------
        >>> m1 = OutputSkyMapIGQU(projection='healpix',
        ...     nside=16, obspix=np.array(range(12*16**2)))
        >>> M = m1.buildP(1000)
        """
        M = np.zeros((3, 3))
        M[0, :] = np.array([self.wm[ipix], self.c[ipix], self.s[ipix]])
        M[1, :] = np.array([self.c[ipix], self.cc[ipix], self.cs[ipix]])
        M[2, :] = np.array([self.s[ipix], self.cs[ipix], self.ss[ipix]])
        return M

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

        Returns
        ----------
        goodpix : 1D arrray
            Array containing the inverse determinant of the weight matrix for
            each sky pixel.

        Examples
        ----------
        >>> m1 = OutputSkyMapIGQU(projection='healpix',
        ...     nside=16, obspix=np.array(range(12*16**2)))
        >>> m1.set_goodpix()
        """
        # We have 4 components: I, Q, U and G.
        inonzero = [pix for pix in range(self.npixsky) if self.nhit[pix] > 10]
        self.goodpix = np.zeros((self.npixsky))
        self.goodpix[inonzero] = 1

        for ipix in inonzero:
            M = self.buildP(ipix)
            det = np.linalg.det(M)
            if det == 0 or np.linalg.cond(M) > 1e2:
                self.goodpix[ipix] = 0

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

        [wm  c  s]   [G]   [dm]
        [c  cc cs]   [Q]   [dc]
        [s  cs ss] * [U] = [ds]

        Work with pair differenced data only.

        Returns
        ----------
        G : 1d array
            Spurious deprojected polarisation signal. 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).

        Examples
        ----------
        Pair difference data
        >>> m1 = OutputSkyMapIGQU(projection='healpix',
        ...     nside=16, obspix=np.array([0, 1, 2, 3]))
        >>> G, Q, U = m1.get_QU()
        """
        if not hasattr(self, "goodpix"):
            self.set_goodpix()
        inonzero = [
            pix for pix in range(self.npixsky) if self.goodpix[pix] != 0
        ]

        G = np.zeros((self.npixsky))
        Q = np.zeros((self.npixsky))
        U = np.zeros((self.npixsky))

        for ipix in inonzero:
            M = self.buildP(ipix)
            GQU = np.dot(np.linalg.inv(M), self.buildV(ipix))
            G[ipix] = GQU[0]
            Q[ipix] = GQU[1]
            U[ipix] = GQU[2]
        return G, Q, U

    def get_IQU(self):
        """
        Solve for the temperature and polarisation maps from
        projected sum and difference timestream maps and weights plus
        a constant projected spurious polarisation component (G):

        [w 0   0  0]   [I]   [d ]
        [0 wm  c  s]   [G]   [dm]
        [0 c  cc cs]   [Q]   [dc]
        [0 cs cs ss] * [U] = [ds]

        Examples
        ----------
        Pair difference data
        >>> m1 = OutputSkyMapIGQU(projection='healpix',
        ...     nside=16, obspix=np.array([0, 1, 2, 3]))
        >>> I, G, 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).
        G : 1d array
            Spurious deprojected polarisation signal. 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()
        G, Q, U = self.get_QU()
        return I, G, Q, U


def shrink_me(dic, based_on):
    """
    Shrink maps to remove unecessary zeros.
    Maps have to be squared (so work only for flat sky).

    Parameters
    ----------
    dic : dictionary
        Contain the maps to shrink. Maps have to be numpy.ndarray.
    based_on : string
        The key refering as to the map in the dic from which we
        will identify zeros. Typically the weights, or a mask.

    Returns
    ----------
    dic : dictionary
        Input dictionary with shrinked maps.

    Examples
    ----------
    >>> dic = {'I': np.array([1, 2, 3, 4, 5, 6, 7, 8, 9,
    ...     10, 11, 12, 13, 14, 15, 16]),
    ...        'w': np.array([0, 0, 0, 0, 0, 1, 1, 0, 0,
    ...     1, 1, 0, 0, 0, 0, 0]), 'n': 16}
    >>> print(dic['I'].reshape((4, 4))) #doctest: +NORMALIZE_WHITESPACE
    [[ 1  2  3  4]
     [ 5  6  7  8]
     [ 9 10 11 12]
     [13 14 15 16]]
    >>> print(dic['w'].reshape((4, 4))) #doctest: +NORMALIZE_WHITESPACE
    [[0 0 0 0]
     [0 1 1 0]
     [0 1 1 0]
     [0 0 0 0]]

    >>> dic = shrink_me(dic, based_on='w')
    >>> print(dic['I'].reshape((2, 2))) #doctest: +NORMALIZE_WHITESPACE
    [[ 6  7]
     [10 11]]
    """
    assert based_on in dic, KeyError(
        "{} not in input dictionary!".format(based_on)
    )

    npixr = int(len(dic[based_on]) ** 0.5)
    halfnpixr = int(npixr / 2)
    mask = dic[based_on].reshape((npixr, npixr))
    idx = np.where(mask > 0)

    minx, maxx = np.min(idx[0]), np.max(idx[0])
    miny, maxy = np.min(idx[1]), np.max(idx[1])
    dx = maxx - minx
    dy = maxy - miny
    dxy = np.max((dx, dy))
    halfdxy = int(dxy / 2)

    for k in dic.keys():
        # Filter out fields which aren't arrays
        if type(dic[k]) == np.ndarray:
            # Filter out fields which are arrays but not like based_on.
            npixr_loc = int(len(dic[k]) ** 0.5)
            if npixr_loc == npixr:
                dic[k] = np.array(
                    [
                        i[halfnpixr - halfdxy - 1: halfnpixr + halfdxy + 1]
                        for i in dic[k].reshape((npixr, npixr))[
                            halfnpixr - halfdxy - 1: halfnpixr + halfdxy + 1
                        ]
                    ]
                ).flatten()

    return dic


def crop_me(dic, based_on, npix_per_row=2 ** 12):
    """
    Crop maps to a chosen size.
    Maps have to be squared (so work only for flat sky).
    npix_per_row has to be smaller than the current number of pixels per row.

    Parameters
    ----------
    dic : dictionary
        Contain the maps to shrink. Maps have to be numpy.ndarray.
    based_on : string
        The key refering as to the map in the dic from which we
        will identify zeros. Typically the weights, or a mask.
    npix_per_row : int
        New number of pixels per row.

    Returns
    ----------
    dic : dictionary
        Input dictionary with shrinked maps.

    Examples
    ----------
    >>> dic = {'I': np.array([1, 2, 3, 4, 5, 6, 7, 8, 9,
    ...     10, 11, 12, 13, 14, 15, 16]), 'n': 16}
    >>> print(dic['I'].reshape((4, 4))) #doctest: +NORMALIZE_WHITESPACE
    [[ 1  2  3  4]
     [ 5  6  7  8]
     [ 9 10 11 12]
     [13 14 15 16]]

    >>> dic = crop_me(dic, based_on='I', npix_per_row=2)
    >>> print(dic['I'].reshape((2, 2))) #doctest: +NORMALIZE_WHITESPACE
    [[ 6  7]
     [10 11]]
    """
    assert based_on in dic, KeyError(
        "{} not in input dictionary!".format(based_on)
    )

    npixr = int(len(dic[based_on]) ** 0.5)
    halfnpixr = int(npixr / 2)
    halfnpix_per_row = int(npix_per_row / 2)

    for k in dic.keys():
        # Filter out fields which aren't arrays
        if type(dic[k]) == np.ndarray:
            # Filter out fields which are arrays but not like based_on.
            npixr_loc = int(len(dic[k]) ** 0.5)
            if npixr_loc == npixr:
                # Check that we have enough pixels to start with
                assert npixr_loc >= npix_per_row, ValueError(
                    "Map too small to be cropped! ({} vs {})".format(
                        npixr_loc,
                        npix_per_row
                    )
                )

                dic[k] = np.array(
                    [
                        i[
                            halfnpixr - halfnpix_per_row: halfnpixr + halfnpix_per_row
                        ]
                        for i in dic[k].reshape((npixr, npixr))[
                            halfnpixr - halfnpix_per_row: halfnpixr + halfnpix_per_row
                        ]
                    ]
                ).flatten()

    return dic


def partial2full(partial_obs, obspix, nside, fill_with=0.0):
    """
    Reconstruct full sky map from a partial observation and a list of observed
    pixels.

    Parameters
    ----------
    partial_obs : 1d array
        Array containg the values of observed pixels.
    obspix : 1d array
        Array containing the healpix indices of observed pixels.
    nside : int
        The resolution of the map (obspix and partial_obs should have the same
        nside).
    fill_with : optional
        Fill the initial array with `fill_with`. Default is 0.0.

    Returns
    ----------
    fullsky : 1d array
        Full sky map of size 12 * nside**2.

    Examples
    ----------
    >>> nside = 16
    >>> data = np.random.rand(10)
    >>> obspix = np.arange(12 * nside**2, dtype=int)[30:40]
    >>> fullsky = partial2full(data, obspix, nside)
    """
    fullsky = np.zeros(12 * nside ** 2) * fill_with
    fullsky[obspix] = partial_obs
    return fullsky


def build_pointing_matrix(
    ra,
    dec,
    nside_in,
    nside_out=None,
    projection="healpix",
    obspix=None,
    ext_map_gal=False,
    xmin=None,
    ymin=None,
    pixel_size=None,
    npix_per_row=None,
    cut_pixels_outside=True,
):
    """
    Given pointing coordinates (RA/Dec), retrieve the corresponding healpix
    pixel index for a full sky map. This acts effectively as an operator
    to go from time domain to map domain.

    If a list of observed pixel in a sky patch provided (obspix),
    the routines returns also local indices of pixels relative
    to where they are in obspix.
    Note that the indexing 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.

    Parameters
    ----------
    ra : float or 1d array
        RA coordinates of the detector in radian.
    dec : float or 1d array
        Dec coordinates of the detector in radian.
    nside_in : int
        Resolution of the input map.
    nside_out : int
        Resolution for the output map.
    obspix : 1d array, optional
        Array with indices of observed pixels for the sky patch (used to make
        the conversion global indices to local indices).
        Should have been built with nside_out. Default is None.
    ext_map_gal : bool, optional
        If True, perform a rotation of the RA/Dec coordinate to Galactic
        coordinates prior to compute healpix indices. Defaut is False.
    xmin : int, optional
        In flat sky projection, the patch is centered on (x0, y0)
        and has a width=w. We define xmin = x0 - w/2. * np.pi/180.
    ymin : int, optional
        In flat sky projection, the patch is centered on (x0, y0)
        and has a width=w. We define ymin = y0 - w/2. * np.pi/180.
    pixel_size : float, optional
        The pixel size for the output maps if projection=flat. In radian.
    npix_per_row : int, optional
        In flat sky projection, it corresponds to the number of pixel per row.
        In other word, this is the square root of the total number of pixels
        in your square patch.
    cut_pixels_outside : bool, optional
        If True assign -1 to pixels not in obspix. If False, the routine
        crashes if there are pixels outside. Default is True.

    Returns
    ----------
    index_global : float or 1d array
        The input pixels seen labeled as if it was a full sky healpix map.
        To be used for the projection map2tod.
    index_local : None or float or 1d array
        The indices of pixels relative to where they are in obspix. None if
        obspix is not provided. To be used for the projection tod2map.

    Examples
    ----------
    >>> index_global, index_local = build_pointing_matrix(0.0, -np.pi/4,
    ...     nside_in=16)
    >>> print(index_global)
    2592

    >>> index_global, index_local = build_pointing_matrix(
    ...     np.array([0.0, 0.0]), np.array([-np.pi/1000, np.pi/1000]),
    ...     nside_in=16, projection='flat', xmin=-30./60.*np.pi/180.,
    ...     ymin=-30./60.*np.pi/180., pixel_size=15. / 60. * np.pi / 180.,
    ...     npix_per_row=4)
    >>> print(index_local)
    [ 9 11]

    >>> index_global, index_local = build_pointing_matrix(
    ... np.array([0.0, 0.0]), np.array([-np.pi/4, np.pi/4]),
    ...  nside_in=16, nside_out=8, obspix=np.array(range(12*16**2)))
    >>> print(index_global, index_local)
    [2592  420] [624 112]
    """
    if nside_out is None:
        nside_out = nside_in

    theta, phi = radec2thetaphi(ra, dec)
    if ext_map_gal:
        r = hp.Rotator(coord=["C", "G"])
        theta, phi = r(theta, phi)

    index_global = hp.ang2pix(nside_in, theta, phi)

    if projection == "healpix" and obspix is not None:
        index_global_out = hp.ang2pix(nside_out, theta, phi)
        index_local = obspix.searchsorted(index_global_out)
        mask1 = index_local < len(obspix)
        loc = mask1
        loc[mask1] = obspix[index_local[mask1]] == index_global_out[mask1]
        outside_pixels = np.invert(loc)

        # Handling annoying cases.
        if np.sum(outside_pixels) and (not cut_pixels_outside):
            msg = """
            Pixels outside patch boundaries.
            Patch width insufficient. To avoid this,
            increase the parameter width while initialising the TOD
            or set cut_pixels_outside to True to get a cropped map.
            """
            raise ValueError(msg)
        elif np.sum(outside_pixels) and cut_pixels_outside:
            if not ("msg_cut" in globals()):
                global msg_cut
                msg_cut = """
                Pixels outside patch boundaries.
                Your output map will be cropped. To avoid this,
                increase the parameter width while initialising the TOD.
                """
                print(msg_cut)
            index_local[outside_pixels] = -1
        else:
            pass

    elif projection == "flat":
        x, y = input_sky.LamCyl(ra, dec)

        xminmap = xmin - pixel_size / 2.0
        yminmap = ymin - pixel_size / 2.0

        ix = np.int_((x - xminmap) / pixel_size)
        iy = np.int_((y - yminmap) / pixel_size)

        index_local = ix * npix_per_row + iy

        outside_pixels = (
            (ix < 0) | (ix >= npix_per_row) | (iy < 0) | (iy >= npix_per_row)
        )

        # Handling annoying cases.
        if np.sum(outside_pixels) and (not cut_pixels_outside):
            msg = """
            Pixels outside patch boundaries.
            Patch width insufficient. To avoid this,
            increase the parameter width while initialising the TOD
            or set cut_pixels_outside to True to get a cropped map.
            """
            raise ValueError(msg)
        elif np.sum(outside_pixels) and cut_pixels_outside:
            if not ("msg_cut_flat" in globals()):
                global msg_cut_flat
                msg_cut_flat = """
                Pixels outside patch boundaries.
                Your output map will be cropped. To avoid this,
                increase the parameter width while initialising the TOD.
                """
                print(msg_cut_flat)
            index_local[outside_pixels] = -1

        else:
            pass
    else:
        index_local = None

    return index_global, index_local


def load_fake_instrument(
    nside=16, nsquid_per_mux=1, fwhm_in2=None, compute_derivatives=False
):
    """
    For test purposes.
    Create instances of HealpixFitsMap, hardware, and
    scanning_strategy to feed TimeOrderedDataPairDiff in tests.

    Returns
    ----------
    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.
    """
    # Add paths to load modules
    sys.path.insert(0, os.path.realpath(os.path.join(os.getcwd(), ".")))
    sys.path.insert(0, os.path.realpath(os.path.join(os.getcwd(), "s4cmb")))
    from s4cmb.input_sky import HealpixFitsMap
    from s4cmb.instrument import Hardware
    from s4cmb.scanning_strategy import ScanningStrategy

    # Create fake inputs

    # Sky
    sky_in = HealpixFitsMap(
        "s4cmb/data/test_data_set_lensedCls.dat",
        do_pol=True,
        fwhm_in=0.0,
        fwhm_in2=fwhm_in2,
        nside_in=nside,
        map_seed=48584937,
        compute_derivatives=compute_derivatives,
        verbose=False,
        no_ileak=False,
        no_quleak=False,
    )

    # Instrument
    inst = Hardware(
        ncrate=1,
        ndfmux_per_crate=1,
        nsquid_per_mux=nsquid_per_mux,
        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=0.2,
        angle_hwp=0.0,
        verbose=False,
    )
    if fwhm_in2 is not None:
        inst.make_dichroic(fwhm=fwhm_in2)

    # Scanning strategy
    scan = ScanningStrategy(
        nces=2,
        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=8.0,
        sky_speed=0.4,
        language="fortran",
    )
    scan.run()

    return inst, scan, sky_in


def noise_ukam(array_noise_level, fsky, nside, tobs):
    """
    Estimate quickly the noise level in map domain [uK.arcmin]
    given the NET of the telescope, and the observation time.

    Note that if you want to form a noise power-spectrum, you just have to
    perform a change of units: uK.arcmin -> uk.rad -> uk
    N_ell = (noise_ukam() * np.pi / 180 / 60)**2

    Parameters
    ----------
    array_noise_level : float
        Noise level for the whole array in uK.sqrt(s).
    fsky : float
        Fraction of sky observed.
    nside : int
        Resolution of the map.
    tobs : float
        Observation time in second.

    Returns
    ----------
    noise : float
        Level of noise in map domain (RMS per pixel) in uK.arcmin.

    Examples
    ----------
    Level of noise in a nside=2048 map after scanning 5 percent of the sky
    for 1 year.
    >>> noise = noise_ukam(8.1, 0.05, nside=2048, tobs=1. * 365 * 24 * 3600)
    >>> print(round(noise, 2), 'uK.arcmin')
    3.93 uK.arcmin
    """
    noise = np.sqrt(
        array_noise_level ** 2 * hp.nside2npix(nside) * fsky / tobs
    )
    return noise * hp.nside2resol(nside, arcmin=True)


if __name__ == "__main__":
    import doctest

    if np.__version__ >= "1.14.0":
        np.set_printoptions(legacy="1.13")
    doctest.testmod()
