/*********************************************************************
*
*  Task FIT
*
* Fourier Integral Transform
* CODE = 0 ->  FT
* CODE = 1 -> IFT
*
*/
#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, HalfN;
double N=0., a1H0H1ON;
double TOTAL=0.;
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("FIT",argv[0])) Zexit(1);     /* Initialize Task */

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

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

   Zencode(2,"FIT     : 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:
      case R_Data:
      case TR_Data:
         break;
      default:
         Zencode(10,"FIT     : Invalid File Type.\n");
         Zexit(1);
      }

   FINIT();          /* Do some initialization */

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

   if (FACTOR<3.) FACTOR = (double)N;
   FACTOR = floor(FACTOR);
   SLOPE = (TRANGE[1]-TRANGE[0])/(FACTOR-1.);
   
   XDataOut.f = 0;

   P3 = 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)
      {
      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,"FIT     : %3d%% Complete\n",P2);
            P2 +=P3;
            }
         }

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

/*
      if (!CODE)
         {
         XDataOut.z.y /= HalfN;
         XDataOut.z.x /= HalfN;
         }
      else
         {
         XDataOut.z.y /= 2.;
         XDataOut.z.x /= 2.;
         }
*/
      Zwrite(OUTSTR,(char *)&XDataOut,X_Data);
      } /* End FOR II */


   Zencode(3,"FIT     : 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,"FIT     : ** 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;
   HalfN = N / 2.;
   Zencode(3,"FIT     : File Contains %G Points\n",N);
   Zencode(3,"FIT     : Mean Value = %G + i(%G)\n",RMEAN, IMEAN);

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

   return;
   }

/***************************************************************
**
** Process "small" data files
*/
void MEMREDUCE()
   {
   short I, Pass1 = 1;
   double h1, h2, LastTime, LastX, LastY, DeltaT, ThisX, ThisY;

   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;
         }
      if (Pass1)
         {
         Pass1 = 0;
         LastTime = TIME;
         TIME *= Omega;
         h1 = cos(TIME);
         h2 = sin(TIME);
         switch (CODE)
            {
            case 0: /* FIT */
               LastX = ((RDATA * h1) + (IDATA * h2));
               LastY = ((IDATA * h1) - (RDATA * h2));
               break;
            case 1: /* IFIT */
               LastX = ((RDATA * h1) - (IDATA * h2));
               LastY = ((IDATA * h1) + (RDATA * h2));
               break;
            }
         }
      else
         {
         DeltaT = TIME - LastTime;
         LastTime = TIME;
         TIME *= Omega;
         h1 = cos(TIME);
         h2 = sin(TIME);
         switch (CODE)
            {
            case 0: /* FIT */
               ThisX = (RDATA * h1) + (IDATA * h2);
               ThisY = (IDATA * h1) - (RDATA * h2);
               XDataOut.z.x +=
                  (DeltaT * (ThisX + LastX) / 2.);
               XDataOut.z.y +=
                  (DeltaT * (ThisY + LastY) / 2.);
               break;
            case 1: /* IFIT */
               ThisX = (RDATA * h1) - (IDATA * h2);
               ThisY = (IDATA * h1) + (RDATA * h2);
               XDataOut.z.x +=
                  (DeltaT * (ThisX + LastX) / 2.);
               XDataOut.z.y +=
                  (DeltaT * (ThisY + LastY) / 2.);
               break;
            }
         LastX = ThisX;
         LastY = ThisY;
         }
      }

   return;
   }

/***************************************************************
**
** Process "large" data files
*/
short DSKREDUCE()
   {
   short Pass1 = 1;
   double h1, h2, LastTime, LastX, LastY, DeltaT, ThisX, ThisY;

   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;
         }
      if (Pass1)
         {
         Pass1 = 0;
         LastTime = TIME;
         h1 = cos(TIME);
         h2 = sin(TIME);
         switch (CODE)
            {
            case 0: /* FIT */
               LastX = ((RDATA * h1) + (IDATA * h2));
               LastY = ((IDATA * h1) - (RDATA * h2));
               break;
            case 1: /* IFIT */
               LastX = ((RDATA * h1) - (IDATA * h2));
               LastY = ((IDATA * h1) + (RDATA * h2));
               break;
            }
         }
      else
         {
         DeltaT = TIME - LastTime;
         LastTime = TIME;
         TIME *= Omega;
         h1 = cos(TIME);
         h2 = sin(TIME);
         switch (CODE)
            {
            case 0: /* FIT */
               ThisX = (RDATA * h1) + (IDATA * h2);
               ThisY = (IDATA * h1) - (RDATA * h2);
               XDataOut.z.x +=
                  (DeltaT * (ThisX + LastX) / 2.);
               XDataOut.z.y +=
                  (DeltaT * (ThisY + LastY) / 2.);
               break;
            case 1: /* IFIT */
               ThisX = (RDATA * h1) - (IDATA * h2);
               ThisY = (IDATA * h1) + (RDATA * h2);
               XDataOut.z.x +=
                  (DeltaT * (ThisX + LastX) / 2.);
               XDataOut.z.y +=
                  (DeltaT * (ThisY + LastY) / 2.);
               break;
            }
         LastX = ThisX;
         LastY = ThisY;
         }
      } /* End WHILE */

   return;
   }

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