#include <math.h>
#include <stdio.h> 
/* compile with
g++ fft.cxx  /cern/pro/lib/libmathlib.a  -lm -lg2c -o fft 

*/
extern "C" void cft_(float *,float *,int&,int&,int&,int&);
 
inline float sqr(float x) { return x*x;} 

int main(int argc, char **argv){
  int a=0,b=0,c=0;
  if ( argc < 1 ) {
      printf("Uporaba: %s <file.dat>\n",argv[0]);
      exit(1);
     
  }
  //-------------------------------------------------------
  //   preberi file
  //-------------------------------------------------------
  FILE *fp=fopen(argv[1],"r");
  float f[10];
  int j=0,i=0,n=0;   
  int ndim  =  400;
  char line[ndim];

#define NDIM 10000
  float x1[NDIM] ,y1[NDIM],x[NDIM],y[NDIM];                                           
  while (fgets(line,ndim,fp)!=NULL) {
    if (sscanf(line,"%f%f%f%f%f",&f[0],&f[1],&f[2],&f[3],&f[4]) != 5) continue;
    if ( j> NDIM ) {
      printf("Povecaj NDIM in se enkrat prevedi program %s\n",argv[0]);
      exit(1);
    }
    x1[j]=f[0];
    y1[j]=f[4];     
    x[j]=f[4];     
    y[j]=0;     
    
    j++;
  }
  fclose(fp);
  //-------------------------------------------------------
  //   FFT
  //-------------------------------------------------------

  n = 1;
  cft_(x,y,j,j,j,n);

  //-------------------------------------------------------
  //   poisci prvi minimum in maksimalno frekvenco
  //-------------------------------------------------------

  float max = 0 ;
  int   imax= 0;
  float xobrat = x1[0];
  float yobrat = -1;
  int   iobrat = 0;
  max = sqr(x[0])+sqr(y[0]);
  int first = 1;
  // zapeljemo se samo do polovice spektra, ker je druga polovica simetricna
  for (i=1;i<(j/2);i++){
    float pow  = sqr(x[i])+sqr(y[i]);
    float dy   = y1[i] - y1[i-1];
    if (pow > max) { 
      max = pow ; 
      imax = i; 
    } 
    if (first) {
      if ( yobrat * dy > 0 ) { 
	yobrat = dy ;iobrat = i; 
      } else {
	first = 0;
      }
    }
  }
  yobrat = y1[iobrat];
  xobrat = x1[iobrat];
  //-------------------------------------------------------
  //   zamenjaj . _ in : v imenu s presledkom
  //-------------------------------------------------------
  char fname[100];
  sprintf(fname,"%s",argv[1]);
  for (int j = 0; j < 100 && fname[j] != '\0' ;j++) {
    fname[j] = (fname[j] == '.') ? ' ' : fname[j];
    fname[j] = (fname[j] == '_') ? ' ' : fname[j];
    fname[j] = (fname[j] == ':') ? ' ' : fname[j];
  }

  //-------------------------------------------------------
  //   izpisi rezultat
  //-------------------------------------------------------

  printf("%s\t%f\t%f\t%d\t%f\n",fname,xobrat,yobrat, imax ,max);
}