<commit_before>/*
 * observation.cpp
 *
 *  Created on: Apr 29, 2013
 *      Author: <NAME>
 */

#include "slsimlib.h"

#ifdef ENABLE_FITS
#include <CCfits/CCfits>
//#include <CCfits>
#endif

#ifdef ENABLE_FFTW
#include "fftw3.h"
#endif

#include <fstream>
#include <limits>


/** * \brief Creates an observation setup that mimics a known instrument
 *
 */
Observation::Observation(Telescope tel_name)
{
	if (tel_name == Euclid_VIS)
	{
		diameter = 119.;
		transmission = 0.30;
		exp_time = 1800.;
		exp_num = 3;
		back_mag = 22.8;
		ron = 5.;
		seeing = 0.18;
		pix_size = .1/60./60./180.*pi;
	}
	if (tel_name == Euclid_Y)
	{
		diameter = 119.;
		transmission = 0.0961;
		exp_time = 264.;
		exp_num = 3;
		back_mag = 22.57;
		ron = 5.;
		seeing = 0.3;
		pix_size = .3/60./60./180.*pi;
	}
	if (tel_name == Euclid_J)
	{
		diameter = 119.;
		transmission = 0.0814;
		exp_time = 270.;
		exp_num = 3;
		back_mag = 22.53;
		ron = 5.;
		seeing = 0.3;
		pix_size = .3/60./60./180.*pi;
	}
	if (tel_name == Euclid_H)
	{
		diameter = 119.;
		transmission = 0.1692;
		exp_time = 162.;
		exp_num = 3;
		back_mag = 22.59;
		ron = 5.;
		seeing = 0.3;
		pix_size = .3/60./60./180.*pi;
	}
	if (tel_name == KiDS_u)
	{
		diameter = 265.;
		transmission = 0.032;
		exp_time = 1000.;
		exp_num = 5;
		back_mag = 22.93;
		ron = 5.;
		seeing = 1.0;
		pix_size = .2/60./60./180.*pi;
	}
	if (tel_name == KiDS_g)
	{
		diameter = 265.;
		transmission = 0.1220;
		exp_time = 900.;
		exp_num = 5;
		back_mag = 22.29;
		ron = 5.;
		seeing = 0.8;
		pix_size = .2/60./60./180.*pi;
	}
	if (tel_name == KiDS_r)
	{
		diameter = 265.;
		transmission = 0.089;
		exp_time = 1800.;
		exp_num = 5;
		back_mag = 21.40;
		ron = 5.;
		seeing = 0.7;
		pix_size = .2/60./60./180.*pi;
	}
	if (tel_name == KiDS_i)
	{
		diameter = 265.;
		transmission = 0.062;
		exp_time = 1200.;
		exp_num = 5;
		back_mag = 20.64;
		ron = 5.;
		seeing = 1.1;
		pix_size = .2/60./60./180.*pi;
	}

	mag_zeropoint = 2.5*log10(diameter*diameter*transmission*pi/4./hplanck) - 48.6;
	telescope = true;
}

/** *  Creates a custom observation setup with parameters decided by the user.
 *
 * \param diameter Diameter of telescope (in cm) (Collecting area is pi/4.*diameter^2)
 * \param transmission Total transmission of the telescope (/int T(\lambda)/\lambda d\lambda)
 * \param exp_time Total exposure time
 * \param exp_num Number of exposures
 * \param back_mag Flux due to background and/or sky in mag/arcsec^2
 * \param ron Read-out noise in electrons/pixel
 * \param seeing FWHM in arcsecs of the image
 */
Observation::Observation(float diameter, float transmission, float exp_time, int exp_num, float back_mag, float ron, float seeing):
		diameter(diameter), transmission(transmission), exp_time(exp_time), exp_num(exp_num), back_mag(back_mag), ron(ron), seeing(seeing)
		{
			mag_zeropoint = 2.5*log10(diameter*diameter*transmission*pi/4./hplanck) - 48.6;
			telescope = false;
		}

/**  Creates a custom observation setup with parameters decided by the user. Allows for the use of a psf fits image.
 *
 * \param diameter Diameter of telescope (in cm) (Collecting area is pi/4.*diameter^2)
 * \param transmission Total transmission of the telescope (/int T(\lambda)/\lambda d\lambda)
 * \param exp_time Total exposure time
 * \param exp_num Number of exposures
 * \param back_mag Flux due to background and/or sky in mag/arcsec^2
 * \param ron Read-out noise in electrons/pixel
 * \param psf_file Input PSF image
 * \param oversample Oversampling rate of the PSF image
 */
Observation::Observation(float diameter, float transmission, float exp_time, int exp_num, float back_mag, float ron, std::string psf_file, float oversample):
		diameter(diameter), transmission(transmission), exp_time(exp_time), exp_num(exp_num), back_mag(back_mag), ron(ron), oversample(oversample)
		{
	mag_zeropoint = 2.5*log10(diameter*diameter*transmission*pi/4./hplanck) - 48.6;

#ifdef ENABLE_FITS

	std::auto_ptr<CCfits::FITS> fp (new CCfits::FITS (psf_file.c_str(), CCfits::Read));
	CCfits::PHDU *h0=&fp->pHDU();
	int side_psf = h0->axis(0);
	int N_psf = side_psf*side_psf;
	map_psf.resize(N_psf);
	h0->read(map_psf);

#else
		std::cout << "Please enable the preprocessor flag ENABLE_FITS !" << std::endl;
		exit(1);
#endif
		telescope = false;
		}

/// Reads in and sets the PSF from a fits file. If the pixel size of the fits is different (smaller) than the one of the telescope, it must be specified.
void Observation::setPSF(std::string psf_file, float os)
{
#ifdef ENABLE_FITS
    std::auto_ptr<CCfits::FITS> fp (new CCfits::FITS (psf_file.c_str(), CCfits::Read));
	CCfits::PHDU *h0=&fp->pHDU();
	int side_psf = h0->axis(0);
	int N_psf = side_psf*side_psf;
	map_psf.resize(N_psf);
	h0->read(map_psf);
    
#else
    std::cout << "Please enable the preprocessor flag ENABLE_FITS !" << std::endl;
    exit(1);
#endif
    
    oversample = os;

}

/**  \brief Converts the input map to a realistic image
 *
 * \param map Input map in photons/(cm^2*Hz)
 * \param psf Decides if the psf smoothing is applied
 * \param noise Decides if noise is added
 * \param unit Decides units of output (if flux, output is in 10**(-0.4*mag)) 
 */
PixelMap Observation::Convert (PixelMap &map, bool psf, bool noise, long *seed, unitType unit)
{
	if (telescope == true && fabs(map.getResolution()-pix_size) > std::numeric_limits<double>::epsilon())
	{
		std::cout << "The resolution of the input map is different from the one of the simulated instrument in Observation::Convert!" << std::endl;
		throw std::runtime_error("The resolution of the input map is different from the one of the simulated instrument!");
	}
	PixelMap outmap = PhotonToCounts(map);
	if (psf == true)  outmap = ApplyPSF(outmap);
	if (noise == true) outmap = AddNoise(outmap,seed);
    
    if (unit == flux)
    {
        double counts_to_flux = pow(10,-0.4*mag_zeropoint);
        outmap.Renormalize(counts_to_flux);
    }
	return outmap;
}

/// Converts an observed image to the units of the lensing simulation
PixelMap Observation::Convert_back (PixelMap &map)
{
	PixelMap outmap(map);
	double Q = pow(10,0.4*(mag_zeropoint+48.6))*hplanck;
	outmap.Renormalize(1./Q);
	return outmap;
}


/** * \brief Smooths the image with a PSF map.
*
*/
PixelMap Observation::ApplyPSF(PixelMap &pmap)
{
	if (map_psf.size() == 0)
	{
		if (seeing > 0.)
		{
			PixelMap outmap(pmap);
			outmap.smooth(seeing/2.355);
			return outmap;
		}
		else
		{
			return pmap;
		}
	}
	else
	{
#ifdef ENABLE_FITS
#ifdef ENABLE_FFTW

        PixelMap outmap(pmap);

        // calculates normalisation of psf
        int N_psf = map_psf.size();
        int side_psf = sqrt(N_psf);
        double map_norm = 0.;
        for (int i = 0; i < N_psf; i++)
        {
            map_norm += map_psf[i];
        }
        fftw_plan p_psf;

        // creates plane for fft of map, sets properly input and output data, then performs fft
        fftw_plan p;
        long Npix = outmap.getNpixels();
        long Npix_zeropad = Npix + side_psf;
        std::complex<double>* out=new std::complex<double> [Npix_zeropad*(Npix_zeropad/2+1)];

        // add zero-padding
        double* in_zeropad = new double[Npix_zeropad*Npix_zeropad];
        for (int i = 0; i < Npix_zeropad*Npix_zeropad; i++)
        {
            long ix = i/Npix_zeropad;
            long iy = i%Npix_zeropad;
            if (ix >= side_psf/2 && ix <= (Npix_zeropad-1)-side_psf/2 && iy >= side_psf/2 && iy <= (Npix_zeropad-1)-side_psf/2)
                in_zeropad[i] = outmap[(ix-side_psf/2)*Npix+(iy-side_psf/2)];
            else
                in_zeropad[i] = 0.;
        }
        
        p = fftw_plan_dft_r2c_2d(Npix_zeropad,Npix_zeropad,in_zeropad, reinterpret_cast<fftw_complex*>(out), FFTW_ESTIMATE);
        fftw_execute(p);

        // creates plane for perform backward fft after convolution, sets output data
        fftw_plan p2;
        double* out2 = new double[Npix_zeropad*Npix_zeropad];
        p2 = fftw_plan_dft_c2r_2d(Npix_zeropad,Npix_zeropad,reinterpret_cast<fftw_complex*>(out), out2, FFTW_ESTIMATE);


        // arrange psf data for fft, creates plane, then performs fft
        long psf_big_zeropad_Npixels = static_cast<int>((Npix+side_psf)*oversample);
        double* psf_big_zeropad = new double[psf_big_zeropad_Npixels*psf_big_zeropad_Npixels];
        std::complex<double>* out_psf=new std::complex<double> [psf_big_zeropad_Npixels*(psf_big_zeropad_Npixels/2+1)];
        p_psf = fftw_plan_dft_r2c_2d(psf_big_zeropad_Npixels,psf_big_zeropad_Npixels,psf_big_zeropad, reinterpret_cast<fftw_complex*>(out_psf), FFTW_ESTIMATE);
        long ix, iy;
        for (int i = 0; i < psf_big_zeropad_Npixels*psf_big_zeropad_Npixels; i++)
        {
            ix = i/psf_big_zeropad_Npixels;
            iy = i%psf_big_zeropad_Npixels;
            if(ix<side_psf/2 && iy<side_psf/2)
                psf_big_zeropad[i] = map_psf[(ix+side_psf/2)*side_psf+(iy+side_psf/2)]/map_norm;
            else if(ix<side_psf/2 && iy>=psf_big_zeropad_Npixels-side_psf/2)
                psf_big_zeropad[i] = map_psf[(ix+side_psf/2)*side_psf+(iy-(psf_big_zeropad_Npixels-side_psf/2))]/map_norm;
            else if(ix>=psf_big_zeropad_Npixels-side_psf/2 && iy<side_psf/2)
                psf_big_zeropad[i] = map_psf[(ix-(psf_big_zeropad_Npixels-side_psf/2))*side_psf+(iy+side_psf/2)]/map_norm;
            else if(ix>=psf_big_zeropad_Npixels-side_psf/2 && iy>=psf_big_zeropad_Npixels-side_psf/2)
                psf_big_zeropad[i] = map_psf[(ix-(psf_big_zeropad_Npixels-side_psf/2))*side_psf+(iy-(psf_big_zeropad_Npixels-side_psf/2))]/map_norm;
            else
                psf_big_zeropad[i] = 0.;
        }
        fftw_execute(p_psf);

        // performs convolution in Fourier space, and transforms back to real space
        for (unsigned long i = 0; i < Npix_zeropad*(Npix_zeropad/2+1); i++)
        {
            ix = i/(Npix_zeropad/2+1);
            iy = i%(Npix_zeropad/2+1);
            if (ix>Npix_zeropad/2)
                out[i] *= out_psf[(psf_big_zeropad_Npixels-(Npix_zeropad-ix))*(psf_big_zeropad_Npixels/2+1)+iy];
            else
                out[i] *= out_psf[ix*(psf_big_zeropad_Npixels/2+1)+iy];
        }
        fftw_execute(p2);

        // translates array of data in (normalised) counts map
        for (unsigned long i = 0; i < Npix_zeropad*Npix_zeropad; i++)
        {
            ix = i/Npix_zeropad;
            iy = i%Npix_zeropad;
            if (ix >= side_psf/2 && ix <= (Npix_zeropad-1)-side_psf/2 && iy >= side_psf/2 && iy <= (Npix_zeropad-1)-side_psf/2)
            {
                int ii = (ix-side_psf/2)*Npix+(iy-side_psf/2);
                outmap.AssignValue(ii,out2[i]/double(Npix_zeropad*Npix_zeropad));
            }
        }
        return outmap;
#else
		std::cout << "Please enable the preprocessor flag ENABLE_FFTW !" << std::endl;
		exit(1);
#endif

#else
		std::cout << "Please enable the preprocessor flag ENABLE_FITS !" << std::endl;
		exit(1);
#endif
	}
}

/// Applies realistic noise (read-out + Poisson) on an image
PixelMap Observation::AddNoise(PixelMap &pmap,long *seed)
{
	PixelMap outmap(pmap);
	double Q = pow(10,0.4*(mag_zeropoint+48.6));
	double res_in_arcsec = outmap.getResolution()*180.*60.*60/pi;
	double back_mean = pow(10,-0.4*(48.6+back_mag))*res_in_arcsec*res_in_arcsec*Q*exp_time;
	double rms, noise;
	double norm_map;
	for (unsigned long i = 0; i < outmap.getNpixels()*outmap.getNpixels(); i++)
	{
		norm_map = outmap[i]*exp_time;
		if (norm_map+back_mean > 500.)
		{
			rms = sqrt(pow(exp_num*ron,2)+norm_map+back_mean);
			noise = gasdev(seed)*rms;
			outmap.AssignValue(i,double(norm_map+noise)/exp_time);
		}
		else
		{
			int k = 0;
			double p = 1.;
			double L = exp(-(norm_map+back_mean));
			while	(p > L)
			{
				k++;
				p *= ran2(seed);
			}
			outmap.AssignValue(i,double(k-1-back_mean)/exp_time);
		}
	}
	return outmap;<commit_msg>Little correction to read-out noise term in Observation::AddNoise<commit_after>/*
 * observation.cpp
 *
 *  Created on: Apr 29, 2013
 *      Author: <NAME>
 */

#include "slsimlib.h"

#ifdef ENABLE_FITS
#include <CCfits/CCfits>
//#include <CCfits>
#endif

#ifdef ENABLE_FFTW
#include "fftw3.h"
#endif

#include <fstream>
#include <limits>


/** * \brief Creates an observation setup that mimics a known instrument
 *
 */
Observation::Observation(Telescope tel_name)
{
	if (tel_name == Euclid_VIS)
	{
		diameter = 119.;
		transmission = 0.30;
		exp_time = 1800.;
		exp_num = 3;
		back_mag = 22.8;
		ron = 5.;
		seeing = 0.18;
		pix_size = .1/60./60./180.*pi;
	}
	if (tel_name == Euclid_Y)
	{
		diameter = 119.;
		transmission = 0.0961;
		exp_time = 264.;
		exp_num = 3;
		back_mag = 22.57;
		ron = 5.;
		seeing = 0.3;
		pix_size = .3/60./60./180.*pi;
	}
	if (tel_name == Euclid_J)
	{
		diameter = 119.;
		transmission = 0.0814;
		exp_time = 270.;
		exp_num = 3;
		back_mag = 22.53;
		ron = 5.;
		seeing = 0.3;
		pix_size = .3/60./60./180.*pi;
	}
	if (tel_name == Euclid_H)
	{
		diameter = 119.;
		transmission = 0.1692;
		exp_time = 162.;
		exp_num = 3;
		back_mag = 22.59;
		ron = 5.;
		seeing = 0.3;
		pix_size = .3/60./60./180.*pi;
	}
	if (tel_name == KiDS_u)
	{
		diameter = 265.;
		transmission = 0.032;
		exp_time = 1000.;
		exp_num = 5;
		back_mag = 22.93;
		ron = 5.;
		seeing = 1.0;
		pix_size = .2/60./60./180.*pi;
	}
	if (tel_name == KiDS_g)
	{
		diameter = 265.;
		transmission = 0.1220;
		exp_time = 900.;
		exp_num = 5;
		back_mag = 22.29;
		ron = 5.;
		seeing = 0.8;
		pix_size = .2/60./60./180.*pi;
	}
	if (tel_name == KiDS_r)
	{
		diameter = 265.;
		transmission = 0.089;
		exp_time = 1800.;
		exp_num = 5;
		back_mag = 21.40;
		ron = 5.;
		seeing = 0.7;
		pix_size = .2/60./60./180.*pi;
	}
	if (tel_name == KiDS_i)
	{
		diameter = 265.;
		transmission = 0.062;
		exp_time = 1200.;
		exp_num = 5;
		back_mag = 20.64;
		ron = 5.;
		seeing = 1.1;
		pix_size = .2/60./60./180.*pi;
	}

	mag_zeropoint = 2.5*log10(diameter*diameter*transmission*pi/4./hplanck) - 48.6;
	telescope = true;
}

/** *  Creates a custom observation setup with parameters decided by the user.
 *
 * \param diameter Diameter of telescope (in cm) (Collecting area is pi/4.*diameter^2)
 * \param transmission Total transmission of the telescope (/int T(\lambda)/\lambda d\lambda)
 * \param exp_time Total exposure time
 * \param exp_num Number of exposures
 * \param back_mag Flux due to background and/or sky in mag/arcsec^2
 * \param ron Read-out noise in electrons/pixel
 * \param seeing FWHM in arcsecs of the image
 */
Observation::Observation(float diameter, float transmission, float exp_time, int exp_num, float back_mag, float ron, float seeing):
		diameter(diameter), transmission(transmission), exp_time(exp_time), exp_num(exp_num), back_mag(back_mag), ron(ron), seeing(seeing)
		{
			mag_zeropoint = 2.5*log10(diameter*diameter*transmission*pi/4./hplanck) - 48.6;
			telescope = false;
		}

/**  Creates a custom observation setup with parameters decided by the user. Allows for the use of a psf fits image.
 *
 * \param diameter Diameter of telescope (in cm) (Collecting area is pi/4.*diameter^2)
 * \param transmission Total transmission of the telescope (/int T(\lambda)/\lambda d\lambda)
 * \param exp_time Total exposure time
 * \param exp_num Number of exposures
 * \param back_mag Flux due to background and/or sky in mag/arcsec^2
 * \param ron Read-out noise in electrons/pixel
 * \param psf_file Input PSF image
 * \param oversample Oversampling rate of the PSF image
 */
Observation::Observation(float diameter, float transmission, float exp_time, int exp_num, float back_mag, float ron, std::string psf_file, float oversample):
		diameter(diameter), transmission(transmission), exp_time(exp_time), exp_num(exp_num), back_mag(back_mag), ron(ron), oversample(oversample)
		{
	mag_zeropoint = 2.5*log10(diameter*diameter*transmission*pi/4./hplanck) - 48.6;

#ifdef ENABLE_FITS

	std::auto_ptr<CCfits::FITS> fp (new CCfits::FITS (psf_file.c_str(), CCfits::Read));
	CCfits::PHDU *h0=&fp->pHDU();
	int side_psf = h0->axis(0);
	int N_psf = side_psf*side_psf;
	map_psf.resize(N_psf);
	h0->read(map_psf);

#else
		std::cout << "Please enable the preprocessor flag ENABLE_FITS !" << std::endl;
		exit(1);
#endif
		telescope = false;
		}

/// Reads in and sets the PSF from a fits file. If the pixel size of the fits is different (smaller) than the one of the telescope, it must be specified.
void Observation::setPSF(std::string psf_file, float os)
{
#ifdef ENABLE_FITS
    std::auto_ptr<CCfits::FITS> fp (new CCfits::FITS (psf_file.c_str(), CCfits::Read));
	CCfits::PHDU *h0=&fp->pHDU();
	int side_psf = h0->axis(0);
	int N_psf = side_psf*side_psf;
	map_psf.resize(N_psf);
	h0->read(map_psf);
    
#else
    std::cout << "Please enable the preprocessor flag ENABLE_FITS !" << std::endl;
    exit(1);
#endif
    
    oversample = os;

}

/**  \brief Converts the input map to a realistic image
 *
 * \param map Input map in photons/(cm^2*Hz)
 * \param psf Decides if the psf smoothing is applied
 * \param noise Decides if noise is added
 * \param unit Decides units of output (if flux, output is in 10**(-0.4*mag)) 
 */
PixelMap Observation::Convert (PixelMap &map, bool psf, bool noise, long *seed, unitType unit)
{
	if (telescope == true && fabs(map.getResolution()-pix_size) > std::numeric_limits<double>::epsilon())
	{
		std::cout << "The resolution of the input map is different from the one of the simulated instrument in Observation::Convert!" << std::endl;
		throw std::runtime_error("The resolution of the input map is different from the one of the simulated instrument!");
	}
	PixelMap outmap = PhotonToCounts(map);
	if (psf == true)  outmap = ApplyPSF(outmap);
	if (noise == true) outmap = AddNoise(outmap,seed);
    
    if (unit == flux)
    {
        double counts_to_flux = pow(10,-0.4*mag_zeropoint);
        outmap.Renormalize(counts_to_flux);
    }
	return outmap;
}

/// Converts an observed image to the units of the lensing simulation
PixelMap Observation::Convert_back (PixelMap &map)
{
	PixelMap outmap(map);
	double Q = pow(10,0.4*(mag_zeropoint+48.6))*hplanck;
	outmap.Renormalize(1./Q);
	return outmap;
}


/** * \brief Smooths the image with a PSF map.
*
*/
PixelMap Observation::ApplyPSF(PixelMap &pmap)
{
	if (map_psf.size() == 0)
	{
		if (seeing > 0.)
		{
			PixelMap outmap(pmap);
			outmap.smooth(seeing/2.355);
			return outmap;
		}
		else
		{
			return pmap;
		}
	}
	else
	{
#ifdef ENABLE_FITS
#ifdef ENABLE_FFTW

        PixelMap outmap(pmap);

        // calculates normalisation of psf
        int N_psf = map_psf.size();
        int side_psf = sqrt(N_psf);
        double map_norm = 0.;
        for (int i = 0; i < N_psf; i++)
        {
            map_norm += map_psf[i];
        }
        fftw_plan p_psf;

        // creates plane for fft of map, sets properly input and output data, then performs fft
        fftw_plan p;
        long Npix = outmap.getNpixels();
        long Npix_zeropad = Npix + side_psf;
        std::complex<double>* out=new std::complex<double> [Npix_zeropad*(Npix_zeropad/2+1)];

        // add zero-padding
        double* in_zeropad = new double[Npix_zeropad*Npix_zeropad];
        for (int i = 0; i < Npix_zeropad*Npix_zeropad; i++)
        {
            long ix = i/Npix_zeropad;
            long iy = i%Npix_zeropad;
            if (ix >= side_psf/2 && ix <= (Npix_zeropad-1)-side_psf/2 && iy >= side_psf/2 && iy <= (Npix_zeropad-1)-side_psf/2)
                in_zeropad[i] = outmap[(ix-side_psf/2)*Npix+(iy-side_psf/2)];
            else
                in_zeropad[i] = 0.;
        }
        
        p = fftw_plan_dft_r2c_2d(Npix_zeropad,Npix_zeropad,in_zeropad, reinterpret_cast<fftw_complex*>(out), FFTW_ESTIMATE);
        fftw_execute(p);

        // creates plane for perform backward fft after convolution, sets output data
        fftw_plan p2;
        double* out2 = new double[Npix_zeropad*Npix_zeropad];
        p2 = fftw_plan_dft_c2r_2d(Npix_zeropad,Npix_zeropad,reinterpret_cast<fftw_complex*>(out), out2, FFTW_ESTIMATE);


        // arrange psf data for fft, creates plane, then performs fft
        long psf_big_zeropad_Npixels = static_cast<int>((Npix+side_psf)*oversample);
        double* psf_big_zeropad = new double[psf_big_zeropad_Npixels*psf_big_zeropad_Npixels];
        std::complex<double>* out_psf=new std::complex<double> [psf_big_zeropad_Npixels*(psf_big_zeropad_Npixels/2+1)];
        p_psf = fftw_plan_dft_r2c_2d(psf_big_zeropad_Npixels,psf_big_zeropad_Npixels,psf_big_zeropad, reinterpret_cast<fftw_complex*>(out_psf), FFTW_ESTIMATE);
        long ix, iy;
        for (int i = 0; i < psf_big_zeropad_Npixels*psf_big_zeropad_Npixels; i++)
        {
            ix = i/psf_big_zeropad_Npixels;
            iy = i%psf_big_zeropad_Npixels;
            if(ix<side_psf/2 && iy<side_psf/2)
                psf_big_zeropad[i] = map_psf[(ix+side_psf/2)*side_psf+(iy+side_psf/2)]/map_norm;
            else if(ix<side_psf/2 && iy>=psf_big_zeropad_Npixels-side_psf/2)
                psf_big_zeropad[i] = map_psf[(ix+side_psf/2)*side_psf+(iy-(psf_big_zeropad_Npixels-side_psf/2))]/map_norm;
            else if(ix>=psf_big_zeropad_Npixels-side_psf/2 && iy<side_psf/2)
                psf_big_zeropad[i] = map_psf[(ix-(psf_big_zeropad_Npixels-side_psf/2))*side_psf+(iy+side_psf/2)]/map_norm;
            else if(ix>=psf_big_zeropad_Npixels-side_psf/2 && iy>=psf_big_zeropad_Npixels-side_psf/2)
                psf_big_zeropad[i] = map_psf[(ix-(psf_big_zeropad_Npixels-side_psf/2))*side_psf+(iy-(psf_big_zeropad_Npixels-side_psf/2))]/map_norm;
            else
                psf_big_zeropad[i] = 0.;
        }
        fftw_execute(p_psf);

        // performs convolution in Fourier space, and transforms back to real space
        for (unsigned long i = 0; i < Npix_zeropad*(Npix_zeropad/2+1); i++)
        {
            ix = i/(Npix_zeropad/2+1);
            iy = i%(Npix_zeropad/2+1);
            if (ix>Npix_zeropad/2)
                out[i] *= out_psf[(psf_big_zeropad_Npixels-(Npix_zeropad-ix))*(psf_big_zeropad_Npixels/2+1)+iy];
            else
                out[i] *= out_psf[ix*(psf_big_zeropad_Npixels/2+1)+iy];
        }
        fftw_execute(p2);

        // translates array of data in (normalised) counts map
        for (unsigned long i = 0; i < Npix_zeropad*Npix_zeropad; i++)
        {
            ix = i/Npix_zeropad;
            iy = i%Npix_zeropad;
            if (ix >= side_psf/2 && ix <= (Npix_zeropad-1)-side_psf/2 && iy >= side_psf/2 && iy <= (Npix_zeropad-1)-side_psf/2)
            {
                int ii = (ix-side_psf/2)*Npix+(iy-side_psf/2);
                outmap.AssignValue(ii,out2[i]/double(Npix_zeropad*Npix_zeropad));
            }
        }
        return outmap;
#else
		std::cout << "Please enable the preprocessor flag ENABLE_FFTW !" << std::endl;
		exit(1);
#endif

#else
		std::cout << "Please enable the preprocessor flag ENABLE_FITS !" << std::endl;
		exit(1);
#endif
	}
}

/// Applies realistic noise (read-out + Poisson) on an image
PixelMap Observation::AddNoise(PixelMap &pmap,long *seed)
{
	PixelMap outmap(pmap);
	double Q = pow(10,0.4*(mag_zeropoint+48.6));
	double res_in_arcsec = outmap.getResolution()*180.*60.*60/pi;
	double back_mean = pow(10,-0.4*(48.6+back_mag))*res_in_arcsec*res_in_arcsec*Q*exp_time;
	double rms, noise;
	double norm_map;
	for (unsigned long i = 0; i < outmap.getNpixels()*outmap.getNpixels(); i++)
	{
		norm_map = outmap[i]*exp_time;
		if (norm_map+back_mean > 500.)
		{
			rms = sqrt(exp_num*ron*ron+norm_map+back_mean);
			noise = gasdev(seed)*rms;
			outmap.AssignValue(i,double(norm_map+noise)/exp_time);
		}
		else
		{
			int k = 0;
			double p = 1.;
			double L = exp(-(norm_map+back_mean));
			while	(p > L)
			{
				k++;
				p *= ran2(seed);
			}
			outmap.AssignValue(i,double(k-1-back_mean)/exp_time);
		}
	}
	return outmap;
{'fraction_non_alphanumeric': 0.1360518453057274, 'fraction_numerical': 0.032277069393595086, 'mean_word_length': 3.1422346173958515, 'pattern_counts': {'":': 0, '<': 103, '<?xml version=': 0, '>': 73, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1661412', 'n_tokens_mistral': 10324, 'n_tokens_neox': 9379, 'n_words': 1984}