#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"
#include "TComplex.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;
    TVirtualFFT *fft = TVirtualFFT::FFT(3, n, "R2C ES K");
    std::vector<int> ipoint(3);
    //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];
               ipoint[0] = i;
               ipoint[1] = j;
               ipoint[2] = k;
               fft->SetPoint(ipoint.data(), W[i][j][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
   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];
               ipoint[0] = i;
               ipoint[1] = j;
               ipoint[2] = k;
               int index = i*128*128+j*128+k;
               fft->GetPointComplex(ipoint.data(), re[index], Im[index] );
            }
        }
    }
   //fft->GetPointsComplex(re,Im);
    
    //
    TVirtualFFT *fft_back = TVirtualFFT::FFT(3, n, "C2R ES K");
    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];
               ipoint[0] = i;
               ipoint[1] = j;
               ipoint[2] = k;
               int index2=128*128*i+128*j+k;
               fft_back->SetPoint(ipoint.data(),re[index2],Im[index2] );
            }
        }
    }
    fft_back->Transform();
    double *back = new double[128*128*128];
    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];
               ipoint[0] = i;
               ipoint[1] = j;
               ipoint[2] = k;
               int index = i*128*128+j*128+k;
               back[index]=fft_back->GetPointReal(ipoint.data());
            }
        }
    }


   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]=back[i];
    }
   for(int i=0;i<128;i++){
        for(int j=0;j<128;j++){
            double value=0;
            for(int k=0;k<128;k++){
                value+=W_FFT_R[i][j][k];
            }
            density_FFT_XY->SetBinContent(i+1,j+1,value);//Get the two-dimension real
        }
    }
    TCanvas *c1 =new TCanvas("c1","c1",500,500);
    c1->cd();
    gPad->SetLogz();
    density_FFT_XY->Draw("colz");
    c1->SaveAs("XY_FFT_Re.png");
}