/*********************************************************************
*
*  Task DFT
*
* Discrete Fourier Transform
* See Ronald N. Bracewell "The Fourier Transform and its Applications"
*
* CODE 0 -> DFT
* CODE 1 -> IDFT
*
* F(Nu) = 2/N SUM(f(Tau) exp(-i 2 pi Tau (Nu/N))
* summed from 0 to N - 1 (Fundamental to 2*Nyquist)
*
* The Inverse transform i -> -i and 2/N is changed to 1/2
*
* Default outclass = DFT
*
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <signal.h>

#include "TISAN.H"

#define TwoPi 6.28318530718
#define    Pi 3.14159265359

void InitializeDFT(void);
void MEMREDUCE(double);
short DSKREDUCE(double);

char BUFFER[Maxbuff*sizeof(struct TRData)];
char INFILE[_MAX_PATH], OUTFILE[_MAX_PATH];
char TMPFILE[_MAX_PATH];
double RVAL, IVAL=0., TIME, DATA;
double N=0., Imean=0., Rmean=0.;
FILE *INSTR, *OUTSTR;
struct FILEHDR FileHeader;
short ERRFLAG=0, FLAG=0, NUMDAT, NDAT, SMALL=0;
struct RData  *RDataPntr;
struct XData  *XDataPntr;
struct XData  XDataOut;
double Nyquist, Fundamental, HalfN;

main(argc,argv)
short argc;
char *argv[];
   {
   double Nu;
   short P1, P2=0, P3;

   signal(SIGINT,BREAKREQ);

   if (Ztskinit("DFT",argv[0])) Zexit(1);
   if ((CODE < 0) || (CODE > 1))
      {
      Zencode(10,"DFT     : Invalid CODE Specification\n");
      Zexit(1);
      }

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

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

   InitializeDFT();

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

   P3 = (short)Max(10.,100./N);

   RDataPntr =  (struct RData *)BUFFER;
   XDataPntr =  (struct XData *)BUFFER;
   XDataOut.f = 0;

   for (Nu=0.; Nu < N; ++Nu)
      {
      P1 = (short)(100. * Nu/N);
      if (P1>=P2)
         {
         Zencode(3,"DFT     : %3d%% Complete\n",P2);
         P2 += P3;
         }

      XDataOut.z.x = XDataOut.z.y = 0.;
      if (SMALL)
         MEMREDUCE(Nu);      /* Data Fit into Memory */
      else
         if (ERRFLAG = DSKREDUCE(Nu)) break;  /* Disk Processing */

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

   if (!ERRFLAG)
      {
      Zencode(3,"DFT     : 100%% Complete\n");

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

/***************************************************************
**
** Process "small" data files
*/
void MEMREDUCE(Nu)
double Nu;
   {
   short I;
   double Tau;
/*
** Determine the output values
*/
   RDataPntr =  (struct RData *)BUFFER;
   XDataPntr =  (struct XData *)BUFFER;
   for (I=0; I<NUMDAT; ++I)
      {
      switch (FileHeader.type)
         {
         case R_Data:
            RVAL = RDataPntr->y;
            FLAG = RDataPntr->f;
            ++RDataPntr;
            break;
         case X_Data:
            RVAL = XDataPntr->z.x;
            IVAL = XDataPntr->z.y;
            FLAG = XDataPntr->f;
            ++XDataPntr;
            break;
         }
      if (!FLAG)
         {
         Tau = TwoPi * (double)I * Nu / N;
         if (!CODE)
            {
            XDataOut.z.x += (RVAL * cos(Tau) - IVAL * sin(Tau));
            XDataOut.z.y += (IVAL * cos(Tau) + RVAL * sin(Tau));
            }
         else
            {
            XDataOut.z.x += (RVAL * cos(Tau) + IVAL * sin(Tau));
            XDataOut.z.y += (IVAL * cos(Tau) - RVAL * sin(Tau));
            }
         }
      }

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


   return;
   }

/***************************************************************
**
** Process "large" data files
*/
short DSKREDUCE(Nu)
double Nu;
   {
   double Tau, II=0.;
/*
** Determine the output values
*/
   if (!Zgethead(INSTR,(struct FILEHDR *)NULL)) return(1);
   while (Zread(INSTR,BUFFER,FileHeader.type))
      {
      ++II;
      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 = XDataPntr->f;
            break;
         }
      if (!FLAG)
         {
         Tau = TwoPi * II * Nu / N;
         if (!CODE)
            {
            XDataOut.z.x += (RVAL * cos(Tau) - IVAL * sin(Tau));
            XDataOut.z.y += (IVAL * cos(Tau) + RVAL * sin(Tau));
            }
         else
            {
            XDataOut.z.x += (RVAL * cos(Tau) + IVAL * sin(Tau));
            XDataOut.z.y += (IVAL * cos(Tau) - RVAL * sin(Tau));
            }
         }
      }

   if (ferror(INSTR)) return(1);

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

   return(0);
   }

/************************************************************
**
** 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 InitializeDFT()
   {
   short I;
   double GoodN = 0.;

   while (NDAT = Zgetdat(INSTR,BUFFER,FileHeader.type))
      {
      N += (double)(NUMDAT = NDAT);
      RDataPntr =  (struct RData *)BUFFER;
      XDataPntr =  (struct XData *)BUFFER;
      if (SMALL < 2) ++SMALL;
      for (I=0; I<NUMDAT; ++I)          /* Iterate over data buffer */
         {
         switch (FileHeader.type)
            {
            case R_Data:
               RVAL = RDataPntr->y;
               FLAG = RDataPntr->f;
               ++RDataPntr;             /* Advance to next value */
               break;
            case X_Data:
               RVAL = XDataPntr->z.x;
               IVAL = XDataPntr->z.y;
               FLAG = XDataPntr->f;
               ++XDataPntr;             /* Advance to next value */
               break;
            default:
               Zencode(10,"DFT     : Invalid File Type.\n");
               Zexit(1);
            }
         if (!FLAG)                 /* Only use good data */
            {
            Imean += IVAL;
            Rmean += RVAL;
            ++GoodN;
            } /* End IF */
         } /* End For */
      } /* End WHILE */

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

   if (N < 2.)
      {
      Zencode(10,"DFT     : ** ERROR ** File Contains %G Point(s)\n",N);
      Zexit(1);
      }

   Rmean /= GoodN;
   Imean /= GoodN;
/*
** If the data was read in 1 pass then the entire data base can
** be processed in memory (i.e. it is SMALL), otherwise the
** file must be processed from disk, which takes a great deal
** more time.
*/
   SMALL = (SMALL > 1) ? 0 : 1;

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

   Zencode(3,"DFT     : File Contains %G Total Points\n",N);
   Zencode(3,"DFT     : And %G Valid Points\n",GoodN);
   Zencode(3,"DFT     : Mean Value =  %G + i(%G)\n",Rmean,Imean);
   Zencode(3,"DFT     : Fundamental Frequency = %G\n",Fundamental );
   Zencode(3,"DFT     : Nyquist Frequency = %G\n",Nyquist);

   return;
   }

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