// file: $PDSP/util/compute_fft/v1.0/compute_fft.cc
//
// this program is intended to compute the fft using a user specified algorithm
//

// system include files
//
#include <assert.h>
#include <memory.h>
#include <time.h>
#include <math.h>
#include <stdio.h>

// isip include files
//
#include <fourier_transform.h>

// main program
//

void Metrics(int Nraw,double *xraw, double *Xw){
  
// nive's declarations
//
  float Diff, EDiff, AreaOver=0.0, SumAreaOver=0.0, Overshoot=0.0;
  float SumOvershoot=0.0, h=2*M_PI/N, Areal, Aideal, A, Sreal, Sideal, S;
  float_8 MaxDXw=0.0, DXw=0.0;
  int first=1, xlast;
  double *x;
  FILE *fp;
  assert(sizeof(double) == sizeof(float_8));


  N = ceil(log(Nraw)/log(2));  
  x = calloc (sizeof(double)*N);
  y = calloc (sizeof(double)*N);
  Xw=calloc(sizeof(double)*N);
  memcpy(x,xraw,sizeof(float)*Nraw);
  
  // define default values
  //
  int_4 i,j;

  // create an object
  //
  Fourier_transform ft;

  // set the algorithm (nive-replace 16384 with order N everywhere)
  //
  ft.set_cc(N, 4, 0);   
  ft.compute_cc(y, x);
  
  
  
  for (i = 0; i < (2*N-1); i += 2) {
    j=i/2;
    Xw[j] = sqrt(y[i]*y[i] + y[i+1]*y[i+1]);
    if(j!=0) {
      DXw=fabs(Xw[j]-Xw[j-1]);
      if(DXw>MaxDXw) MaxDXw=DXw;
    }
    
    if(j<(N/2) && Xw[j]>2*M_PI*j/N){
      xlast=j; 
      Diff = Xw[j]-(2*M_PI*j/N);
      EDiff = Xw[j]*Xw[j] - (2*M_PI*j/N)*(2*M_PI*j/N);
      if(first==1){
	SumAreaOver = Diff;
	SumOvershoot = EDiff;
	first=0;

      }
      else{
	if(j==(N/2)-1){
	  SumAreaOver += Diff ; 
	  SumOvershoot += EDiff;
	  AreaOver = (2*M_PI/N)*0.5*SumAreaOver;
	  SumOvershoot = (2*M_PI/N)*0.5*SumOvershoot;
	  
	  if(AreaOver!=0) Overshoot += SumOvershoot/AreaOver;
	}
	else{
	  SumAreaOver += 2*Diff ; 
	 SumOvershoot += 2*EDiff;
	}
      }
    }
    else{
      if(first==0){		
	SumAreaOver -= Diff; 
	AreaOver = (2*M_PI/N)*0.5*SumAreaOver;
	SumOvershoot -= EDiff;
	SumOvershoot = (2*M_PI/N)*0.5*SumOvershoot;
	if(AreaOver!=0) Overshoot += SumOvershoot/AreaOver;
	first=1;
      }
    }  
  }

  xlast++;
  
  for(int i=xlast; i<(N/2); i++) Areal += 2*Xw[i];
  Areal=Xw[xlast-1] + Xw[N/2];
  Areal=h*0.5*Areal;
    
  for(int i=xlast; i<(N/2); i++) Aideal += 2*2*M_PI*i/N;
  Aideal+= 2*2*M_PI*i/N;
  Aideal=h*0.5*Aideal;
  
  A=Aideal-Areal;
  
  for(int i=xlast; i<N/2; i++) Sreal += 2*Xw[i]*Xw[i];
  Sreal=Xw[xlast-1]*Xw[xlast-1] + Xw[N/2]*Xw[N/2];
  Sreal=h*0.5*Sreal;
  
  for(int i=xlast; i<N/2; i++) Sideal += 2*(2*M_PI*i/N)*(2*M_PI*i/N);
  Sideal += 2*(2*M_PI*i/N)*(2*M_PI*i/N);
  Sideal=h*0.5*Sideal;
  
  S=(Sideal-Sreal)/A;



  fp=fopen("Eval.dat",w);
  fprintf(fp,S,OverShoot,MaxDXw);
  
}