#include "TF1.h"
#include "TGraph.h"
#include "TCanvas.h"
#include "TH1.h"
#include "TGraphErrors.h"
#include "TMultiGraph.h"
#include <cmath>
#include "TH2.h"
#include <stdio.h>
#include <iostream>
#include <vector>
#include "TLegend.h"
#include "TH3.h"
#include "TVirtualFFT.h"
using namespace std;

double Pi=3.1415926;

double HernquistMass(double r){
    double rou;
    rou=2*r*2/pow(r+2,3);
    return rou;
}

double MCMC(double x0){
    double y=x0+gRandom->Uniform(-1,1);
    double ratio=min(1.0,HernquistMass(y)/HernquistMass(x0));
    double random =gRandom->Uniform(0,1);
    double nextX;
    if(random<=ratio){
        nextX=y;
    }
    else{
        nextX=x0;
    }
    return nextX;
}

void FFT3(){
    //Create the input data
    TH3D *density = new TH3D("density","density",128,-15,15,128,-15,15,128,-15,15);
    double *Could;
    vector<vector<vector<double>>> W(128,vector<vector<double>>(128,vector<double>(128,0)));
    double costheta,phi,r_x,r_y,r_z,r;
    r=5;

    int num=0,index_x,index_y,index_z;
    for(int i=0;num<300000;i++){
        if(r<=15){
            costheta=gRandom->Uniform(-1,1);
            phi=gRandom->Uniform(0,2*Pi);
            r_x=r*sin(acos(costheta))*cos(phi);
            r_y=r*sin(acos(costheta))*sin(phi);
            r_z=r*costheta;
            density->Fill(r_x,r_y,r_z);
            num++;
        }
        r=MCMC(r);   
    }
    for(int i=0;i<128;i++){
        for(int j=0;j<128;j++){
            for(int k=0;k<128;k++){
                W[i][j][k]=density->GetBinContent(i+1,j+1,k+1);
            }
        }
    }


    int *n= new int[3];//sizes of each dimension
    n[0]=128;
    n[1]=128;
    n[2]=128;
    double *A = new double[128*128*128];//input data
    for(int i=0;i<128;i++){
        for(int j=0;j<128;j++){
            for(int k=0;k<128;k++){
                A[128*128*i+128*j+k]=W[i][j][k];
            }
        }
    }
   TVirtualFFT *fft = TVirtualFFT::FFT(3, n, "R2C ES K");
   fft->SetPoints(A);
   fft->Transform();//FFT
   double *re = new double[128*128*128];//real part
   double *Im = new double[128*128*128];//imaginary part
   fft->GetPointsComplex(re,Im);

   TH2D *density_FFT_XY = new TH2D("density__FFT_XY","FFT_XY",128,-15,15,128,-15,15);
   vector<vector<vector<double>>> W_FFT_R(128,vector<vector<double>>(128,vector<double>(128,0)));//Get the result from double *re
    for(int i=0;i<128*128*128;i++){
        index_x=i/(128*128);
        index_y=(i-128*128*index_x)/128;
        index_z=i-128*128*index_x-128*index_y;
        W_FFT_R[index_x][index_y][index_z]=re[i];
    }
   for(int i=0;i<128;i++){
        for(int j=0;j<128;j++){
            density_FFT_XY->SetBinContent(i+1,j+1,W_FFT_R[i][j][63]);//Get the two-dimension real part at k=64
        }
    }
    TCanvas *c1 =new TCanvas("c1","c1",500,500);
    c1->cd();
    gPad->SetLogz();
    density_FFT_XY->Draw("colz");
    c1->SaveAs("XY_FFT_Re.png");
}