/*********************************************************************
*
*  Task PGRAM
*
* Task to Calculate the Periodogram of the Data within the angular
* frequency range TRANGE with FACTOR output points.
* Factor must be greater than 3 and less than or euqal to 2048
*
* CODE = 1 -> Scale by the variance
* CODE = 2 -> Change to % Probability Level
*
* Defaults are 2pi/T to (2pi/min(t))/2 and N points
*
* This routine buffers the data and bypasses the disk I/O
* for small data bases.
*
*/
#include <errno.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <io.h>
#include <time.h>
#include <dos.h>
#include <process.h>
#include <ctype.h>
#include <signal.h>

#include "TISAN.H"

#define Twopi 6.28318530718

double TAU(double);
void PXW(double,double);
void VARIANCE(void);

char IBUFFER[Maxbuff*sizeof(struct TRData)]; /* Largest Data Type */
char INFILE[_MAX_PATH], OUTFILE[_MAX_PATH], TMPFILE[_MAX_PATH];
double N0, SIGMA2, MEAN;
FILE *INSTR, *OUTSTR;
struct FILEHDR FileHeader;
short NUMPNT, SMALL=0, NUMDAT, NDAT;
struct RData  RDATBUFF[Maxbuff];
short  FLAG=0;
double DATA, TIME;

main(argc,argv)
short argc;
char *argv[];
   {
   short I, ERRFLAG=0, FFLAG=0, P1, P2=0, P3;
   double FIRST, OMEGA, SLOPE;
   double OSTART, OSTOP, TCNT=0.;
   double SUM=0., SUM2=0., DELT=0., FT, DT;
   struct RData  *RDataPntr;
   struct TRData *TRDataPntr;

   signal(SIGINT,BREAKREQ);      /* Process ^C Interrupt */

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

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

   Zencode(2,"PGRAM   : Opening Input File '%s'\n",INFILE);
   if ((INSTR = Zopen(INFILE,O_readb)) == NULL) Zexit(1);

   if (!Zgethead(INSTR,&FileHeader)) Zexit(1);
   if ((FileHeader.type != R_Data) &&
       (FileHeader.type != TR_Data))
      {
      Zencode(10,"PGRAM   : Invalid File Type.\n");
      Zexit(1);
      }

   OSTART = TRANGE[0] * Twopi;     /* Omega Start and Stop Ranges */
   OSTOP  = TRANGE[1] * Twopi;
/*
** Find Size of File, MEAN and SIGMA2.
** File is loaded into IBUFFER and NUMDAT is set.
** If the file is small, (i.e. it fit in memory on 1 read)
** then IBUFFER and NUMDAT are not changed ever again.
**
*/
   while (NDAT = Zgetdat(INSTR,IBUFFER,FileHeader.type))
      {
      NUMDAT = NDAT;
      RDataPntr =  (struct RData *)IBUFFER;
      TRDataPntr = (struct TRData *)IBUFFER;
      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;
               DATA = RDataPntr->y;
               FLAG = RDataPntr->f;
               ++RDataPntr;             /* Advance to next value */
               break;
            case TR_Data:
               TIME = TRDataPntr->t;
               DATA = TRDataPntr->y;
               ++TRDataPntr;            /* Advance to next value */
               break;
            }
         if (!FLAG)                 /* Only use good data */
            {
            ++N0;
            if (!FFLAG)
               {
               FIRST = FT = TIME;
               FFLAG = 1;
               }
            else
               {
               if (!DELT) 
                  DELT = TIME - FT;
               else
                  {
                  DT = TIME - FT;
                  if (DT) DELT = Min(DELT,DT);
                  }
               FT = TIME;
               }
            SUM += DATA;
            SUM2 += Square(DATA);
            }
         }
      } /* End WHILE */

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

   SMALL = (SMALL > 1) ? 0 : 1;
   if (N0 < 2.)
      {
      Zencode(10,"PGRAM   : Insufficient Data\n");
      Zexit(1);
      }

   SIGMA2 = ((N0 * SUM2) - Square(SUM))/(N0*(N0-1.));
   MEAN = SUM/N0;
   DELT = 1./(2.*DELT);
   if (OSTART >= OSTOP) OSTART = OSTOP = 0.;
   if (!OSTOP) OSTOP  = Twopi*DELT;
   if (!OSTART) OSTART = Twopi/(TIME-FIRST);

   Zencode(3,"PGRAM   : File Contains %G Valid Points\n",N0);
   Zencode(3,"PGRAM   : Variance = %G\n",SIGMA2);
   Zencode(3,"PGRAM   : Mean = %G\n",MEAN);
   Zencode(3,"PGRAM   : Fundamental Frequency = %G\n",1/(TIME-FIRST));
   Zencode(3,"PGRAM   : Nyquist Frequency = %G\n",DELT);

   if (FACTOR < 3.)           /* Set Defaults for funny FACTOR values */
      NUMPNT = (short)N0;
   else
      NUMPNT = (short)FACTOR;

   NUMPNT = Min(NUMPNT,Maxbuff);

   P3 = Max(10,100./NUMPNT);        /* Used to find percentage complete */
   SLOPE = (OSTOP-OSTART)/(double)(NUMPNT-1);   /* Omega value slope */

   for (I=0;I<NUMPNT;++I)
      {
      P1 = (short)((double)100*(double)I/(double)NUMPNT);
      if (P1>=P2)
         {
         Zencode(3,"PGRAM   : %3d%% Complete\n",P2);
         P2 += P3;
         }
      OMEGA = OSTART + (double)I * SLOPE;
      PXW(TAU(OMEGA),OMEGA);
      } /* End For */

   Zencode(3,"PGRAM   : 100%% Complete\n");

   VARIANCE();

   Zencode(2,"PGRAM   : Opening Scratch File '%s'\n",TMPFILE);
   if ((OUTSTR = Zopen(TMPFILE,O_writeb)) == NULL) Zexit(1);

   FileHeader.b = OSTART/Twopi;     /* Build Output header */
   FileHeader.m = SLOPE/Twopi;
   FileHeader.type = R_Data;
   if (Zputhead(OUTSTR,&FileHeader))
      ERRFLAG=1;
   else
      ERRFLAG = Zputdat(NUMPNT,OUTSTR,(char *)RDATBUFF,R_Data);

   Zclose(INSTR);
   Zclose(OUTSTR);
   ERRFLAG = Znameout(OUTFILE,TMPFILE,ERRFLAG);
   Zexit(ERRFLAG);
   }

/*********************************************************************
*
* Function to Calculate Tau's
*
*/
double TAU(OMEGA)
double OMEGA;
   {
   double ARG, V1=0.,V2=0.,TCNT=0.;
   short J;
   struct RData  *RDataPntr;
   struct TRData *TRDataPntr;

   RDataPntr =  (struct RData *)IBUFFER;
   TRDataPntr = (struct TRData *)IBUFFER;
   OMEGA *= 2.;

   if (SMALL)        /* Data is memory resident */
      {
      for (J=0;J<NUMDAT;++J)
         {
         switch (FileHeader.type)
            {
            case R_Data:
               TIME = TCNT * FileHeader.m + FileHeader.b;
               ++TCNT;
               FLAG = RDataPntr->f;
               break;
            case TR_Data:
               TIME = TRDataPntr->t;
               ++TRDataPntr;
               break;
            }
         if (!FLAG)
            {
            ARG = OMEGA*TIME;
            V1 += sin(ARG);
            V2 += cos(ARG);
            }
         }  /* End For */
      }  /* End IF */
   else                    /* Read data from disk */
      {
      if (!Zgethead(INSTR,(struct FILEHDR *)NULL)) Zexit(1);
      while (Zread(INSTR,IBUFFER,FileHeader.type))
         {
         switch (FileHeader.type)
            {
            case R_Data:
               TIME = TCNT * FileHeader.m + FileHeader.b;
               ++TCNT;
               FLAG = RDataPntr->f;
               break;
            case TR_Data:
               TIME = TRDataPntr->t;
               break;
            }
         if (!FLAG)
            {
            ARG = OMEGA*TIME;
            V1 += sin(ARG);
            V2 += cos(ARG);
            }
         } /* End WHILE */
      if (ferror(INSTR))
         {
         Zerror();
         Zexit(1);
         }
      }    /* End ELSE */
   return(atan2(V1,V2)/OMEGA);
   }

/*********************************************************************
*
* Function to Calculate Periodogram
*
*/
void PXW(TAUV,OMEGA)
double TAUV, OMEGA;
   {
   double V1=0., V2=0., V3=0., V4=0., V, ARG, TCNT=0.;
   short J;
   static short M=0;
   struct RData  *RDataPntr;
   struct TRData *TRDataPntr;

   RDataPntr =  (struct RData *)IBUFFER;
   TRDataPntr = (struct TRData *)IBUFFER;

   if (SMALL)
      {
      for (J=0;J<NUMDAT;++J)
         {
         switch (FileHeader.type)
            {
            case R_Data:
               TIME = TCNT * FileHeader.m + FileHeader.b;
               DATA = RDataPntr->y - MEAN;
               FLAG = RDataPntr->f;
               ++RDataPntr;
               ++TCNT;
               break;
            case TR_Data:
               TIME = TRDataPntr->t;
               DATA = TRDataPntr->y - MEAN;
               ++TRDataPntr;
               break;
            }
         if (!FLAG)
            {
            ARG = OMEGA * (TIME - TAUV);
            V = cos(ARG);
            V1 += DATA * V;
            V3 += Square(V);
            V = sin(ARG);
            V2 += DATA * V;
            V4 += Square(V);
            }
         }  /* End For */
      } /* End IF */
   else
      {
      if (!Zgethead(INSTR,(struct FILEHDR *)NULL)) Zexit(1);
      while (Zread(INSTR,IBUFFER,FileHeader.type))
         {
         switch (FileHeader.type)
            {
            case R_Data:
               TIME = TCNT * FileHeader.m + FileHeader.b;
               DATA = RDataPntr->y - MEAN;
               FLAG = RDataPntr->f;
               ++TCNT;
               break;
            case TR_Data:
               TIME = TRDataPntr->t;
               DATA = TRDataPntr->y - MEAN;
               break;
            }
         if (!FLAG)
            {
            ARG = OMEGA * (TIME - TAUV);
            V = cos(ARG);
            V1 += DATA * V;
            V3 += Square(V);
            V = sin(ARG);
            V2 += DATA * V;
            V4 += Square(V);
            }
         } /* End While */
      if (ferror(INSTR))
         {
         Zerror();
         Zexit(1);
         }
      } /* End ELSE */

   RDATBUFF[M].y = (V1*(V1/V3) + V2*(V2/V4))/2.;
   ++M;
   return;
   }

/*********************************************************************
*
* Function to Scale Pgram by the Variance if CODE=1
* or as probability if CODE=2
*
*/
void VARIANCE()
   {
   short I;
   double NI, M;

   M = NI = -6.362 + 1.193*N0 + .00098*N0*N0;

   for (I=0;I<NUMPNT;++I)
      {
      if (CODE) RDATBUFF[I].y /= SIGMA2;

      if (CODE == 2)
         RDATBUFF[I].y = pow(1. - exp(-RDATBUFF[I].y),NI);

      RDATBUFF[I].f = 0;
      }

   if (CODE)
      {
      DATA = 0.;
      while (M > 1.)
         {
         DATA += (1./M);
         --M;
         }
      if (CODE == 2) DATA = pow(1. - exp(-DATA),NI);
      Zencode(3,"PGRAM   : Expected Noise Peak = %G\n",DATA);
      }

   return;
   }

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