/*
 
 convolution_kernel.cc - This file is part of MUSIC -
 a code to generate multi-scale initial conditions 
 for cosmological simulations 
 
 Copyright (C) 2010  Oliver Hahn
 
*/

#include "general.hh"
#include "densities.hh"
#include "convolution_kernel.hh"

#if defined(FFTW3) && defined( SINGLE_PRECISION)
//#define fftw_complex fftwf_complex
typedef fftw_complex fftwf_complex;
#endif

double T0 = 1.0;

namespace convolution{

	std::map< std::string, kernel_creator *>& 
	get_kernel_map()
	{
		static std::map< std::string, kernel_creator* > kernel_map;
		return kernel_map;
	}
	
	template< typename real_t >
	void perform( kernel * pk, void *pd, bool shift )
	{
		//return;
		
		parameters cparam_ = pk->cparam_;
		double fftnorm = pow(2.0*M_PI,1.5)/sqrt(cparam_.lx*cparam_.ly*cparam_.lz)/sqrt((double)cparam_.nx*(double)cparam_.ny*(double)cparam_.nz);
		
		fftw_complex	*cdata,*ckernel;
		fftw_real		*data;
		
		data		= reinterpret_cast<fftw_real*>(pd);
		cdata		= reinterpret_cast<fftw_complex*>(data);
		ckernel		= reinterpret_cast<fftw_complex*>( pk->get_ptr() );
		
		
		
		std::cout << "   - Performing density convolution... (" 
			<< cparam_.nx <<  ", " << cparam_.ny << ", " << cparam_.nz << ")\n";
		
		
		LOGUSER("Performing kernel convolution on (%5d,%5d,%5d) grid",cparam_.nx ,cparam_.ny ,cparam_.nz );
		LOGUSER("Performing forward FFT...");
#ifdef FFTW3
	#ifdef SINGLE_PRECISION
		fftwf_plan plan, iplan;
		plan = fftwf_plan_dft_r2c_3d( cparam_.nx, cparam_.ny, cparam_.nz, data, cdata, FFTW_ESTIMATE);
		iplan = fftwf_plan_dft_c2r_3d(cparam_.nx, cparam_.ny, cparam_.nz, cdata, data, FFTW_ESTIMATE);
		
		fftwf_execute(plan);
	#else
		fftw_plan plan, iplan;
		plan = fftw_plan_dft_r2c_3d( cparam_.nx, cparam_.ny, cparam_.nz, data, cdata, FFTW_ESTIMATE);
		iplan = fftw_plan_dft_c2r_3d(cparam_.nx, cparam_.ny, cparam_.nz, cdata, data, FFTW_ESTIMATE);
		
		fftw_execute(plan);
	#endif
#else
		rfftwnd_plan	iplan, plan;
		
		plan  = rfftw3d_create_plan( cparam_.nx, cparam_.ny, cparam_.nz,
									FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE|FFTW_IN_PLACE);
		
		iplan = rfftw3d_create_plan( cparam_.nx, cparam_.ny, cparam_.nz,
									FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE|FFTW_IN_PLACE);
		
#ifndef SINGLETHREAD_FFTW		
		rfftwnd_threads_one_real_to_complex( omp_get_max_threads(), plan, data, NULL );
#else
		rfftwnd_one_real_to_complex( plan, data, NULL );
#endif
		
#endif
		//..... need a phase shift for baryons for SPH
		double dstag = 0.0;
        
		if( shift )
		{	
			double boxlength = pk->pcf_->getValue<double>("setup","boxlength");
            double stagfact  = pk->pcf_->getValueSafe<double>("setup","baryon_staggering",0.5);
			int lmax = pk->pcf_->getValue<int>("setup","levelmax");	
			double dxmax = boxlength/(1<<lmax);
			double dxcur = cparam_.lx/cparam_.nx;
			//std::cerr << "Performing staggering shift for SPH\n";
			LOGUSER("Performing staggering shift for SPH");
			dstag = stagfact * 2.0 * M_PI/cparam_.nx * dxmax/dxcur;
		}
		
		//.............................................
		
		#pragma omp parallel for
		for( int i=0; i<cparam_.nx; ++i )
			for( int j=0; j<cparam_.ny; ++j )
				for( int k=0; k<cparam_.nz/2+1; ++k )
				{
					size_t ii = (size_t)(i*cparam_.ny + j) * (size_t)(cparam_.nz/2+1) + (size_t)k;
					
					double kx,ky,kz;
					
					kx = (double)i;
					ky = (double)j;
					kz = (double)k;
					
					if( kx > cparam_.nx/2 ) kx -= cparam_.nx;
					if( ky > cparam_.ny/2 ) ky -= cparam_.ny;
					
					double arg = (kx+ky+kz)*dstag;
					std::complex<double> carg( cos(arg), sin(arg) );
					
					std::complex<double> 
						ccdata(RE(cdata[ii]),IM(cdata[ii])), 
						cckernel(RE(ckernel[ii]),IM(ckernel[ii]));
					
					ccdata = ccdata * cckernel *fftnorm * carg;
					
					RE(cdata[ii]) = ccdata.real();
					IM(cdata[ii]) = ccdata.imag();
				}

		LOGUSER("Performing backward FFT...");
		
#ifdef FFTW3
	#ifdef SINGLE_PRECISION
		fftwf_execute(iplan);
		fftwf_destroy_plan(plan);
		fftwf_destroy_plan(iplan);
	#else	
		fftw_execute(iplan);
		fftw_destroy_plan(plan);
		fftw_destroy_plan(iplan);
		
	#endif
#else
		#ifndef SINGLETHREAD_FFTW		
		rfftwnd_threads_one_complex_to_real( omp_get_max_threads(), iplan, cdata, NULL);
		#else		
		rfftwnd_one_complex_to_real(iplan, cdata, NULL);
		#endif
		
		rfftwnd_destroy_plan(plan);
		rfftwnd_destroy_plan(iplan);
#endif
	}
	
	
	template void perform<double>( kernel* pk, void *pd, bool shift );
	template void perform<float>( kernel* pk, void *pd, bool shift );
	
	/*****************************************************************************************\
	 ***    SPECIFIC KERNEL IMPLEMENTATIONS      *********************************************
	\*****************************************************************************************/	
	

	template< typename real_t >
	class kernel_k : public kernel
	{
	protected:
		std::vector<real_t> kdata_;
		
		void compute_kernel( tf_type type );
		
	public:
		kernel_k( config_file& cf, transfer_function* ptf, refinement_hierarchy& refh, tf_type type )
		: kernel( cf, ptf, refh, type )
		{	}
		
		kernel* fetch_kernel( int ilevel, bool isolated=false );
		
		void *get_ptr()
		{	return reinterpret_cast<void*> (&kdata_[0]);	}
		
		~kernel_k()
		{	deallocate();	}
		
		void deallocate()
		{	std::vector<real_t>().swap( kdata_ );	}
		
	};
	
	template< typename real_t >
	kernel* kernel_k<real_t>::fetch_kernel( int ilevel, bool isolated )
	{
		double 
			boxlength	= pcf_->getValue<double>("setup","boxlength"),
			nspec		= pcf_->getValue<double>("cosmology","nspec"),
			pnorm		= pcf_->getValue<double>("cosmology","pnorm"),
			fac			= pow(boxlength,3)/pow(2.0*M_PI,3);
		
		TransferFunction_k *tfk = new TransferFunction_k(type_,ptf_,nspec,pnorm);
		
		int 
			nx = prefh_->size(prefh_->levelmax(),0),
			ny = prefh_->size(prefh_->levelmax(),1),
			nz = prefh_->size(prefh_->levelmax(),2);
		
		cparam_.nx = nx;
		cparam_.ny = ny;
		cparam_.nz = nz;
		cparam_.lx = boxlength;
		cparam_.ly = boxlength;
		cparam_.lz = boxlength;
		cparam_.pcf = pcf_;
		
		kdata_.assign( (size_t)nx*(size_t)ny*(size_t)(nz+2), 0.0 );
		
		fftw_complex *kdata = reinterpret_cast<fftw_complex*> ( this->get_ptr() );
		
		unsigned nzp = (nz/2+1);
		fac =1.0;
		
		double kfac = 2.0*M_PI/boxlength, ksum = 0.0;
		size_t kcount = 0;
		
		#pragma omp parallel for reduction(+:ksum,kcount)
		for( int i=0; i<nx; ++i )
			for( int j=0; j<ny; ++j )
				for( int k=0; k<nz/2+1; ++k )
				{
					double kx,ky,kz;
					
					kx = (double)i;
					ky = (double)j;
					kz = (double)k;
					
					if( kx > nx/2 ) kx -= nx;
					if( ky > ny/2 ) ky -= ny;
					
					size_t q = ((size_t)i*ny+(size_t)j)*nzp+(size_t)k;
					
					RE(kdata[q]) = fac*tfk->compute(kfac*sqrt(kx*kx+ky*ky+kz*kz));
					IM(kdata[q]) = 0.0;
					
					if( k==0 || k==nz/2 )
					{
						ksum  += RE(kdata[q]);
						kcount++;
					}else{
						ksum  += 2.0*(RE(kdata[q]));
						kcount+=2;
					}
				}
		
		delete tfk;
		
		return this;
	}
	
	
	////////////////////////////////////////////////////////////////////////////
	
	template< typename real_t >
	class kernel_real_cached : public kernel
	{
	protected:
		std::vector<real_t> kdata_;
		void precompute_kernel( transfer_function* ptf, tf_type type, const refinement_hierarchy& refh );
		
	public:
		kernel_real_cached( config_file& cf, transfer_function* ptf, refinement_hierarchy& refh, tf_type type )
		: kernel( cf, ptf, refh, type )
		{
			precompute_kernel(ptf, type, refh);
		}
		
		kernel* fetch_kernel( int ilevel, bool isolated=false );
		
		void *get_ptr()
		{	return reinterpret_cast<void*> (&kdata_[0]);	}
		
		~kernel_real_cached()
		{	
			deallocate();
		}
		
		void deallocate()
		{
			std::vector<real_t>().swap( kdata_ );	
		}
		
	};
	
	template< typename real_t >
	kernel* kernel_real_cached<real_t>::fetch_kernel( int ilevel, bool isolated )
	{
		char cachefname[128];
		sprintf(cachefname,"temp_kernel_level%03d.tmp",ilevel);
		FILE *fp = fopen(cachefname,"r");
		
		
		std::cout << " - Fetching kernel for level " << ilevel << std::endl;
		
		LOGUSER("Loading kernel for level %3d from file \'%s\'...",ilevel,cachefname);
		
		if( fp == NULL )
		{	
			LOGERR("Could not open kernel file \'%s\'.",cachefname);
			throw std::runtime_error("Internal error: cached convolution kernel does not exist on disk!");
		}
		
		unsigned nx,ny,nz;
		size_t nread = 0;
		
		nread = fread(	reinterpret_cast<void*> (&nx), sizeof(unsigned), 1, fp );
		nread = fread(	reinterpret_cast<void*> (&ny), sizeof(unsigned), 1, fp );
		nread = fread(	reinterpret_cast<void*> (&nz), sizeof(unsigned), 1, fp );
		
		kdata_.assign((size_t)nx*(size_t)ny*(size_t)nz,0.0);
		
		for( size_t ix=0;ix<nx; ++ix )
		{	
			const size_t sz = ny*nz;
			nread = fread( reinterpret_cast<void*>(&kdata_[(size_t)ix * sz]), sizeof(fftw_real), sz, fp );
			assert( nread == sz );
		}

		fclose(fp);
		
		//... set parameters
		
		double boxlength = pcf_->getValue<double>("setup","boxlength");
		double dx = boxlength / (1<<ilevel);
		
		cparam_.nx = nx;
		cparam_.ny = ny;
		cparam_.nz = nz-2;
		cparam_.lx = dx * cparam_.nx;
		cparam_.ly = dx * cparam_.ny;
		cparam_.lz = dx * cparam_.nz;
		cparam_.pcf = pcf_;
		
		fftw_real		*rkernel	= reinterpret_cast<fftw_real*>( &kdata_[0] );
		
#ifdef FFTW3
	#ifdef SINGLE_PRECISION
		fftwf_complex		*kkernel	= reinterpret_cast<fftwf_complex*> (&rkernel[0]);
		fftwf_plan plan = fftwf_plan_dft_r2c_3d( cparam_.nx, cparam_.ny, cparam_.nz, rkernel, kkernel, FFTW_ESTIMATE);
		fftwf_execute(plan);
		fftwf_destroy_plan(plan);
	#else	
		fftw_complex		*kkernel	= reinterpret_cast<fftw_complex*> (&rkernel[0]);
		fftw_plan plan = fftw_plan_dft_r2c_3d( cparam_.nx, cparam_.ny, cparam_.nz, rkernel, kkernel, FFTW_ESTIMATE);
		fftw_execute(plan);
		fftw_destroy_plan(plan);
	#endif
#else
		rfftwnd_plan	plan		= rfftw3d_create_plan( cparam_.nx, cparam_.ny, cparam_.nz,
														  FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE|FFTW_IN_PLACE);
		
		
	#ifndef SINGLETHREAD_FFTW		
		rfftwnd_threads_one_real_to_complex( omp_get_max_threads(), plan, rkernel, NULL );
	#else
		rfftwnd_one_real_to_complex( plan, rkernel, NULL );
	#endif
		
		rfftwnd_destroy_plan( plan );
#endif
		return this;
	}
	
	
	template< typename real_t >
	void kernel_real_cached<real_t>::precompute_kernel( transfer_function* ptf, tf_type type, const refinement_hierarchy& refh )
	{
		//... compute kernel for finest level
		int nx,ny,nz;
		double dx,lx,ly,lz;
		
		double
			boxlength	= pcf_->getValue<double>("setup","boxlength"),
			boxlength2  = 0.5*boxlength;
		
		int
			levelmax	= refh.levelmax(),
			levelmin	= refh.levelmin();
		
		LOGUSER("Precomputing transfer function kernels...");
			
		nx = refh.size(refh.levelmax(),0);
		ny = refh.size(refh.levelmax(),1);
		nz = refh.size(refh.levelmax(),2);
		
		if( levelmax != levelmin )
		{
			nx *= 2;
			ny *= 2;
			nz *= 2;
		}
		
		dx = boxlength / (1<<refh.levelmax());
		lx = dx * nx;
		ly = dx * ny;
		lz = dx * nz;
		
		double 
		kny			= M_PI/dx,
		fac			= lx*ly*lz/pow(2.0*M_PI,3)/((double)nx*(double)ny*(double)nz),
		
		nspec		= pcf_->getValue<double>("cosmology","nspec"),
		pnorm		= pcf_->getValue<double>("cosmology","pnorm");
		
		bool
			bperiodic	= pcf_->getValueSafe<bool>("setup","periodic_TF",true),
			deconv = pcf_->getValueSafe<bool>("setup","deconvolve",true);
		//		bool deconv_baryons = true;//pcf_->getValueSafe<bool>("setup","deconvolve_baryons",false) || do_SPH;
		bool bsmooth_baryons = false;//type==baryon && !deconv_baryons;
									 //bool bbaryons = type==baryon | type==vbaryon;
		bool kspacepoisson = ((pcf_->getValueSafe<bool>("poisson","fft_fine",true)|
							   pcf_->getValueSafe<bool>("poisson","kspace",false)));// & !(type==baryon&!do_SPH);//&!baryons ;
		
		std::cout << "   - Computing transfer function kernel...\n";
		
		TransferFunction_real *tfr = 
						new TransferFunction_real( boxlength, 1<<levelmax, type, ptf,nspec,pnorm,
								  0.25*dx,2.0*boxlength,kny, (int)pow(2,levelmax+2));		
		
		fftw_real *rkernel = new fftw_real[(size_t)nx*(size_t)ny*((size_t)nz+2)], *rkernel_coarse;
		
		#pragma omp parallel for
		for( int i=0; i<nx; ++i )
			for( int j=0; j<ny; ++j )
				for( int k=0; k<nz; ++k )
				{
					size_t q=((size_t)(i)*ny + (size_t)(j)) * (size_t)(nz+2) + (size_t)(k);
					rkernel[q] = 0.0;
					
				}
		
		LOGUSER("Computing fine kernel (level %d)...", levelmax);
		
		int ref_fac = (deconv&&kspacepoisson)? 2 : 0;
		const int ql = -ref_fac/2+1, qr = ql+ref_fac;
		const double rf8 = pow(ref_fac,3);
		const double dx05 = 0.5*dx, dx025 = 0.25*dx;
		
		if( bperiodic  )
		{		
			
			
			#pragma omp parallel for
			for( int i=0; i<=nx/2; ++i )
				for( int j=0; j<=ny/2; ++j )
					for( int k=0; k<=nz/2; ++k )
					{
						int iix(i), iiy(j), iiz(k);
						double rr[3], rr2;
						
						if( iix > (int)nx/2 ) iix -= nx;
						if( iiy > (int)ny/2 ) iiy -= ny;
						if( iiz > (int)nz/2 ) iiz -= nz;
						
						
						//... speed up 8x by copying data to other octants
						size_t idx[8];
						
						idx[0] = ((size_t)(i)*ny + (size_t)(j)) * 2*(nz/2+1) + (size_t)(k);
						idx[1] = ((size_t)(nx-i)*ny + (size_t)(j)) * 2*(nz/2+1) + (size_t)(k);
						idx[2] = ((size_t)(i)*ny + (size_t)(ny-j)) * 2*(nz/2+1) + (size_t)(k);
						idx[3] = ((size_t)(nx-i)*ny + (size_t)(ny-j)) * 2*(nz/2+1) + (size_t)(k);
						idx[4] = ((size_t)(i)*ny + (size_t)(j)) * 2*(nz/2+1) + (size_t)(nz-k);
						idx[5] = ((size_t)(nx-i)*ny + (size_t)(j)) * 2*(nz/2+1) + (size_t)(nz-k);
						idx[6] = ((size_t)(i)*ny + (size_t)(ny-j)) * 2*(nz/2+1) + (size_t)(nz-k);
						idx[7] = ((size_t)(nx-i)*ny + (size_t)(ny-j)) * 2*(nz/2+1) + (size_t)(nz-k);
						
						if(i==0||i==nx/2){ idx[1]=idx[3]=idx[5]=idx[7]=(size_t)-1;}
						if(j==0||j==ny/2){ idx[2]=idx[3]=idx[6]=idx[7]=(size_t)-1;}
						if(k==0||k==nz/2){ idx[4]=idx[5]=idx[6]=idx[7]=(size_t)-1;}
						
						double val = 0.0;
						
						for( int ii=-1; ii<=1; ++ii )
							for( int jj=-1; jj<=1; ++jj )
								for( int kk=-1; kk<=1; ++kk )
								{
									rr[0] = ((double)iix ) * dx + ii*boxlength;
									rr[1] = ((double)iiy ) * dx + jj*boxlength;
									rr[2] = ((double)iiz ) * dx + kk*boxlength;
									
									if( rr[0] > -boxlength && rr[0] <= boxlength
									   && rr[1] > -boxlength && rr[1] <= boxlength
									   && rr[2] > -boxlength && rr[2] <= boxlength )
									{
										if( ref_fac > 0 )
										{
											double rrr[3];
											register double rrr2[3];
											for( int iii=ql; iii<qr; ++iii )
											{       
												rrr[0] = rr[0]+(double)iii*dx05 - dx025;
												rrr2[0]= rrr[0]*rrr[0];
												for( int jjj=ql; jjj<qr; ++jjj )
												{
													rrr[1] = rr[1]+(double)jjj*dx05 - dx025;
													rrr2[1]= rrr[1]*rrr[1];
													for( int kkk=ql; kkk<qr; ++kkk )
													{
														rrr[2] = rr[2]+(double)kkk*dx05 - dx025;
														rrr2[2]= rrr[2]*rrr[2];
														rr2 = rrr2[0]+rrr2[1]+rrr2[2];
														val += tfr->compute_real(rr2)/rf8;
													}
												}
											}
										}
										else
										{
											rr2 = rr[0]*rr[0]+rr[1]*rr[1]+rr[2]*rr[2];
											val += tfr->compute_real(rr2);	
										}

										
									}
								}
						
						val *= fac;
						
						for(int q=0;q<8;++q)
							if(idx[q]!=(size_t)-1)  
								rkernel[idx[q]] = val;	
					}
			
		}else{
			#pragma omp parallel for
			for( int i=0; i<nx; ++i )
				for( int j=0; j<ny; ++j )
					for( int k=0; k<nz; ++k )
					{
						int iix(i), iiy(j), iiz(k);
						double rr[3], rr2;
						
						if( iix > (int)nx/2 ) iix -= nx;
						if( iiy > (int)ny/2 ) iiy -= ny;
						if( iiz > (int)nz/2 ) iiz -= nz;
						
						//size_t idx = ((size_t)i*ny + (size_t)j) * 2*(nz/2+1) + (size_t)k;
						
						rr[0] = ((double)iix ) * dx;
						rr[1] = ((double)iiy ) * dx;
						rr[2] = ((double)iiz ) * dx;
						
						//rkernel[idx] = 0.0;
						
						//rr2 = rr[0]*rr[0]+rr[1]*rr[1]+rr[2]*rr[2];
						
						//... speed up 8x by copying data to other octants
						size_t idx[8];
						
						idx[0] = ((size_t)(i)*ny + (size_t)(j)) * 2*(nz/2+1) + (size_t)(k);
						idx[1] = ((size_t)(nx-i)*ny + (size_t)(j)) * 2*(nz/2+1) + (size_t)(k);
						idx[2] = ((size_t)(i)*ny + (size_t)(ny-j)) * 2*(nz/2+1) + (size_t)(k);
						idx[3] = ((size_t)(nx-i)*ny + (size_t)(ny-j)) * 2*(nz/2+1) + (size_t)(k);
						idx[4] = ((size_t)(i)*ny + (size_t)(j)) * 2*(nz/2+1) + (size_t)(nz-k);
						idx[5] = ((size_t)(nx-i)*ny + (size_t)(j)) * 2*(nz/2+1) + (size_t)(nz-k);
						idx[6] = ((size_t)(i)*ny + (size_t)(ny-j)) * 2*(nz/2+1) + (size_t)(nz-k);
						idx[7] = ((size_t)(nx-i)*ny + (size_t)(ny-j)) * 2*(nz/2+1) + (size_t)(nz-k);
						
						if(i==0||i==nx/2){ idx[1]=idx[3]=idx[5]=idx[7]=(size_t)-1;}
						if(j==0||j==ny/2){ idx[2]=idx[3]=idx[6]=idx[7]=(size_t)-1;}
						if(k==0||k==nz/2){ idx[4]=idx[5]=idx[6]=idx[7]=(size_t)-1;}
						
						double val = 0.0;//(fftw_real)tfr->compute_real(rr2)*fac;
						
						if( ref_fac > 0 )
						{
							double rrr[3];
							register double rrr2[3];
							for( int iii=ql; iii<qr; ++iii )
							{       
								rrr[0] = rr[0]+(double)iii*dx05 - dx025;
								rrr2[0]= rrr[0]*rrr[0];
								for( int jjj=ql; jjj<qr; ++jjj )
								{
									rrr[1] = rr[1]+(double)jjj*dx05 - dx025;
									rrr2[1]= rrr[1]*rrr[1];
									for( int kkk=ql; kkk<qr; ++kkk )
									{
										rrr[2] = rr[2]+(double)kkk*dx05 - dx025;
										rrr2[2]= rrr[2]*rrr[2];
										rr2 = rrr2[0]+rrr2[1]+rrr2[2];
										val += tfr->compute_real(rr2)/rf8;
									}
								}
							}
						}
						else
						{
							rr2 = rr[0]*rr[0]+rr[1]*rr[1]+rr[2]*rr[2];
							val = tfr->compute_real(rr2);	
						}
						
						
						//if( rr2 <= boxlength2*boxlength2 )
						//rkernel[idx] += (fftw_real)tfr->compute_real(rr2)*fac;
						val *= fac;	
						
						for(int q=0;q<8;++q)
							if(idx[q]!=(size_t)-1)  
								rkernel[idx[q]] = val;
						
					}
		}
		
		rkernel[0] = tfr->compute_real(0.0)*fac;
		
		/*************************************************************************************/
		/*************************************************************************************/
		/******* perform deconvolution *******************************************************/
		
		
		
		//bool baryons = type==baryon||type==vbaryon;
		if( deconv )
		{
				
			LOGUSER("Deconvolving fine kernel...");
			std::cout << " - Deconvolving density kernel...\n";
			
			
			
			double fftnorm = 1.0/((size_t)nx*(size_t)ny*(size_t)nz);
			double k0 = rkernel[0];
			
			fftw_complex	*kkernel	= reinterpret_cast<fftw_complex*>( &rkernel[0] );
			
			//... subtract white noise component before deconvolution
			if(!bsmooth_baryons)
				rkernel[0] = 0.0;
			
#ifdef FFTW3
	#ifdef SINGLE_PRECISION
			fftwf_plan 
				plan  = fftwf_plan_dft_r2c_3d(nx,ny,nz, rkernel, kkernel, FFTW_ESTIMATE),
				iplan = fftwf_plan_dft_c2r_3d(nx,ny,nz, kkernel, rkernel, FFTW_ESTIMATE);			
			
			fftwf_execute(plan);
	#else
			fftw_plan 
			plan  = fftw_plan_dft_r2c_3d(nx,ny,nz, rkernel, kkernel, FFTW_ESTIMATE),
			iplan = fftw_plan_dft_c2r_3d(nx,ny,nz, kkernel, rkernel, FFTW_ESTIMATE);			
			
			fftw_execute(plan);
	#endif
#else
			rfftwnd_plan	plan		= rfftw3d_create_plan( nx,ny,nz, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE|FFTW_IN_PLACE),
							iplan		= rfftw3d_create_plan( nx,ny,nz, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE|FFTW_IN_PLACE);
			
	#ifndef SINGLETHREAD_FFTW		
			rfftwnd_threads_one_real_to_complex( omp_get_max_threads(), plan, rkernel, NULL );
	#else
			rfftwnd_one_real_to_complex( plan, rkernel, NULL );
	#endif
#endif		
			
			
			
			
			if( deconv )
			{
				
				double ksum = 0.0;
				size_t kcount = 0;
				double kmax = 0.5*M_PI/std::max(nx,std::max(ny,nz));
				
				#pragma omp parallel for reduction(+:ksum,kcount)
				for( int i=0; i<nx; ++i )
					for( int j=0; j<ny; ++j )
						for( int k=0; k<nz/2+1; ++k )
						{
							double kx,ky,kz;
							
							kx = (double)i;
							ky = (double)j;
							kz = (double)k;
							
							if( kx > nx/2 ) kx -= nx;
							if( ky > ny/2 ) ky -= ny;
							
							
							double kkmax = kmax;
							size_t q  = ((size_t)i*ny+(size_t)j)*(size_t)(nz/2+1)+(size_t)k;
							
							if( !bsmooth_baryons )
							{
								if( kspacepoisson )
								{
									//... Use child average response function to emulate sub-sampling
									double ipix = cos(kx*kkmax)*cos(ky*kkmax)*cos(kz*kkmax);
									
									RE(kkernel[q]) /= ipix;
									IM(kkernel[q]) /= ipix;
									
								}else{
									
									//... Use piecewise constant response function (NGP-kernel)
									//... for finite difference methods
									kkmax = kmax;
									double ipix = 1.0;
									if( i > 0 )
										ipix /= sin(kx*2.0*kkmax)/(kx*2.0*kkmax);
									if( j > 0 )
										ipix /= sin(ky*2.0*kkmax)/(ky*2.0*kkmax);
									if( k > 0 )
										ipix /= sin(kz*2.0*kkmax)/(kz*2.0*kkmax);
									
									RE(kkernel[q]) *= ipix;
									IM(kkernel[q]) *= ipix;
									
								}
								
							}
#if 1					
							else{
								//... if smooth==true, convolve with 
								//... NGP kernel to get CIC smoothness
								
								//kkmax *= 2.0;
								
								double ipix = 1.0;
								if( i > 0 )
									ipix /= sin(kx*2.0*kkmax)/(kx*2.0*kkmax);
								if( j > 0 )
									ipix /= sin(ky*2.0*kkmax)/(ky*2.0*kkmax);
								if( k > 0 )
									ipix /= sin(kz*2.0*kkmax)/(kz*2.0*kkmax);
								
								RE(kkernel[q]) /= ipix;
								IM(kkernel[q]) /= ipix;

							}
#endif
							
							//... store k-space average
							if( k==0 || k==nz/2 )
							{
								ksum  += RE(kkernel[q]);
								kcount++;
							}else{
								ksum  += 2.0*(RE(kkernel[q]));
								kcount+=2;
							}
							
						}
				
				double dk;
				
				//... re-add white noise component for finest grid
				dk = k0-ksum/kcount;
				
				//... set white noise component to zero if smoothing is enabled
				//if( false )//cparam_.smooth )
				if( bsmooth_baryons )
					dk = 0.0;
				
				
				
				//... enforce the r=0 component by adjusting the k-space mean
	#pragma omp parallel for reduction(+:ksum,kcount)
				for( int i=0; i<nx; ++i )
					for( int j=0; j<ny; ++j )
						for( int k=0; k<(nz/2+1); ++k )
						{
							size_t q  = ((size_t)i*ny+(size_t)j)*(nz/2+1)+(size_t)k;
							
							RE(kkernel[q]) += dk;
							
							RE(kkernel[q]) *= fftnorm;
							IM(kkernel[q]) *= fftnorm;
						}
			}		
			
			
			
#ifdef FFTW3
	#ifdef SINGLE_PRECISION
			fftwf_execute(iplan);
			fftwf_destroy_plan(plan);
			fftwf_destroy_plan(iplan);
	#else
			fftw_execute(iplan);
			fftw_destroy_plan(plan);
			fftw_destroy_plan(iplan);
	#endif
#else
	#ifndef SINGLETHREAD_FFTW		
			rfftwnd_threads_one_complex_to_real( omp_get_max_threads(), iplan, kkernel, NULL );
	#else
			rfftwnd_one_complex_to_real( iplan, kkernel, NULL );
	#endif
			rfftwnd_destroy_plan(plan);
			rfftwnd_destroy_plan(iplan);
#endif
		}
		
		/*************************************************************************************/
		/*************************************************************************************/
		/*************************************************************************************/
		
		
		char cachefname[128];
		sprintf(cachefname,"temp_kernel_level%03d.tmp",levelmax);
		LOGUSER("Storing kernel in temp file \'%s\'.",cachefname);

		FILE *fp = fopen(cachefname,"w+");
		unsigned q = nx;
		fwrite(	reinterpret_cast<void*> (&q), sizeof(unsigned), 1, fp );
		q = ny;
		fwrite(	reinterpret_cast<void*> (&q), sizeof(unsigned), 1, fp );
		q = 2*(nz/2+1);
		fwrite(	reinterpret_cast<void*> (&q), sizeof(unsigned), 1, fp );
		
		
		for( int ix=0; ix<nx; ++ix )
		{
			size_t sz = ny*2*(nz/2+1);
			//fwrite( reinterpret_cast<void*>(&rkernel[0]), sizeof(fftw_real), nx*ny*2*(nz/2+1), fp );			
			fwrite( reinterpret_cast<void*>(&rkernel[(size_t)ix * sz]), sizeof(fftw_real), sz, fp );			
		}

		
		
		
		fclose(fp);
		
		//... average and fill for other levels
		for( int ilevel=levelmax-1; ilevel>=levelmin; ilevel-- )
		{
			LOGUSER("Computing coarse kernel (level %d)...",ilevel);
			
			int nxc, nyc, nzc;
			double dxc,lxc,lyc,lzc;
			
			nxc = refh.size(ilevel,0);
			nyc = refh.size(ilevel,1);
			nzc = refh.size(ilevel,2);
			
			if( ilevel!=levelmin )
			{
				nxc *= 2;
				nyc *= 2;
				nzc *= 2;
			}
			
			
			dxc = boxlength / (1<<ilevel);
			lxc = dxc * nxc;
			lyc = dxc * nyc;
			lzc = dxc * nzc;
			
			rkernel_coarse = new fftw_real[(size_t)nxc*(size_t)nyc*2*((size_t)nzc/2+1)];
			fac			= lxc*lyc*lzc/pow(2.0*M_PI,3)/((double)nxc*(double)nyc*(double)nzc);
			
			if( bperiodic  )
			{		
				#pragma omp parallel for
				for( int i=0; i<=nxc/2; ++i )
					for( int j=0; j<=nyc/2; ++j )
						for( int k=0; k<=nzc/2; ++k )
						{
							int iix(i), iiy(j), iiz(k);
							double rr[3], rr2;
							
							if( iix > (int)nxc/2 ) iix -= nxc;
							if( iiy > (int)nyc/2 ) iiy -= nyc;
							if( iiz > (int)nzc/2 ) iiz -= nzc;
							
							//... speed up 8x by copying data to other octants
							size_t idx[8];
							
							idx[0] = ((size_t)(i)*nyc + (size_t)(j)) * 2*(nzc/2+1) + (size_t)(k);
							idx[1] = ((size_t)(nxc-i)*nyc + (size_t)(j)) * 2*(nzc/2+1) + (size_t)(k);
							idx[2] = ((size_t)(i)*nyc + (size_t)(nyc-j)) * 2*(nzc/2+1) + (size_t)(k);
							idx[3] = ((size_t)(nxc-i)*nyc + (size_t)(nyc-j)) * 2*(nzc/2+1) + (size_t)(k);
							idx[4] = ((size_t)(i)*nyc + (size_t)(j)) * 2*(nzc/2+1) + (size_t)(nzc-k);
							idx[5] = ((size_t)(nxc-i)*nyc + (size_t)(j)) * 2*(nzc/2+1) + (size_t)(nzc-k);
							idx[6] = ((size_t)(i)*nyc + (size_t)(nyc-j)) * 2*(nzc/2+1) + (size_t)(nzc-k);
							idx[7] = ((size_t)(nxc-i)*nyc + (size_t)(nyc-j)) * 2*(nzc/2+1) + (size_t)(nzc-k);
							
							if(i==0||i==nxc/2){ idx[1]=idx[3]=idx[5]=idx[7]=(size_t)-1;}
							if(j==0||j==nyc/2){ idx[2]=idx[3]=idx[6]=idx[7]=(size_t)-1;}
							if(k==0||k==nzc/2){ idx[4]=idx[5]=idx[6]=idx[7]=(size_t)-1;}
							
							double val = 0.0;
							
							for( int ii=-1; ii<=1; ++ii )
								for( int jj=-1; jj<=1; ++jj )
									for( int kk=-1; kk<=1; ++kk )
									{
										rr[0] = ((double)iix ) * dxc + ii*boxlength;
										rr[1] = ((double)iiy ) * dxc + jj*boxlength;
										rr[2] = ((double)iiz ) * dxc + kk*boxlength;
										
										if( rr[0] > -boxlength && rr[0] < boxlength
										   && rr[1] > -boxlength && rr[1] < boxlength
										   && rr[2] > -boxlength && rr[2] < boxlength )
										{
											rr2 = rr[0]*rr[0]+rr[1]*rr[1]+rr[2]*rr[2];
											val += tfr->compute_real(rr2);	
										}
									}
							
							val *= fac;
							
							for(int qq=0;qq<8;++qq)
								if(idx[qq]!=(size_t)-1)  
									rkernel_coarse[idx[qq]] = val;	
						}
				
			}else{
				#pragma omp parallel for
				for( int i=0; i<nxc; ++i )
					for( int j=0; j<nyc; ++j )
						for( int k=0; k<nzc; ++k )
						{
							int iix(i), iiy(j), iiz(k);
							double rr[3], rr2;
							
							if( iix > (int)nxc/2 ) iix -= nxc;
							if( iiy > (int)nyc/2 ) iiy -= nyc;
							if( iiz > (int)nzc/2 ) iiz -= nzc;
							
							size_t idx = ((size_t)i*nyc + (size_t)j) * 2*(nzc/2+1) + (size_t)k;
							
							rr[0] = ((double)iix ) * dxc;
							rr[1] = ((double)iiy ) * dxc;
							rr[2] = ((double)iiz ) * dxc;
							
							rkernel_coarse[idx] = 0.0;
							
							rr2 = rr[0]*rr[0]+rr[1]*rr[1]+rr[2]*rr[2];
							if( fabs(rr[0])<=boxlength2||fabs(rr[1])<=boxlength2||fabs(rr[2])<=boxlength2 )
								rkernel_coarse[idx] += (fftw_real)tfr->compute_real(rr2)*fac;
							
						}
			}

			LOGUSER("Averaging fine kernel to coarse kernel...");

			//... copy averaged and convolved fine kernel to coarse kernel
			#pragma omp parallel for
			for( int ix=0; ix<nx; ix+=2 )
				for( int iy=0; iy<ny; iy+=2 )
					for( int iz=0; iz<nz; iz+=2 )
					{
						int iix(ix/2),iiy(iy/2),iiz(iz/2);
						if( ix>nx/2 ) iix+=nxc-nx/2;
						if( iy>ny/2 ) iiy+=nyc-ny/2;
						if( iz>nz/2 ) iiz+=nzc-nz/2;
						
						if( ix==nx/2||iy==ny/2||iz==nz/2 ) continue;
						
						for( int i=0; i<=1; ++i )
							for( int j=0; j<=1; ++j )
								for( int k=0; k<=1; ++k )
									if(i==0&&k==0&&j==0 )
										rkernel_coarse[ ACC_RC(iix,iiy,iiz) ] =
										0.125 * ( rkernel[ACC_RF(ix-i,iy-j,iz-k)] + rkernel[ACC_RF(ix-i+1,iy-j,iz-k)] 
												 +rkernel[ACC_RF(ix-i,iy-j+1,iz-k)] + rkernel[ACC_RF(ix-i,iy-j,iz-k+1)]
												 +rkernel[ACC_RF(ix-i+1,iy-j+1,iz-k)] + rkernel[ACC_RF(ix-i+1,iy-j,iz-k+1)]
												 +rkernel[ACC_RF(ix-i,iy-j+1,iz-k+1)] + rkernel[ACC_RF(ix-i+1,iy-j+1,iz-k+1)]);
						
									else
									{	
										
										rkernel_coarse[ ACC_RC(iix,iiy,iiz) ] +=
										0.125 * ( rkernel[ACC_RF(ix-i,iy-j,iz-k)] + rkernel[ACC_RF(ix-i+1,iy-j,iz-k)] 
												 +rkernel[ACC_RF(ix-i,iy-j+1,iz-k)] + rkernel[ACC_RF(ix-i,iy-j,iz-k+1)]
												 +rkernel[ACC_RF(ix-i+1,iy-j+1,iz-k)] + rkernel[ACC_RF(ix-i+1,iy-j,iz-k+1)]
												 +rkernel[ACC_RF(ix-i,iy-j+1,iz-k+1)] + rkernel[ACC_RF(ix-i+1,iy-j+1,iz-k+1)]);
									}
						
					}

			
			sprintf(cachefname,"temp_kernel_level%03d.tmp",ilevel);
			LOGUSER("Storing kernel in temp file \'%s\'.",cachefname);
			fp = fopen(cachefname,"w+");
			q = nxc;
			fwrite(	reinterpret_cast<void*> (&q), sizeof(unsigned), 1, fp );
			q = nyc;
			fwrite(	reinterpret_cast<void*> (&q), sizeof(unsigned), 1, fp );
			q = 2*(nzc/2+1);
			fwrite(	reinterpret_cast<void*> (&q), sizeof(unsigned), 1, fp );
			
			for( int ix=0; ix<nxc; ++ix )
			{
				size_t sz = nyc*2*(nzc/2+1);
				//fwrite( reinterpret_cast<void*>(&rkernel_coarse[0]), sizeof(fftw_real), nxc*nyc*2*(nzc/2+1), fp );
				fwrite( reinterpret_cast<void*>(&rkernel_coarse[ix*sz]), sizeof(fftw_real), sz, fp );
			}
			
			fclose(fp);
			
			delete[] rkernel;
			
			//... prepare for next iteration
			nx = nxc;
			ny = nyc;
			nz = nzc;
			lx = lxc;
			ly = lyc;
			lz = lzc;
			dx = dxc;
			rkernel = rkernel_coarse;
		}
		
		//... clean up
		delete[] rkernel;
	}
	
}


/**************************************************************************************/
/**************************************************************************************/


namespace
{
	convolution::kernel_creator_concrete< convolution::kernel_real_cached<double> > creator_d("tf_kernel_real_double");
	convolution::kernel_creator_concrete< convolution::kernel_real_cached<float> > creator_f("tf_kernel_real_float");
	convolution::kernel_creator_concrete< convolution::kernel_k<double> > creator_kd("tf_kernel_k_double");
	convolution::kernel_creator_concrete< convolution::kernel_k<float> > creator_kf("tf_kernel_k_float");
}