/* standard begin
*/

#include <stdio.h>
#include <spec.h>

#define  MAXPEAKS 10

float x,y,*spc,*err,*tim,*tr,*ti;

help()
{
  printf("Fast fourier transformation from input spectrum\n");
  printf("The resulting spectrum is allready a power spectrum\n");
  printf("option     -iv        inverse spectrum\n");
  printf("           -rf name   save real part to file\n");
  printf("           -if name   save imaginary part to file\n");
  printf("           -pc n      append n peaks to comment\n");
  printf("           -pp n      print n peaks to stdout\n");
  printf("           -s         suppress printing of powerspectrum to stdout\n");
}

main(argc,argv)
int argc;
char *argv[];
{
int kmax2,n,max,inv,peaks[MAXPEAKS];
char z[1024],comment[1024];

   inv=0;
   if(checkopt(argc,argv,"-iv",z)) inv=1;

   tr = (float *)calloc(_MAXSPCLEN+2,sizeof(float));
   ti = (float *)calloc(_MAXSPCLEN+2,sizeof(float));
   if(ti==NULL) {
      printf("sorry, not enough memory\n");
      exit(-1);
   }
   for(n=0;n<=_MAXSPCLEN;n++) {  /* should not be neccessary, but good style */
      tr[n]=0; ti[n]=0;
   }
   spc= (float *)calloc(_MAXSPCLEN+2,sizeof(float));
   err= (float *)calloc(_MAXSPCLEN+2,sizeof(float));
   tim= (float *)calloc(_MAXSPCLEN+2,sizeof(float));
   if(tim==NULL) {
      printf("sorry, not enough memory\n");
      exit(-1);
   }

   max=readspec(argv[1],spc,err,tim,comment);
   for(n=0;n<=max;n++) tr[n]=spc[n];
   kmax2=four(inv,max);
   for(n=0;n<=max;n++) {         /* now recalculate time spectrum */
      x = (float) n;
      tim[n] = (6.2831853 * ((float) n)) / (((float) kmax2) * _tica);
   }
   for(n=0;n<=max;n++) spc[n]=sqrt((tr[n]*tr[n])+(ti[n]*ti[n]));
   fpeaks(peaks,MAXPEAKS,spc,kmax2/2);
   if(checkopt(argc,argv,"-pc",z)) {
      n=atoi(z);
      peaktoa(z,n,peaks,kmax2,spc);
      strcat(comment,": ");
      strcat(comment,z);
   }
   if(!checkopt(argc,argv,"-s",z)) writespec("",spc,err,kmax2/2,2,comment);
   if(checkopt(argc,argv,"-rf",z)) writespec(z,tr,err,kmax2/2,2,comment);
   if(checkopt(argc,argv,"-if",z)) writespec(z,ti,err,kmax2/2,2,comment);
   if(checkopt(argc,argv,"-pp",z)) {
      n=atoi(z);
      peaktoa(z,n,peaks,kmax2,spc);
      printf("%s\n",z);
   }
   if(checkopt(argc,argv,"-o",z)) {            /* write modified time array */
      strcat(z,".tim");
      writespec(z,tim,err,kmax2/2,2,comment);
   }
   free(tr); free(ti); free(spc); free(err); free(tim);
   exit(0);
}

/* ------------------------------------------------------------------------
   Fast Fourier transform
   ------------------------------------------------------------------------ */
log2(n)
int n;
{
int erg;

   erg=0;
   while((1<<erg)<n) erg++;
   return(erg);
}

four(inv,max)
int inv;
{
int   k,ii,i,j,it,i0,i1,i2,i3,
      kk,j0,kk1,l,l1,jp1,kp1,j1,j2,
      kmax2;
float z,zr,zi,wr,wi,ws,ur[15],ui[15],um,th;

   it=max;
   um=0.5;
   for(i=1;i<=15;i++) {
      um = 0.5 * um;
      th = 6.283185 * um;
      ur[i] = cos(th); ui[i]=sin(th);
   }
   um = 1 - (2*inv);
   i0 = log2(it); kmax2 = 1 << i0;
   it = kmax2;
   ii=i0 ; i1=it/2; i3=1;
   do {
      k=0; i2=i1+i1;
      do {
         wr=1; wi=0; kk=k; j0=i0;
         while(kk!=0) {
            do {
               j0=j0-1; kk1=kk; kk=kk/2;
            } while((kk1-(2*kk)) == 0);
            ws=(wr*ur[j0]) - (wi*ui[j0]);
            wi=(wr*ui[j0]) + (wi*ur[j0]);
            wr = ws;
         }
         wi=wi*um; j=0;
         do {
            l=j*i2+k+1; l1=l+i1;
            zr=tr[l]+tr[l1];
            zi=ti[l]+ti[l1];
            z=wr*(tr[l]-tr[l1]) - (wi*(ti[l]-ti[l1]));
            ti[l1]=wr*(ti[l]-ti[l1]) + (wi*(tr[l]-tr[l1]));
            tr[l]=zr; tr[l1]=z; ti[l]=zi; j=j+1;
         } while((j-i3)<0);
         k=k+1;
      } while((k-i1)<0);
      i3 = i3 + i3; i0 = i0 - 1; i1 = i1 / 2;
   } while(i1 > 0);
   j=1; um=1;
   if(inv==1) um=1/it;
   do {
      k=0; j1=j;
      for(i=1;i<=ii;i++) {
         j2 = j1 / 2; k=(2*(k-j2))+j1; j1=j2;
      }
      jp1=j+1;
      if((k-j) >= 0) {
         if((k-j) == 0) {
            tr[jp1] = um * tr[jp1];
            ti[jp1] = um * ti[jp1];
         } else {
            kp1 = k + 1;
            zr = tr[jp1]; zi=ti[jp1];
            tr[jp1] = um * tr[kp1];
            ti[jp1] = um * ti[kp1];
            tr[kp1] = zr * um;
            ti[kp1] = zi * um;
         }
      }
      j = jp1;
   } while((j-it+1)<0);
   tr[1] = um * tr[1];
   ti[1] = um * ti[1];
   return(kmax2);
}

/* -----------------------------------------------------------------------
   find peaks in powerspectrum
   ----------------------------------------------------------------------- */

fpeaks(peaks,pmax,spc,smax)
float spc[];
int peaks[],pmax,smax;
{
int n,sptr;

   sptr=0;
   for(n=0;n<=pmax;n++) peaks[n]= -1 ;             /* mark peaks as unset */
   for(;;) {
      while(spc[sptr] <= spc[sptr+1]) {             /* find next maximum */
         sptr=sptr+1;
         if(sptr>smax) return(0);
      }
      peaks[pmax]=sptr;                            /* smallest peak redefined */
      peaksort(peaks,pmax,spc,smax);               /* sort peaks for height  */
      while(spc[sptr] > spc[sptr+1]) {             /* find next minimum */
         sptr=sptr+1;
         if(sptr>smax) return(0);
      }
   }
}

peaksort(peaks,pmax,spc,smax)
float spc[];
int peaks[],pmax,smax;
{
int m,p1,p2,flag;

   flag=0;
   if(peaks[pmax]==-1) return(0);
   while(flag == 0) {
      flag=1;
      for(m=0;m<pmax;m++) {
         p1=peaks[m];
         p2=peaks[m+1];
         if(p1 == -1) {
            peaks[m]=pmax;
	    flag=0;
            continue;
         } 
         if(p2 == -1) {
	    peaks[m+1]=pmax;
	    flag=0;
	    continue;
         }
         if(spc[p1]<spc[p2]) {
            peaks[m]=p2;
            peaks[m+1]=p1;
	    flag=0;
         }
      }
   }
}

/* -------------------------------------------------------------- */

float peakmean(peak,spc)
int peak;
float spc[];
{
int n,left,right;
float x,y,y1,y2,sum,lsum,mean;

   left=peak;                              /* finding left side of peak */
   for(n=0;n<=20;n++) {
      if(left==0) break;
      y1=spc[left];
      y2=spc[left-1];
      if(y2>y1) break;
      y = y1 / 10.0;
      x = fabs(y1-y2);
      if(x<y) break;
      left = left - 1;
   }
   right=peak;                             /* finding right side of peak */
   for(n=0;n<=20;n++) {
      y1=spc[right];
      y2=spc[right+1];
      if(y2>y1) break;
      y = y1 / 10.0;
      x = fabs(y1-y2);
      if(x<y) break;
      right = right + 1;
   }
   sum = 0.0;                                /* peak integral */
   for(n=left;n<=right;n++) sum = sum + spc[n];
   sum = sum / 2;
   lsum = 0.0;
   n=left;
   while(lsum < sum) {                       /* finding half integral */
      y1 = lsum;
      lsum = lsum + spc[n++];
   }
   n = n - 1;
   y1 = y1 - sum;
   y2 = lsum - sum;
   if(y1==y2) {
      x = peak;
      return(x);
   }
   x = y1 / (y1 - y2);
   mean = (n-1) + x;
   return(mean);
}

peaktoa(s,max,peaks,kmax2,spc)
char s[];
int max,peaks[],kmax2;
float spc[];
{
char z[80];
int i,n,m,ix;
float x,x1,x2;

   strcpy(s,"");
   for(n=0;n<max;n++) {
      i=peaks[n];                                 /* get maximum */
      x=peakmean(i,spc);                          /* mean instead of maximum */
      ix=x;
      x1=tim[ix]; x2=tim[ix+1];                   /* interpolation between    */
      x = x1 + ((x-ix) * (x2-x1));                /* tim[i] and tim[i+1]      */
      sprintf(z," %E",x);
      if(fabs(x)>0.001) sprintf(z," %-1.5f",x);
      if(fabs(x)>1.0) sprintf(z," %-3.2f",x);
      for(m=0;m<n;m++) {
	 if(peaks[m] == i) return(0);
      }
      strcat(s,z);
   }
}