/*********************************************************************
*
*  Task DCDFT
*
* Date Compenstaed Discrete Fourier Transform
* CODE = 0 ->  DCDFT
* CODE = 1 -> IDCDFT
* CODE = 2 -> FILTER
* CODE = 3 -> FILTER information only
*
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <signal.h>

#include "TISAN.H"

#define TwoPi 6.28318530718

void FINIT(void);
void MEMREDUCE(void);
short DSKREDUCE(void);

FILE *INSTR, *OUTSTR;
double Omega, Nu, SLOPE, II, TIME, RDATA, IDATA=0., TCNT;
double N=0., a1H0H1ON;
double H0H1, H0H2, H1H1, H2H2, H1H2, h1H2, a1, a2, h1, h2;
double d1, d2;
double NSQ, H0H1SQ, H0H2SQ, H1, H2, TOTAL=0, SQNO2, SQNO8;
short  ERRFLAG=0, NUMDAT, SMALL=0, P1, P2=0, P3, FLAG=0, NDAT;
char INFILE[_MAX_PATH], OUTFILE[_MAX_PATH], TMPFILE[_MAX_PATH];

struct FILEHDR FileHeader;
struct RData  *RDataPntr;
struct TRData *TRDataPntr;
struct XData  *XDataPntr;
struct TXData *TXDataPntr;
struct XData  XDataOut;
char BUFFER[Maxbuff*sizeof(struct TXData)];

main(argc,argv)
short argc;
char *argv[];
   {
   signal(SIGINT,BREAKREQ);         /* Handle ^C Interrupt */

   if (Ztskinit("DCDFT",argv[0])) Zexit(1);     /* Initialize Task */

   if ((CODE<0) || (CODE>3))
      {
      Zencode(10,"DCDFT   : Invalid Code Specification\n");
      Zexit(1);
      }

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

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

   switch (FileHeader.type)
      {
      case X_Data:
      case TX_Data:
         if (CODE > 1)
            {
            Zencode(10,"DCDFT   : Can't Filter Complex Files.\n");
            Zexit(1);
            }
      case R_Data:
      case TR_Data:
         break;
      default:
         Zencode(10,"DCDFT   : Invalid File Type.\n");
         Zexit(1);
      }

   FINIT();          /* Do some initialization */

   if (CODE != 3)    /* CODE=3 does not create a file, all others do */
      {
      Zencode(2,"DCDFT   : Opening Scratch File '%s'\n",TMPFILE);
      if ((OUTSTR = Zopen(TMPFILE,O_writeb)) == NULL) Zexit(1);
      if (Zputhead(OUTSTR,&FileHeader))
         {
         ERRFLAG=1;
         goto EXIT;
         }
      }

   if (CODE > 1)
      {
      TRANGE[0] = Nu;
      SLOPE=0.;
      FACTOR=1.;
      }
   else
      {
      if (FACTOR<3.) FACTOR = (double)N;
      FACTOR = floor(FACTOR);
      SLOPE = (TRANGE[1]-TRANGE[0])/(FACTOR-1.);
      }
   
   XDataOut.f = 0;

   P3 = (short)Max(10.,100./FACTOR);
   RDataPntr =  (struct RData *)BUFFER;
   TRDataPntr = (struct TRData *)BUFFER;
   XDataPntr =  (struct XData *)BUFFER;
   TXDataPntr = (struct TXData *)BUFFER;
   for (II=0;II<FACTOR;++II)
      {
      H0H1=H0H2=H1H1=H2H2=H1H2=h1H2=0.;   /* Initialize Iterative Values */
      XDataOut.z.x = XDataOut.z.y = 0.;
      Omega = TwoPi*(II*SLOPE+TRANGE[0]);  /* Calculate next Omega */
      if (CODE < 2)                        /* Display % COMPLETE message */
         {
         P1 = 100.*II/FACTOR;
         if (P1>=P2)
            {
            Zencode(3,"DCDFT   : %3d%% Complete\n",P2);
            P2 +=P3;
            }
         }

      if (SMALL)
         MEMREDUCE();      /* Data Fit into Memory */
      else
         if (ERRFLAG = DSKREDUCE()) goto EXIT;  /* Disk Processing */

      if (CODE > 1)   /* Either Filter or just display results */
         {
         d2 = a2*XDataOut.z.y;
         d1 = a1*XDataOut.z.x + d2*Square(a1)*(H0H1*H0H2/N - H1H2);
         Zencode(5,"DCDFT   : Amplitude = %G\n",sqrt(d1*d1+d2*d2));
         if (d1 || d2) 
            Zencode(5,"DCDFT   : Phase = %G\n",atan2(d2,d1));
         else
            Zencode(5,"DCDFT   : Phase is undefined.\n");
         }
      else
         {
         if (!CODE)
            {
            XDataOut.z.y /= SQNO2;
            XDataOut.z.x /= SQNO2;
            }
         else
            {
            XDataOut.z.y *= SQNO8;
            XDataOut.z.x *= SQNO8;
            }
         Zwrite(OUTSTR,(char *)&XDataOut,X_Data);
         }
      } /* End FOR II */


   if (CODE < 2) /* File is Complex DCDFT so update header information */
      {
      Zencode(3,"DCDFT   : 100%% Complete\n");
      FileHeader.type = X_Data;        /* Output File is Complex */
      FileHeader.m = SLOPE;
      FileHeader.b = TRANGE[0];
      Zputhead(OUTSTR,&FileHeader);
      }

EXIT:
   Zclose(INSTR);
   if (CODE != 3)
      {
      Zclose(OUTSTR);
      ERRFLAG = Znameout(OUTFILE,TMPFILE,ERRFLAG);
      }
   Zexit(ERRFLAG);
   }

/************************************************************
**
** Initialize file stuff
*/
void FINIT()
   {
   double FIRST, FT, DELT=0., RMEAN=0., IMEAN=0.;
   short I, FFLAG=0;

   TCNT = 0.;

   while (NDAT = Zgetdat(INSTR,BUFFER,FileHeader.type))
      {
      NUMDAT = NDAT;
      RDataPntr =  (struct RData *)BUFFER;       /* Setup some data pointers */
      TRDataPntr = (struct TRData *)BUFFER;
      XDataPntr =  (struct XData *)BUFFER;
      TXDataPntr = (struct TXData *)BUFFER;
      if (SMALL < 2) ++SMALL;
      for (I=0;I<NUMDAT;++I)           /* Iterate over data buffer */
         {
         switch (FileHeader.type)
            {
            case R_Data:
               TIME = TCNT * FileHeader.m + FileHeader.b;
               ++TCNT;
               RDATA = RDataPntr->y;
               FLAG = RDataPntr->f;
               ++RDataPntr;             /* Advance to next vale */
               break;
            case TR_Data:
               TIME = TRDataPntr->t;
               RDATA = TRDataPntr->y;
               ++TRDataPntr;            /* Advance to next value */
               break;
            case X_Data:
               TIME = TCNT * FileHeader.m + FileHeader.b;
               ++TCNT;
               RDATA = XDataPntr->z.x;
               IDATA = XDataPntr->z.y;
               FLAG = XDataPntr->f;
               ++XDataPntr;             /* Advance to next vale */
               break;
            case TX_Data:
               TIME = TXDataPntr->t;
               RDATA = TXDataPntr->z.x;
               IDATA = TXDataPntr->z.y;
               ++TXDataPntr;            /* Advance to next value */
               break;
            }
         if (!FLAG)                 /* Only use good data */
            {
            ++N;
            if (!FFLAG)
               {
               FIRST = FT = TIME;
               FFLAG = 1;
               }
            else
               {
               if (!DELT) 
                  DELT = TIME - FT;
               else
                  DELT = Min(DELT,TIME - FT);
               FT = TIME;
               }
            RMEAN += RDATA;
            IMEAN += IDATA;
            } /* End IF */
         } /* End For */
      } /* End WHILE */

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

   if (N<2)
      {
      Zencode(10,"DCDFT   : ** ERROR ** File Contains %G Point(s)\n",N);
      Zexit(1);
      }
/*
** 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;

   RMEAN /= N;
   IMEAN /= N;
   SQNO2 = sqrt(N/2.);
   SQNO8 = sqrt(N/8.);
   Zencode(3,"DCDFT   : File Contains %G Points\n",N);
   Zencode(3,"DCDFT   : Mean Value = %G + i(%G)\n",RMEAN, IMEAN);

   if (CODE > 1) Nu = TRANGE[0];

   if (TRANGE[0] >= TRANGE[1])
      {
      TRANGE[0] = 1./(TIME-FIRST);
      TRANGE[1] = 1./DELT;
      }
   Zencode(3,"DCDFT   : Fundamental Frequency = %G\n",1./(TIME-FIRST));
   Zencode(3,"DCDFT   : Nyquist Frequency = %G\n",1./(2. * DELT));

   return;
   }

/***************************************************************
**
** Process "small" data files
*/
void MEMREDUCE()
   {
   short I;
/*
** First we must determine the inner product sets
** H0H1, H0H2, a1 and a2
*/
   TRDataPntr = (struct TRData *)BUFFER;
   TXDataPntr = (struct TXData *)BUFFER;
   TCNT = 0.;
   for (I=0;I<NUMDAT;++I)
      {
      switch (FileHeader.type)
         {
         case R_Data:
            TIME = TCNT * FileHeader.m + FileHeader.b;
            ++TCNT;
            break;
         case TR_Data:
            TIME = TRDataPntr->t;
            ++TRDataPntr;
            break;
         case X_Data:
            TIME = TCNT * FileHeader.m + FileHeader.b;
            ++TCNT;
            break;
         case TX_Data:
            TIME = TXDataPntr->t;
            ++TXDataPntr;
            break;
         }
      TIME *= Omega;
      H1 = cos(TIME);
      H2 = sin(TIME);

      H0H1 += H1;
      H0H2 += H2;
      H1H1 += (H1*H1);
      H2H2 += (H2*H2);
      H1H2 += (H1*H2);
      }

   H0H1SQ = Square(H0H1);
   H0H2SQ = Square(H0H2);

   a1 = H1H1 - H0H1SQ/N;
   if (a1 > 0.)
      a1 = 1./sqrt(a1);
   else
      a1 = 0.;

   a2 = H2H2 - H0H2SQ/N - Square(a1)*(Square(H1H2) +
               H0H1SQ*H0H2SQ/(N*N) - 2*H0H1*H0H2*H1H2/N);

   if (a2 > 0.)
      a2 = 1./sqrt(a2);
   else
      a2 = 0.;

   a1H0H1ON = a1 * H0H1/N;
/*
** Now we must determine the inner product set h1H2
*/
   TRDataPntr = (struct TRData *)BUFFER;
   TXDataPntr = (struct TXData *)BUFFER;
   TCNT = 0.;
   for (I=0;I<NUMDAT;++I)
      {
      switch (FileHeader.type)
         {
         case R_Data:
            TIME = TCNT * FileHeader.m + FileHeader.b;
            ++TCNT;
            break;
         case TR_Data:
            TIME = TRDataPntr->t;
            ++TRDataPntr;
            break;
         case X_Data:
            TIME = TCNT * FileHeader.m + FileHeader.b;
            ++TCNT;
            break;
         case TX_Data:
            TIME = TXDataPntr->t;
            ++TXDataPntr;
            break;
         }
      TIME *= Omega;
      h1H2 += (a1 * cos(TIME) - a1H0H1ON) * sin(TIME);
      }
/*
** Now determine the output values
*/
   RDataPntr =  (struct RData *)BUFFER;
   TRDataPntr = (struct TRData *)BUFFER;
   XDataPntr =  (struct XData *)BUFFER;
   TXDataPntr = (struct TXData *)BUFFER;
   TCNT = 0.;
   for (I=0;I<NUMDAT;++I)
      {
      switch (FileHeader.type)
         {
         case R_Data:
            TIME = TCNT * FileHeader.m + FileHeader.b;
            RDATA = RDataPntr->y;
            ++TCNT;
            ++RDataPntr;
            break;
         case TR_Data:
            TIME = TRDataPntr->t;
            RDATA = TRDataPntr->y;
            ++TRDataPntr;
            break;
         case X_Data:
            TIME = TCNT * FileHeader.m + FileHeader.b;
            RDATA = XDataPntr->z.x;
            IDATA = XDataPntr->z.y;
            ++TCNT;
            ++XDataPntr;
            break;
         case TX_Data:
            TIME = TXDataPntr->t;
            RDATA = TXDataPntr->z.x;
            IDATA = TXDataPntr->z.y;
            ++TXDataPntr;
            break;
         }
      TIME *= Omega;
      h1 = a1 * (cos(TIME) - H0H1/N);
      h2 = a2*(sin(TIME) - H0H2/N - h1*h1H2);
      switch (CODE)
         {
         case 0: /* DCDFT */
            XDataOut.z.x += ((RDATA * h1) + (IDATA * h2));
            XDataOut.z.y += ((IDATA * h1) - (RDATA * h2));
            break;
         case 1: /* IDCDFT */
            XDataOut.z.x += ((RDATA * h1) - (IDATA * h2));
            XDataOut.z.y += ((IDATA * h1) + (RDATA * h2));
            break;
         default: /* Filter */
            XDataOut.z.x += (RDATA * h1);
            XDataOut.z.y += (RDATA * h2);
         }
      }

   if (CODE==2)    /* Filter Data */
      {
      RDataPntr =  (struct RData *)BUFFER;
      TRDataPntr = (struct TRData *)BUFFER;
      TCNT = 0.;
      for (I=0;I<NUMDAT;++I)
         {
         switch (FileHeader.type)
            {
            case R_Data:
               TIME = TCNT * FileHeader.m + FileHeader.b;
               ++TCNT;
               break;
            case TR_Data:
               TIME = TRDataPntr->t;
               break;
            }
         TIME *= Omega;
         h1 = a1 * (cos(TIME) - H0H1/N);
         h2 = a2*(sin(TIME) - H0H2/N - h1*h1H2);
         switch (FileHeader.type)
            {
            case R_Data:
               RDataPntr->y -= (XDataOut.z.x*h1 + XDataOut.z.y*h2);
               ++RDataPntr;
               break;
            case TR_Data:
               TRDataPntr->y -= (XDataOut.z.x*h1 + XDataOut.z.y*h2);
               ++TRDataPntr;
               break;
            }
         } /* END for I */
      Zputdat(NUMDAT,OUTSTR,(char *)BUFFER,FileHeader.type);
      }
   return;
   }

/***************************************************************
**
** Process "large" data files
*/
short DSKREDUCE()
   {
/*
** First we must determine the inner product sets
** H0H1, H0H2, a1 and a2
*/
   TCNT = 0.;
   if (!Zgethead(INSTR,(struct FILEHDR *)NULL)) return(1);
   while (Zread(INSTR,BUFFER,FileHeader.type))
      {
      switch (FileHeader.type)
         {
         case R_Data:
            TIME = TCNT * FileHeader.m + FileHeader.b;
            ++TCNT;
            break;
         case TR_Data:
            TIME = TRDataPntr->t;
            break;
         case X_Data:
            TIME = TCNT * FileHeader.m + FileHeader.b;
            ++TCNT;
            break;
         case TX_Data:
            TIME = TXDataPntr->t;
            break;
         }
      TIME *= Omega;
      H1 = cos(TIME);
      H2 = sin(TIME);
      H0H1 += H1;
      H0H2 += H2;
      H1H1 += (H1*H1);
      H2H2 += (H2*H2);
      H1H2 += (H1*H2);
      }

   H0H1SQ = Square(H0H1);
   H0H2SQ = Square(H0H2);

   a1 = H1H1 - H0H1SQ/N;
   if (a1 > 0.)
      a1 = 1./sqrt(a1);
   else
      a1 = 0.;
   a2 = H2H2 - H0H2SQ/N - Square(a1)*(Square(H1H2) +
               H0H1SQ*H0H2SQ/(N*N) - 2*H0H1*H0H2*H1H2/N);
   if (a2 > 0.)
      a2 = 1./sqrt(a2);
   else
      a2 = 0.;

   a1H0H1ON = a1 * H0H1/N;
/*
** Now we must determine the inner product set h1H2
*/
   TCNT = 0.;
   if (!Zgethead(INSTR,(struct FILEHDR *)NULL)) return(1);
   while (Zread(INSTR,BUFFER,FileHeader.type))
      {
      switch (FileHeader.type)
         {
         case R_Data:
            TIME = TCNT * FileHeader.m + FileHeader.b;
            ++TCNT;
            break;
         case TR_Data:
            TIME = TRDataPntr->t;
            break;
         case X_Data:
            TIME = TCNT * FileHeader.m + FileHeader.b;
            ++TCNT;
            break;
         case TX_Data:
            TIME = TXDataPntr->t;
            break;
         }
      TIME *= Omega;
      h1H2 += (a1 * cos(TIME) - a1H0H1ON) * sin(TIME);
      }
/*
** Now determine the output values
*/
   TCNT = 0.;
   if (!Zgethead(INSTR,(struct FILEHDR *)NULL)) return(1);
   while (Zread(INSTR,BUFFER,FileHeader.type))
      {
      switch (FileHeader.type)
         {
         case R_Data:
            TIME = TCNT * FileHeader.m + FileHeader.b;
            RDATA = RDataPntr->y;
            ++TCNT;
            break;
         case TR_Data:
            TIME = TRDataPntr->t;
            RDATA = TRDataPntr->y;
            break;
         case X_Data:
            TIME = TCNT * FileHeader.m + FileHeader.b;
            RDATA = XDataPntr->z.x;
            IDATA = XDataPntr->z.y;
            ++TCNT;
            break;
         case TX_Data:
            TIME = TXDataPntr->t;
            RDATA = TXDataPntr->z.x;
            IDATA = TXDataPntr->z.y;
            break;
         }
      TIME *= Omega;
      h1 = a1 * (cos(TIME) - H0H1/N);
      h2 = a2*(sin(TIME) - H0H2/N - h1*h1H2);
      switch (CODE)
         {
         case 0: /* DCDFT */
            XDataOut.z.x += ((RDATA * h1) + (IDATA * h2));
            XDataOut.z.y += ((IDATA * h1) - (RDATA * h2));
            break;
         case 1: /* IDCDFT */
            XDataOut.z.x += ((RDATA * h1) - (IDATA * h2));
            XDataOut.z.y += ((IDATA * h1) + (RDATA * h2));
            break;
         default: /* Filter */
            XDataOut.z.x += (RDATA * h1);
            XDataOut.z.y += (RDATA * h2);
         }
      } /* End WHILE */

   if (CODE==2)    /* Filter Data */
      {
      if (!Zgethead(INSTR,(struct FILEHDR *)NULL)) return(1);
      TCNT = 0.;
      while (Zread(INSTR,BUFFER,FileHeader.type))
         {
         switch (FileHeader.type)
            {
            case R_Data:
               TIME = TCNT * FileHeader.m + FileHeader.b;
               ++TCNT;
               break;
            case TR_Data:
               TIME = TRDataPntr->t;
               break;
            }
         TIME *= Omega;
         h1 = a1 * (cos(TIME) - H0H1/N);
         h2 = a2*(sin(TIME) - H0H2/N - h1*h1H2);
         switch (FileHeader.type)
            {
            case R_Data:
               RDataPntr->y -= (XDataOut.z.x*h1 + XDataOut.z.y*h2);
               break;
            case TR_Data:
               TRDataPntr->y -= (XDataOut.z.x*h1 + XDataOut.z.y*h2);
               break;
            }
         Zwrite(OUTSTR,(char *)BUFFER,FileHeader.type);
         } /* END while */
      }
   return;
   }

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