fft.c

/*
 * BHAC tools
 * Copyright 2024
 * Author: Jordy Davelaar
 */

#include "definitions.h"
#include "functions.h"
#include "global_vars.h"
#include "model_definitions.h"
#include "model_functions.h"
#include "model_global_vars.h"
#include <fftw3.h>
#include <stdlib.h>
#include <math.h>
#include <omp.h>
void perform_fft(double TIME_INIT){
  	fftw_plan p;

  	fftw_complex *quant;
  	fftw_complex *outx, *outy, *outz ;

	int ind;

	double *unif;

	double intg[UNI2];
	double freq[UNI2];

	quant = (fftw_complex *) fftw_malloc(UNI1*UNI2*UNI3*sizeof(fftw_complex));
	outx = (fftw_complex *) fftw_malloc(UNI1*UNI2*UNI3*sizeof(fftw_complex));
	outy = (fftw_complex *) fftw_malloc(UNI1*UNI2*UNI3*sizeof(fftw_complex));
	outz = (fftw_complex *) fftw_malloc(UNI1*UNI2*UNI3*sizeof(fftw_complex));
	unif = (double*)malloc(UNI1*UNI2*UNI3*sizeof(double));

	uniform_data_array(unif,5);

        get_quant(unif, quant );

	fftw_init_threads();
	int current_num_threads = omp_get_num_threads();

	fftw_plan_with_nthreads(current_num_threads);

  	/* forward Fourier transform, save the result in 'out' */
#if(THREED)
     	p = fftw_plan_dft_3d(UNI1, UNI2, UNI3, quant, outx, FFTW_FORWARD, FFTW_ESTIMATE);
#else
  	p = fftw_plan_dft_2d(UNI1, UNI2, quant, outx, FFTW_FORWARD, FFTW_ESTIMATE);
#endif
  	fftw_execute(p);

        uniform_data_array(unif,6);

        get_quant(unif, quant );

        /* forward Fourier transform, save the result in 'out' */
#if(THREED)
	p = fftw_plan_dft_3d(UNI1, UNI2, UNI3, quant, outy, FFTW_FORWARD, FFTW_ESTIMATE);
#else
     	p = fftw_plan_dft_2d(UNI1, UNI2, quant, outy, FFTW_FORWARD, FFTW_ESTIMATE);
#endif
      	fftw_execute(p);

        uniform_data_array(unif,7);

	get_quant(unif, quant );

        /* forward Fourier transform, save the result in 'out' */
#if(THREED)
	p = fftw_plan_dft_3d(UNI1, UNI2, UNI3, quant, outz, FFTW_FORWARD, FFTW_ESTIMATE);
#else
     	p = fftw_plan_dft_2d(UNI1, UNI2, quant, outz, FFTW_FORWARD, FFTW_ESTIMATE);
#endif
      	fftw_execute(p);

	fftfreq(freq,UNI2,YMAX-YMIN);

	for(int i=0;i<UNI2; i++)
		intg[i]=0;
double prod;
#pragma omp parallel for shared(intg,quant) private(prod) schedule(static, 1)
        for(int i=0;i<UNI1; i++){
	        for(int j=0;j<UNI2; j++){
        		for(int k=0;k<UNI3; k++){
		           //intg[i]+=(fabs(outx[i*UNI2*UNI3+j*UNI3+k][0])+fabs(outy[i*UNI2*UNI3+j*UNI3+k][0])+fabs(outz[i*UNI2*UNI3+j*UNI3+k][0]))/((double) (UNI2*UNI3));
#if(THREED)
                               	ind = i*UNI2*UNI3+j*UNI2+k;
#else
                               	ind = i*UNI2+j;
#endif
				prod = outx[ind][0]*outx[ind][0]+outx[ind][1]*outx[ind][1];
				prod += outy[ind][0]*outy[ind][0]+outy[ind][1]*outy[ind][1];
				prod += outz[ind][0]*outz[ind][0]+outz[ind][1]*outz[ind][1];
				intg[j]+=prod/(UNI1*UNI3);
			}
		}
        }


	write_fft_file(freq,intg,TIME_INIT);

	fftw_free(outx);
	fftw_free(outy);
	fftw_free(outz);
	fftw_free(quant);

  	fftw_destroy_plan(p);

}


void fftfreq(double* freq, int N, double size){
	if(N%2==0){
		for(int i=0;i<N/2;i++){
			freq[i] = i/size;
		}

		for(int i=N/2;i<N;i++){
			freq[i] = (-N/2 + (i-N/2))/size;
		}
	}
	else{
		for(int i=0;i<=(N-1)/2;i++){
			freq[i] = i/size;
		}

		for(int i=(N-1)/2;i<N;i++){
			freq[i] = (-(N-1)/2 + (i-(N-1)/2))/size;
		}

	}
}

double window_func(double x){
	double k=3.0;
	double m=8.0;
	double enumarator = k * (fabs(x) - 0.5 * ( XMAX + XMIN));
	double denom = XMAX - XMIN;

	return  exp(-powf(enumarator/denom,m));

}

void get_quant(double* unif, fftw_complex *quant ){
	double x;
	int ind;

#pragma omp parallel for shared(quant) private(x,ind) schedule(static, 1)
        for(int i=0;i<UNI1; i++){
               	for(int j=0;j<UNI2; j++){
#if(THREED)
                        for(int k=0;k<UNI3; k++){
#endif
	                       	x =  (XMAX-XMIN)/((double) UNI1) * i + XMIN + 0.5 * (XMAX-XMIN)/((double) UNI1);

#if(THREED)
                               	ind = i*UNI2*UNI3+j*UNI2+k;
#else
                               	ind = i*UNI2+j;
#endif

        	               	quant[ind][0]=unif[ind]*window_func(x);
                	        quant[ind][1]=0.0;
#if(THREED)
 		      	}
#endif
		}
	}
}