/*********************************************************************
*
*  Task FFT
*
* Fast Fourier Transform
* See Ronald N. Bracewell "The Fourier Transform and its Applications"
*
* CODE 0 -> FFT
* CODE 1 -> IFFT
*
* Default outclass = FFT
*
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <signal.h>

#include "TISAN.H"

void InitializeFFT(void);
void FOUR1(int);

char BUFFER[sizeof(struct TRData)];
char INFILE[_MAX_PATH], OUTFILE[_MAX_PATH];
char TMPFILE[_MAX_PATH];
unsigned long N=0L;
double Imean=0., Rmean=0.;
FILE *INSTR, *OUTSTR;
struct FILEHDR FileHeader;
short ERRFLAG=0, NUMDAT, NDAT;
struct RData  *RDataPntr;
struct XData  *XDataPntr;
struct XData  XDataOut;
double Nyquist, Fundamental;
double huge Z[32770];

main(argc,argv)
short argc;
char *argv[];
   {
   unsigned int wI;

   signal(SIGINT,BREAKREQ);

   if (Ztskinit("FFT",argv[0])) Zexit(1);
   Zencode(2,"FFT     : This is a Non-Standard Program\n");
   if ((CODE < 0) || (CODE > 1))
      {
      Zencode(10,"FFT     : Invalid CODE Specification\n");
      Zexit(1);
      }

   ZBuildFileName(M_inname,INFILE);
   if (!OUTCLASS[0]) strcpy(OUTCLASS,"FFT");
   ZBuildFileName(M_outname,OUTFILE);
   ZBuildFileName(M_tmpname,TMPFILE);

   Zencode(2,"FFT     : Opening Input File '%s'\n",INFILE);
   if (((INSTR = Zopen(INFILE,O_readb)) == NULL) ||
       (!Zgethead(INSTR,&FileHeader))) Zexit(1);

   InitializeFFT();

   FOUR1(CODE ? -1 : 1);

   Zencode(2,"FFT     : Opening Scratch File '%s'\n",TMPFILE);
   if ((OUTSTR = Zopen(TMPFILE,O_writeb)) == NULL)
      {
      ERRFLAG=1;
      goto EXIT;
      }

   if (Zputhead(OUTSTR,&FileHeader))
      {
      ERRFLAG=1;
      goto EXIT;
      }

   XDataOut.f = 0;
   for (wI = 1; wI <= (unsigned int)(2L * N); wI += 2)
      {
      XDataOut.z.x = Z[wI];
      XDataOut.z.y = Z[wI+1];

      if (!CODE)
         {
         XDataOut.z.y /= (double)(N / 2L);
         XDataOut.z.x /= (double)(N / 2L);
         }
      else
         {
         XDataOut.z.y /= 2.;
         XDataOut.z.x /= 2.;
         }

     if (Zwrite(OUTSTR,(char *)&XDataOut,X_Data))
         {
         ERRFLAG = 1;
         break;
         }
      }

   if (!ERRFLAG)
      {
      FileHeader.type = X_Data;
      FileHeader.b = 0.;
      FileHeader.m = (2. * Nyquist - Fundamental) / (double)(N - 1L);
      if (Zputhead(OUTSTR,&FileHeader)) ERRFLAG = 1;
      }
EXIT:
   Zclose(OUTSTR);
   ERRFLAG = Znameout(OUTFILE,TMPFILE,ERRFLAG);
   Zexit(ERRFLAG);
   }

/************************************************************
**
** Initialize file stuff
** When Done, We know if the entire file can be processed
** in memory and how many data points there are (N).
*/
void InitializeFFT()
   {
   double RVAL, IVAL=0., TIME, DATA;
   long L, lFLAGGED=0L;
   int FLAG;

   RDataPntr =  (struct RData *)BUFFER;
   XDataPntr =  (struct XData *)BUFFER;
   L = 1L;
   while (Zread(INSTR,BUFFER,FileHeader.type))
      {
      ++N;
      switch (FileHeader.type)
         {
         case R_Data:
            RVAL = RDataPntr->y;
            FLAG = RDataPntr->f;
            break;
         case X_Data:
            RVAL = XDataPntr->z.x;
            IVAL = XDataPntr->z.y;
            FLAG = RDataPntr->f;
            break;
         default:
            Zencode(10,"FFT     : Invalid File Type.\n");
            Zexit(1);
         }

      if (FLAG)
         {
         RVAL = IVAL = 0.0;
         ++lFLAGGED;
         }

      if (N < 16385L)
         {
         Z[L]   = RVAL;
         Z[L+1] = IVAL;
         L += 2L;
         }
      Imean += IVAL;
      Rmean += RVAL;
      } /* End WHILE */

   if (ferror(INSTR)) Zexit(1);

   if ((N < 2L) || (N > 16384L))
      {
      Zencode(10,"FFT     : ** ERROR ** File Contains %lu Point(s)\n",N);
      Zexit(1);
      }


   Zencode(3,"FFT     : File Contains %lu Total Points\n",N);
   Zencode(3,"FFT     : File Contains %lu Flagged Points\n",lFLAGGED);

   L = 1L;
   while (L < N) L *= 2L;
   if (L != N)
      {
      N = L;
      Zencode(10,"FFT     : ** WARNING ** Padding to %lu\n",N);
      }

   Rmean /= (double)N;
   Imean /= (double)N;

   Nyquist     = 1./(2. * FileHeader.m);
   Fundamental = 1./((double)(N - 1L) * FileHeader.m);

   Zencode(3,"FFT     : Mean Value =  %G + i(%G)\n",Rmean,Imean);
   Zencode(3,"FFT     : Fundamental Frequency = %G\n",Fundamental );
   Zencode(3,"FFT     : Nyquist Frequency = %G\n",Nyquist);

   Zclose(INSTR);
   return;
   }

void FOUR1(ISIGN)
int ISIGN;
   {
   double WR,WI,WPR,WPI,WTEMP,THETA, TEMPR, TEMPI;
   unsigned int wN, J, I, M, MMAX, ISTEP;

   wN = 2 * (unsigned int)N;
   J = 1;
   for (I=1; I <= wN; I += 2)
      {
      if (J > I)
         {
         TEMPR = Z[J];
         TEMPI = Z[J+1];
         Z[J]   = Z[I];
         Z[J+1] = Z[I+1];
         Z[I]   = TEMPR;
         Z[I+1] = TEMPI;
         }
      M=wN/2;
      while ((M >= 2 ) && (J > M))
         {
         J = J - M;
         M /= 2;
         }
      J += M;
      }

   MMAX = 2;

   while (wN > MMAX)
      {
      ISTEP = 2 * MMAX;
      THETA = 6.28318530717959 / ((double)ISIGN * (double)MMAX);
      WPR = sin(THETA / 2.0);
      WPR *= (-2.0 * WPR);
      WPI = sin(THETA);
      WR = 1.0;
      WI = 0.0;
      for (M=1; M <= MMAX; M += 2)
         {
         for (I=M; I<=wN; I+=ISTEP)
            {
            J = I + MMAX;
            TEMPR = Z[J]   * WR - Z[J+1] * WI;
            TEMPI = Z[J+1] * WR + Z[J]   * WI;
            Z[J]   = Z[I]   - TEMPR;
            Z[J+1] = Z[I+1] - TEMPI;
            Z[I]   = Z[I]   + TEMPR;
            Z[I+1] = Z[I+1] + TEMPI;
            }
         WTEMP = WR;
         WR = WR * WPR - WI    * WPI + WR;
         WI = WI * WPR + WTEMP * WPI + WI;
         }
       MMAX=ISTEP;
       }
    return;
    }

/***************************************************************
**
** Process ^C Interrupt
*/
short BREAKREQ()
     {
     fcloseall();
     unlink(TMPFILE);
     Zexit(3);
     }