// To write out single pulses after correcting for baseline
// To Compile: gcc outstokes.c  nrutil.c -o outstokes -lm 
//

#include <stdio.h>
#include <string.h>
#include <math.h>
#include <stdlib.h>
#include "snglfunc.h"


#define USAGE "\
usage: outstokes dedispfile filelist filestat <parfile>\n\
filelist:\n\
ifile\n\
qfile\n\
ufile\n\
vfile\n\
inp_hdr setup: \n\
bwin -1\n\
ewin -1\n\
nbinavg 1\n\
npuls_avg 1\n"


int main(int argc, char *argv[]) 
{
   FILE                 *fpt, *fpta;
   char                 *outstokesfile, *outpolfile, *outstokesascii; 
   char                 *outhdrfile, *outfold, *parfile;
   char                 *ifile, *qfile, *ufile, *vfile;
   char                 *infile, *statfile, *dedispfile, *inphdr[20];
   char                 *readline, *readval[2];
   float                *dedispi, *dedispq, *dedispu, *dedispv, *dwrt;
   float                *sftdispi, *sftdispq, *sftdispu, *sftdispv;
   double                period, int_time;
   double               *pulstki, *pulstkq, *pulstku, *pulstkv; 
   double               *foldstki, *foldstkq, *foldstku, *foldstkv; 
   double               *pulstkl, *stkppa, *ppaerr; 
   double               *foldarray, *inarrq, *inarru;
   double               *runmean, *statsngl;
   double                minstat[2], rmsu, rmsq;
   int                   nsize, nskip, nbin, npuls, nbinavg, nrun, rbin;
   long long int         fsize;
   int                   win[6], *nullprd;
   int                   ii, jj, sprd, bbase, ebase, npuls_avg, indx;


  
   dedispfile = (char *)calloc(WORDSIZE, sizeof(char)); 
   infile = (char *)calloc(WORDSIZE, sizeof(char)); 
   ifile = (char *)calloc(WORDSIZE, sizeof(char));
   qfile = (char *)calloc(WORDSIZE, sizeof(char));
   ufile = (char *)calloc(WORDSIZE, sizeof(char));
   vfile = (char *)calloc(WORDSIZE, sizeof(char));
   statfile = (char *)calloc(WORDSIZE, sizeof(char));
   outstokesfile = (char *)calloc(WORDSIZE, sizeof(char));
   outstokesascii = (char *)calloc(WORDSIZE, sizeof(char));
   outfold = (char *)calloc(WORDSIZE, sizeof(char));
   outpolfile = (char *)calloc(WORDSIZE, sizeof(char));
   outhdrfile = (char *)calloc(WORDSIZE, sizeof(char));
   parfile = (char *)calloc(WORDSIZE, sizeof(char));

   if (argc < 4)
      cmdlineError (USAGE);

   ii = 1;
   sprintf (dedispfile,"%s", argv[ii++]);
   sprintf (infile,"%s", argv[ii++]);
   sprintf (statfile,"%s", argv[ii++]);

   win[0] = -1;
   win[1] = -1;
   nbinavg = 1;
   npuls_avg = 1;
   if (argc == 5) {
      sprintf (parfile,"%s", argv[ii++]);

      readheader(parfile, inphdr);
      win[0] = atoi(inphdr[0]);
      win[1] = atoi(inphdr[1]);
      nbinavg = atoi(inphdr[2]);
      npuls_avg = atoi(inphdr[3]);
   }

   fprintf (stderr,"okies1\n");

   // Reading the the filenames from infile
   readline = (char *)calloc(LINELEN, sizeof(char));

   if ((fpt = fopen(infile,"r")) == NULL){
      fprintf(stderr, "Error in opening %s file\n", infile);
      exit(1);
   }

   fgets(readline, LINELEN, fpt); 
   splitstr (readline, readval, " ", 1);
   sprintf (ifile,"%s", readval[0]);

   fgets(readline, LINELEN, fpt); 
   splitstr (readline, readval, " ", 1);
   sprintf (qfile,"%s", readval[0]);

   fgets(readline, LINELEN, fpt); 
   splitstr (readline, readval, " ", 1);
   sprintf (ufile,"%s", readval[0]);

   fgets(readline, LINELEN, fpt); 
   splitstr (readline, readval, " ", 1);
   sprintf (vfile,"%s", readval[0]);

   fclose(fpt);


   // Reading header from header file
   readheader(statfile, inphdr);

   period = atof(inphdr[0]);
   int_time = atof(inphdr[8]); 

   fprintf(stderr,"period %lf  int_time %lf\n", period, int_time);

   nbin = floor(period/(int_time*(double)(nbinavg))) + ZERO_lf;
 
   // Using I-Stokes file to read size.

   // Determine file size
   if ((fpt = fopen(ifile,"rb")) == NULL){
      fprintf(stderr, "Error in opening %s file\n", ifile);
      exit(1);
   }

   fseeko(fpt, 0, SEEK_END);
   fsize = ftello(fpt);
   fseeko(fpt, 0, SEEK_SET);

   nsize = fsize/sizeof(float);

   // Allocating memory 
   dwrt    = (float *)calloc((1), sizeof(float));
   dedispi = (float *)calloc(nsize, sizeof(float));
   dedispq = (float *)calloc(nsize, sizeof(float));
   dedispu = (float *)calloc(nsize, sizeof(float));
   dedispv = (float *)calloc(nsize, sizeof(float));

   fprintf (stderr,"Reading file %s\n", ifile);

   // Reading I-Stokes file
   ii = 0;
   while (!feof(fpt)) {
      fread(dwrt, sizeof(float), 1, fpt);
      dedispi[ii++] = dwrt[0];
   }
   fclose(fpt);

   // Reading Q-Stokes file
   if ((fpt = fopen(qfile,"rb")) == NULL){
      fprintf(stderr, "Error in opening %s file\n", qfile);
      exit(1);
   }

   fprintf (stderr,"Reading file %s\n", qfile);

   ii = 0;
   while (!feof(fpt)) {
      fread(dwrt, sizeof(float), 1, fpt);
      dedispq[ii++] = dwrt[0];
   }
   fclose(fpt);

   // Reading U-Stokes file
   if ((fpt = fopen(ufile,"rb")) == NULL){
      fprintf(stderr, "Error in opening %s file\n", ufile);
      exit(1);
   }

   fprintf (stderr,"Reading file %s\n", ufile);

   ii = 0;
   while (!feof(fpt)) {
      fread(dwrt, sizeof(float), 1, fpt);
      dedispu[ii++] = dwrt[0];
   }
   fclose(fpt);

   // Reading V-Stokes file
   if ((fpt = fopen(vfile,"rb")) == NULL){
      fprintf(stderr, "Error in opening %s file\n", vfile);
      exit(1);
   }

   fprintf (stderr,"Reading file %s\n", vfile);

   ii = 0;
   while (!feof(fpt)) {
      fread(dwrt, sizeof(float), 1, fpt);
      dedispv[ii++] = dwrt[0];
   }
   fclose(fpt);

   fprintf(stderr,"nsize %d  %d \n", nsize, ii);


   // Folding dedisp file
   foldstki = (double *)calloc(nbin+1, sizeof(double));
   foldstkq = (double *)calloc(nbin+1, sizeof(double));
   foldstku = (double *)calloc(nbin+1, sizeof(double));
   foldstkv = (double *)calloc(nbin+1, sizeof(double));

   npuls = floor((double)(nsize*int_time/period)) + ZERO_lf;

   foldarray = (double *)calloc(nbin, sizeof(double));

   dedispfold (dedispi, foldarray, nbin, nsize, period, int_time, 0, npuls);
/* 
   dedispfold (dedispi, foldstki, nbin, nsize, period, int_time, 0, npuls);
 
   dedispfold (dedispq, foldstkq, nbin, nsize, period, int_time, 0, npuls);
 
   dedispfold (dedispu, foldstku, nbin, nsize, period, int_time, 0, npuls);
 
   dedispfold (dedispv, foldstkv, nbin, nsize, period, int_time, 0, npuls);
*/ 

   // finding Peak and shifting dedisp to move peak to center.
   sftdispi = (float *)calloc(nsize, sizeof(float));
   sftdispq = (float *)calloc(nsize, sizeof(float));
   sftdispu = (float *)calloc(nsize, sizeof(float));
   sftdispv = (float *)calloc(nsize, sizeof(float));

   nskip = sftpuls (dedispi, sftdispi, foldarray, nbin, nsize, period, 
                      int_time);
  
   // shifting the Q,U,V-Stokes to peak of I-Stokes.
   sftdedisp (dedispq, sftdispq, nsize, nskip);

   sftdedisp (dedispu, sftdispu, nsize, nskip);

   sftdedisp (dedispv, sftdispv, nsize, nskip);

   
   // Folding shifted dedisp file
   dedispfold (sftdispi, foldarray, nbin, nsize-nskip, period, int_time, 0, 
                 npuls);


   // finding base with minimum rms of pulse   
   findbase(foldarray, nbin, 5.0, &bbase, &ebase, minstat); 

   fprintf(stderr,"bbase %d ebase %d nskip %d mean minbase %lf rms minbase %lf\n", bbase, ebase, nskip, minstat[0], minstat[1]);

   
   // finding width of pulse
   if (win[0] == -1 && win[1] == -1) {
      pulswid(foldarray, nbin, minstat, win, 3, 3);

      if (win[0] == 0 || win[1] == 0) {
         win[0] = nbin/2-nbin/20;
         win[1] = nbin/2+nbin/20;   
      }
   }
   fprintf(stderr,"bwin %d ewin %d \n", win[0], win[1]);

   
   // Single pulse stack for I,Q,U,V-Stokes
   npuls = floor((double)(nsize-nskip)*int_time/period) + ZERO_lf;

   pulstki  = (double *)calloc(nbin*npuls+1, sizeof(double));
   pulstkq  = (double *)calloc(nbin*npuls+1, sizeof(double));
   pulstku  = (double *)calloc(nbin*npuls+1, sizeof(double));
   pulstkv  = (double *)calloc(nbin*npuls+1, sizeof(double));
   pulstkl  = (double *)calloc(nbin*npuls+1, sizeof(double));
   stkppa   = (double *)calloc(nbin*npuls+1, sizeof(double));
   ppaerr   = (double *)calloc(nbin*npuls+1, sizeof(double));
   inarru   = (double *)calloc(nbin*npuls+1, sizeof(double));
   inarrq   = (double *)calloc(nbin*npuls+1, sizeof(double));

   // Running Mean of baseline for I,Q,U,V-Stokes
   nrun = ceil((double)(nbin)/(2.0*(double)(win[1]-win[0])));
   rbin = nbin/nrun;

   fprintf (stderr," nrun %d  rbin %d  npuls %d\n", nrun, rbin, npuls);

   runmean = (double *)calloc(rbin*npuls, sizeof(double));


   stkpuls (sftdispi, pulstki, nbin, npuls, nsize-nskip, period, int_time, 0, 
              npuls);
//   runbase (pulstki, runmean, npuls, nbin, rbin, nrun, win[0], win[1]);
//   basesub_pl (pulstki, runmean, npuls, nbin, nrun);


   stkpuls (sftdispq, pulstkq, nbin, npuls, nsize-nskip, period, int_time, 0,
              npuls);
//   runbase (pulstkq, runmean, npuls, nbin, rbin, nrun, win[0], win[1]);
//   basesub_pl (pulstkq, runmean, npuls, nbin, nrun);


   stkpuls (sftdispu, pulstku, nbin, npuls, nsize-nskip, period, int_time, 0,
              npuls);
//   runbase (pulstku, runmean, npuls, nbin, rbin, nrun, win[0], win[1]);
//   basesub_pl (pulstku, runmean, npuls, nbin, nrun);


   stkpuls (sftdispv, pulstkv, nbin, npuls, nsize-nskip, period, int_time, 0,
              npuls);
//   runbase (pulstkv, runmean, npuls, nbin, rbin, nrun, win[0], win[1]);
//   basesub_pl (pulstkv, runmean, npuls, nbin, nrun);


   // Estimating L-Pol, PPA of single pulses and folded profiles

   for (jj=0;jj<nbin;jj++) {
      foldstki[jj] = 0.0;
      foldstkq[jj] = 0.0;
      foldstku[jj] = 0.0;
      foldstkv[jj] = 0.0;
   }

   for (ii=0;ii<npuls;ii++) {

      //Estimating rms in U and Q
      for (jj=bbase;jj<ebase;jj++) {
         inarru[jj-bbase] = pulstku[ii*nbin+jj];
         inarrq[jj-bbase] = pulstkq[ii*nbin+jj];
      }
      meanrms (inarru, ebase-bbase, minstat);
      rmsu = minstat[1];

      meanrms (inarrq, ebase-bbase, minstat);
      rmsq = minstat[1];

      for (jj=0;jj<nbin;jj++) {
         indx = ii*nbin+jj;

         foldstki[jj]+=pulstki[indx]/nbin;
         foldstkq[jj]+=pulstkq[indx]/nbin;
         foldstku[jj]+=pulstku[indx]/nbin;
         foldstkv[jj]+=pulstkv[indx]/nbin;

         pulstkl[indx] = sqrt(pulstkq[indx]*pulstkq[indx]+
                                    pulstku[indx]*pulstku[indx]);

         stkppa[indx] = 0.5*atan2(pulstku[indx],pulstkq[indx])/M_PI*180;

         ppaerr[indx] = 0.5*sqrt(pow(pulstkq[indx]*rmsq,2) +
                        pow(pulstku[indx]*rmsu,2.0))/pow(pulstkl[indx],2.0);
      }
   }

   runbase (pulstkl, runmean, npuls, nbin, rbin, nrun, win[0], win[1]);
   basesub_pl (pulstkl, runmean, npuls, nbin, nrun);

/*
   fpt = fopen("outpolchk.dat","w");

   for (ii=0;ii<npuls;ii++) {
      for (jj=0;jj<nbin;jj++)
         fprintf (fpt,"%d %d %lf %lf %lf %lf\n", ii, jj, pulstki[ii*nbin+jj], pulstkq[ii*nbin+jj], pulstku[ii*nbin+jj], pulstkv[ii*nbin+jj]);
   }
   fclose(fpt);
*/

   //Writing out the stokes file and header
   sprintf (outstokesfile, "%s_iquv.bin", dedispfile);
   sprintf (outstokesascii, "%s_iquv.ascii", dedispfile);

   fpt = fopen(outstokesfile,"wb");
   fpta = fopen(outstokesascii,"w");

   fprintf (fpta, "#Npulse Nbin   I   Q   U   V\n");

   for (ii=0;ii<npuls;ii++) {
      for (jj=0;jj<nbin;jj++) {
         dwrt[0] = pulstki[ii*nbin+jj];
         fwrite (dwrt, sizeof(float), 1, fpt);

         dwrt[0] = pulstkq[ii*nbin+jj];
         fwrite (dwrt, sizeof(float), 1, fpt);

         dwrt[0] = pulstku[ii*nbin+jj];
         fwrite (dwrt, sizeof(float), 1, fpt);

         dwrt[0] = pulstkv[ii*nbin+jj];
         fwrite (dwrt, sizeof(float), 1, fpt);
      }
      for (jj=0;jj<nbin;jj++)
         fprintf (fpta, "%d %d %lf %lf %lf %lf\n", ii, jj, pulstki[ii*nbin+jj], 
               pulstkq[ii*nbin+jj], pulstku[ii*nbin+jj], pulstkv[ii*nbin+jj]);
   }
   fclose(fpt);
   fclose(fpta); 


   sprintf (outfold, "%s_iquv_fold.dat", dedispfile);

   fpt = fopen(outfold,"w");

   fprintf (fpt, "#Nbin   I   Q   U   V\n");

   for (jj=0;jj<nbin;jj++)
      fprintf (fpt, "%d %lf %lf %lf %lf \n", jj, foldstki[jj], foldstkq[jj],
                                              foldstku[jj], foldstkv[jj]);

   fclose(fpt);


   sprintf (outpolfile, "%s_pol.bin", dedispfile);

   fpt = fopen(outpolfile,"wb");

   for (ii=0;ii<npuls;ii++) {
      for (jj=0;jj<nbin;jj++) {
         dwrt[0] = pulstkl[ii*nbin+jj];
         fwrite (dwrt, sizeof(float), 1, fpt);

         dwrt[0] = pulstkv[ii*nbin+jj];
         fwrite (dwrt, sizeof(float), 1, fpt);
         
         dwrt[0] = stkppa[ii*nbin+jj];
         fwrite (dwrt, sizeof(float), 1, fpt);

         dwrt[0] = ppaerr[ii*nbin+jj];
         fwrite (dwrt, sizeof(float), 1, fpt);
      }
   } 
   fclose(fpt);

   sprintf (outhdrfile, "%s_iquvhdr.dat", dedispfile);

   fpt = fopen(outhdrfile, "w");

   fprintf (fpt, "npuls %d\n", npuls);
   fprintf (fpt, "nbin %d\n", nbin);
   fprintf (fpt, "bwin %d\n", win[0]);
   fprintf (fpt, "ewin %d\n", win[1]);
   fprintf (fpt, "npuls_avg %d\n", npuls_avg);
   fprintf (fpt, "nskip %d\n", nskip);
   fclose(fpt);


   return 0;
}