/* ----------------------------------------------------------------
   PWS16SVR written by R. Baines 2/90.
    based on FORTRAN code by L. Granroth.

   This module contains the Begin, Command, Data and End functions
   for internal data processing of raw data read from file.

      Begin   - opens communication.

      Command - instructs server on what data is req'd and how to
                process it.

      Data    - returns one channel of data.

      End     - closes communication.


   Since there is alot of code, the Command function has been
   broken down into several functions to make it more readable
   and to reduce repititous code segments.
   ---------------------------------------------------------------- */


/* Includes */
#include <stdio.h>
#include <stdlib.h>
#include <mem.h>
#include <math.h>
#include <alloc.h>
#include <string.h>
#include <io.h>
#include <fcntl.h>
#include "pwstruct.h"
#include "pwsglobl.h"
#include "proto.h"


/* Macros for this module */
#define sfail 411077030.0
#define big 10.0e30
#define iseed 1313137L
#define dat_length 454


/* Internal Global Vars - some used in PWS16SRB */
int ibin,
    mode,
    ichan = 0,
    numbin,
    maxbin,
    jstat,
    idatmin[NUM_CHAN],
    idatmax[NUM_CHAN],
    jsc,
    itype=1,
    kyr,
    kdoy,
    ndays,
    kyear,
    kleap,
    kountwait,
    index,
    lyr,
    ldoy,
    lhr;

boolean pw_found = FALSE,
        fg_found = FALSE,
        data_found = FALSE,
        cmd_pw,
        found;

float span,
      binwidth,
      offset,
      cal[NUM_CHAN+1][256],
      tonl[8] = { 2.0, 1.0, -1.0, -2.0, -3.0, 1.0, 2.0, 1.0 },
      bw[NUM_CHAN][2] =
         { { 2.99, 2.16 }, { 3.77, 3.58 }, { 7.50, 4.50 }, { 10.6, 10.7 },
	   { 13.3, 13.8 }, { 29.8, 28.8 }, { 59.5, 39.8 }, { 106., 75.9 },
	   { 133., 75.9 }, { 211., 151. }, { 298., 324. }, { 421., 513. },
	   { 943., 832. }, { 2110, 1260 }, { 4210, 2400 }, { 5950, 3800 } };


FILE *fileunit,
     *datunit;

char filename[LLEN];

double tbeg,
       tend,
       tyme=0;

cd_type cd[64];



/* Internal Function Prototypes */
void process_data(int, int*, float*, float*, float*);
void pw_set(int, int);
void pw_range(int, int, int, long, float);
void pw_get(float *);
int getcd(void);
void cd_read(void *, size_t, size_t, FILE *);
void cd_read_rest(size_t, size_t);
int doy_to_md (int);
int leap_doy_to_md (int);

/* External Function Prototype for reading compressed raw data files */
int LZOPEN(char *);
int LZREAD(void *, int *);
int LZCLOSE(void);






/*---------------------------------------------------------------------

  pws16_server_begin by R. Baines 2/90.
  to initialize server connections for the PWS16PC program.

  ---------------------------------------------------------------------*/

void pws16_server_begin(void)
{
  /* Set Initial Values for cmd, Zmin and Zmax */
  printf ("\nPWS16C version %s\n",VERSION);
  cmd_reset();
  puts ("\nEnter plot parameters (or help)");
}


/*---------------------------------------------------------------------

  pws16_server_command written by R. Baines
  Nonzero is returned after a successfully transmitted request and
  acknowledgment, otherwise zero is returned.

  NOTE: ALL processing is done in the LOOP here.

  ---------------------------------------------------------------------*/

int pws16_server_command (void)
{

   float data[NUM_CHAN], sum[NUM_CHAN], peak[NUM_CHAN];
   int num[NUM_CHAN], ichan, numsend;

   if (Cmd.survey)
       error(2);

   if (Cmd.length < 0)
      error(3);

   printf("000 Out Of %3d", Cmd.numbin);
   jstat = 2;
   cmd_pw = !Cmd.fgonly;
   pw_found = pw_found && Cmd.replot;
   fg_found = fg_found && Cmd.replot;
   data_found = pw_found || fg_found;

   span = (float) Cmd.span;
   if ((Cmd.numbin < 1) || (Cmd.numbin > 900))
      numbin = 900;
   else
      numbin = Cmd.numbin;

   maxbin = numbin - 1;
   binwidth = span / (float) numbin;

   if (cmd_pw && !pw_found) {
      pw_set(Cmd.units, Cmd.spacecraft);
      pw_range(Cmd.year, Cmd.day, Cmd.hour, Cmd.minute, span);
      pw_get(data);
      pw_found = jstat == 1;
      if (pw_found) {

         data_found = TRUE;
         ibin = 0;
         for (ichan=0; ichan<NUM_CHAN; ichan++) {
            sum[ichan] = 0.0;
            peak[ichan] = 0.0;
            num[ichan] = 0.0;
            idatmin[ichan] = 32000;
            idatmax[ichan] = -32000;
         }

         while (jstat == 1) {
            process_data(offset/binwidth, num, peak, sum, data);
            pw_get(data);
         }
         process_data(numbin, num, peak, sum, data);

         /* if units 5 or 3 then rms averaging */
         if ((Cmd.units == 3) || (Cmd.units == 5)) 

            for (ichan=0; ichan<NUM_CHAN; ichan++) {

               for (ibin=0; ibin<=maxbin; ibin++) {

                  /* divide 2 for logs */
                  ipks[ichan][ibin] = ipks[ichan][ibin]/2; 
                  iavg[ichan][ibin] = iavg[ichan][ibin]/2;
               }

               idatmin[ichan] = idatmin[ichan]/2;
               idatmax[ichan] = idatmax[ichan]/2;
            }

      } /* end if pw_found */
   } /* end if cmd_pw */

   numsend = numbin + numbin % 2;
   Dat.length = 4 + numsend/2;
   Dat.spare0 = 0;
   Dat.spare1 = 0;

   Iteration = 0;

   return (TRUE);


} /* pws16_server_command */


/*---------------------------------------------------------------------

  pws16_server_data written by R. Baines 2/90
  does nothing really except compute whether this data is
  peaks or averages and the correct position of the channel.

  ---------------------------------------------------------------------*/

int pws16_server_data (void)
{
   int i;

   if (Iteration/2 >= Cmd.nchan) {
      exit_graph();
      return(FALSE);
   }

   Dat.position = Cmd.chan[Iteration/2]-1;
   Dat.datmin = idatmin[Dat.position];
   Dat.datmax = idatmax[Dat.position];
   if (Cmd.peaks && ISODD(Iteration)) {
	  Dat.flag = PEAKS;
      for (i=0; i<900; i++)
         Dat.data[i] = ipks[Dat.position][i];
   }
   else {
      for (i=0; i<900; i++)
         Dat.data[i] = iavg[Dat.position][i];
      Dat.flag = AVERAGES;
   }

   Dat.position++;
   Iteration += (Cmd.peaks) ? 1 : 2;
   return(TRUE);

} /* pws16_server_data */


/*---------------------------------------------------------------------

  pws16_server_end written by R. Baines 2/90
  does nothing.

  ---------------------------------------------------------------------*/

void pws16_server_end (void)
{
  return;
} /* pws16_server_end */




/* ------------------------------------------------------------------

   process_data written by R. Baines 3/90
   this is a segment of code extracted from pws16_server_command to
   make the code more readable and reduce repititous code.

   ------------------------------------------------------------------ */

void process_data(int jbin, int num[], float peak[], float sum[], float data[])
{
   int ichan, j;
   float l;
   static float tcyc[32] =
      { 4, 4, 9.6, 38.4, 76.8, 96, 4, 4, 4, 4, 4, 4, 4, 4, 4,
        4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4 };


   if (jbin != ibin) {

      for (ichan=0; ichan<NUM_CHAN; ichan++) {

         if (num[ichan] > 0) {

            if (Cmd.units == 1) {
	       iavg[ichan][ibin] = NINT(sum[ichan] / (float) num[ichan]);
               ipks[ichan][ibin] = peak[ichan];
            }
            else {
               iavg[ichan][ibin] = NINT(log10(sum[ichan] / (float) num[ichan])* 1000.0);
               ipks[ichan][ibin] = NINT(log10(peak[ichan]) * 1000.0);
            }

            if (ipks[ichan][ibin] > idatmax[ichan])
               idatmax[ichan] = ipks[ichan][ibin];

            if (iavg[ichan][ibin] < idatmin[ichan])
               idatmin[ichan] = iavg[ichan][ibin];

         }
         else /* num[ichan] <= 0 */ {
            iavg[ichan][ibin] = 0;
            ipks[ichan][ibin] = 0;
         }

         if (data[ichan] > 0.0) {
            peak[ichan] = data[ichan];
            sum[ichan] = data[ichan];
            num[ichan] = 1;
         }
         else {
            peak[ichan] = 0.0;
            sum[ichan] = 0.0;
            num[ichan] = 0;
         }
      } /* for */

      /* Fill In iavg and ipks with values - fills in gap */
      if (jbin > ibin+1) {
		 l = ibin + NINT(tcyc[mode]/binwidth);
         if (l > jbin-1)
            l = jbin-1;
         for (j=ibin+1; j<=l; j++)
            for (ichan=0; ichan<NUM_CHAN; ichan++) {
               iavg[ichan][j] = iavg[ichan][ibin];
               ipks[ichan][j] = ipks[ichan][ibin];
            }
         for (j=l+1; j < jbin; j++)
            for (ichan=0; ichan<NUM_CHAN; ichan++) {
               iavg[ichan][j] = 0;
               ipks[ichan][j] = 0;
            }
      }
      ibin = jbin;
      printf("\r%03d", jbin);
   }

   else  /* if (jbin == ibin) */
       for (ichan=0; ichan<NUM_CHAN; ichan++)
          if (data[ichan] > 0.0) {
            if (data[ichan] > peak[ichan])
               peak[ichan] = data[ichan];
            sum[ichan] += data[ichan];
            num[ichan]++;
         }
}




/* ----------------------------------------------------------

   pw_set written by R. Baines 2/90

   reads in calibration tables and initializes them based
   upon the units used in the data request.
   ---------------------------------------------------------- */

void pw_set(int iset, int isc)
{
   int j, ichan;
   int dummy;
   FILE *fileunit;

   /* returns if data to be seen has already been processed */
   if ((iset == itype) && (isc == jsc))
      return;

   strcpy(filename, Cal_Path);
   if (isc == 2)
      strcat(filename,"vg2pwscl.tab");
   else
      strcat(filename,"vg1pwscl.tab");

   itype = iset;
   jsc = isc;

   /* create calibration table */
   setmem (cal[0],256, -99);

   if (!(fileunit = fopen(filename, "r")))
      error(9);


   for (ichan = 0; ichan < 256; ichan++)
      fscanf (fileunit, "%d,%E,%E,%E,%E,%E,%E,%E,%E,%E,%E,%E,%E,%E,%E,%E,%E", &dummy,
	 &cal[1][ichan],  &cal[2][ichan],  &cal[3][ichan],  &cal[4][ichan],
	 &cal[5][ichan],  &cal[6][ichan],  &cal[7][ichan],  &cal[8][ichan],
	 &cal[9][ichan],  &cal[10][ichan], &cal[11][ichan], &cal[12][ichan],
	 &cal[13][ichan], &cal[14][ichan], &cal[15][ichan], &cal[16][ichan]);

   fclose (fileunit);
/*
   for (j=1; j<NUM_CHAN+1; j++)
      for (ichan=0; ichan<256; ichan++)
	 read(fileunit, &cal[j][ichan], sizeof(float));
   close(fileunit);
*/

   /* verify if valid type and process table if neccesary
      else set type to raw data                                 */
   if ((itype > 1) && (itype < 8)) {
      if (itype == 3)
	 for (j=1; j<NUM_CHAN+1; j++)
	    for (ichan=0; ichan<256; ichan++)
	       cal[j][ichan] = cal[j][ichan] * cal[j][ichan];

      else if (itype == 4)
	 for (j=1; j<NUM_CHAN+1; j++)
	    for (ichan=0; ichan<256; ichan++)
	       cal[j][ichan] = cal[j][ichan] / 7.07;

      else if ((itype == 5) || (itype == 6))
	 for (j=1; j<NUM_CHAN+1; j++)
	    for (ichan=0; ichan<256; ichan++)
	       cal[j][ichan] = cal[j][ichan] * cal[j][ichan] / (7.07 * 7.07);

      else if (itype == 7)
	 for (j=1; j<NUM_CHAN+1; j++)
	    for (ichan=0; ichan<256; ichan++)
	       cal[j][ichan] = cal[j][ichan] * cal[j][ichan] / (7.07*7.07*377);
   }
   else
	  itype = 1;

}

/* ----------------------------------------------------------

   pw_range written by R. Baines 2/90

   sets initial values to tbeg and tend.  These vars are used
   to determine if the raw data is from the correct time to
   start processing (compare to TBEG) and to stop processing
   (compare to TEND).  Gets first segment of data in most
   probably in CD format.

   ---------------------------------------------------------- */

void pw_range(int iyr, int idoy, int ihr, long imn, float span)
{

   float secl;

   tbeg = time2sec(iyr, idoy, ihr, (int) imn, 0.0);
   tend = tbeg + (float) span;
   sec2time(tend, &lyr, &ldoy, &lhr, &(int)imn, &secl);

   kyr = iyr;
   kdoy = idoy;

   getcd();
}

/* ----------------------------------------------------------

   pw_get written by R. Baines 2/90

   get data from cd read in.  Fix upper eight channels.
   compute time of cd chunk. Return array of values to
   pws_server_command.  This function is called MANY times.
   ---------------------------------------------------------- */

void pw_get(float data[])
{

   int ichan;

   /* see if time to quit */
   if (tyme > tend) {
      cd_read_rest(sizeof(cd_type), 64);
      if (Compress) {
	 LZCLOSE();
         Compress = FALSE;
      }
      else
	 fclose(datunit);

      jstat = 2;
      offset = -1;
      return;
   }

   offset = tyme - tbeg;
   mode = (int) cd[index].mode;

   /* read in data from cd */
   if (cd[index].ef) {
      if (itype == 1)
         for (ichan=0; ichan<NUM_CHAN; ichan++)
            data[ichan] = cd[index].data[ichan];

      else if (itype < 6)
         for (ichan=0; ichan<NUM_CHAN; ichan++)
            data[ichan] = cal[ichan][cd[index].data[ichan]];

      else
         for (ichan=0; ichan<NUM_CHAN; ichan++)
            data[ichan] = cal[ichan][cd[index].data[ichan]]/ bw[ichan][(jsc-1)];

      jstat = 1;
      index++;

      if (index > 64) {
         index = 0;
	 cd_read(cd, sizeof(cd_type), 64, datunit);
      }

      if ((cd[index].org.ef.time == 0.0) && (cd[index].org.ef.mfds == 0))
         if (!getcd()) {
            jstat=2;
            return;
         }

      tyme = cd[index].org.ef.time;
   }
   else {
   /* Fix the upper 8 channels on voyager 2 */
   /* note: this approx. correction is good only for part of '78 and later */

      if ((jsc == 2) && (tyme >= sfail))
         for (ichan=8; ichan<NUM_CHAN; ichan++) {
            if (cd[index].data[ichan] > 0) {

	       if (cd[index].data[ichan] < 64)
                  cd[index].data[ichan] = 64;

               if (cd[index].data[ichan] <= 72)
		  cd[index].data[ichan] = tonl[ichan-8] - 530.4 + 8.6 * (float) cd[index].data[ichan];
               else
		  cd[index].data[ichan] = tonl[ichan-8] + 20.113 + 0.99 * (float) cd[index].data[ichan];
            }
         }

      if (itype == 1)
         for (ichan=0; ichan<NUM_CHAN; ichan++)
	    data[ichan] = max(0, cd[index].data[ichan]);

      else if (itype < 6)
         for (ichan=0; ichan<NUM_CHAN; ichan++)
            data[ichan] = cal[ichan+1][max(0, cd[index].data[ichan])];

      else
         for (ichan=0; ichan<NUM_CHAN; ichan++)
            data[ichan] = cal[ichan+1][max(0, cd[index].data[ichan])] / bw[ichan][(jsc-1)];

      jstat = 1;
      index++;

      if (index >= 64) {
         index = 0;
	 cd_read(&cd[0], sizeof(cd_type), 64, datunit);
      }

      if ((cd[index].org.ef.time == 0.0) && (cd[index].org.ef.mfds == 0))
         if (!getcd()) {
            jstat=2;
            return;
         }


      tyme = (double) (cd[index].org.cd.millisec / 1000.0)
             + (double) cd[index].org.cd.second
             + (cd[index].org.cd.hour -24) * (double) HOUR_SEC
             + kleap * (double) DAY_SEC
             + kyear * (double) YEAR_SEC;
   }

}


/* ------------------------------------------------------------

   getcd written by R. Baines 2/90

   gets a segment of raw data in cd format from a data file(s).
   data files may be in compressed format.

   ------------------------------------------------------------ */

int getcd(void)
{
   char filename[LLEN], name[13];

   /* Search for first file in range */
   if (datunit)
      if (Compress) {
         LZCLOSE();
         Compress = FALSE;
      }
      else
         fclose(datunit);

   datunit = NULL;

   while (!datunit)
   {
      ndays = (kyr % 4 == 0) ? LEAP_YEAR_DAYS : YEAR_DAYS;
      sprintf(name, "pws%0d%.3d.", kyr, kdoy);
      strcpy(filename, Data_Path[Cmd.spacecraft-1]);
      strcat(filename, name);
      strcat(filename, Extension);

      /* Attempt to open raw data file */
      datunit = fopen(filename, "rb");

      /* if unsuccesful and we have a extension for compressed files
         then attempt to open COMPRESSED raw data file.               */

      if (!datunit)
      {
	 strcpy(filename, Data_Path[Cmd.spacecraft-1]);
	 strcat(filename, name);
	 strcat(filename, ZExtension);
	 if (LZOPEN(filename))
	    Compress = TRUE;
      }

      if (!datunit)
      {
	 ndays = (kyr % 4 == 0) ? LEAP_YEAR_DAYS : YEAR_DAYS;
	 sprintf(name, "pw%0d%0d%.3d.", Cmd.spacecraft, kyr, kdoy);
	 strcpy(filename, Data_Path[Cmd.spacecraft-1]);
	 strcat(filename, name);
	 strcat(filename, Extension);
         datunit = fopen(filename, "rb");
      }

      if (!datunit) {
	 strcpy(filename, Data_Path[Cmd.spacecraft-1]);
	 strcat(filename, name);
	 strcat(filename, ZExtension);
	 if (LZOPEN(filename))
	    Compress = TRUE;

      }

      if (!datunit)
      {
	 ndays = (kyr % 4 == 0) ? LEAP_YEAR_DAYS : YEAR_DAYS;
         if (ndays == YEAR_DAYS)
            sprintf(name, "T%0d%.4d.", kyr, doy_to_md (kdoy));
         else
            sprintf(name, "T%0d%.4d.", kyr, leap_doy_to_md (kdoy));
	 strcpy(filename, Data_Path[Cmd.spacecraft-1]);
	 strcat(filename, name);
	 strcat(filename, Extension);
         datunit = fopen(filename, "rb");
      }

      if (!datunit) {
	 strcpy(filename, Data_Path[Cmd.spacecraft-1]);
	 strcat(filename, name);
	 strcat(filename, ZExtension);
	 if (LZOPEN(filename))
	    Compress = TRUE;

      }

      kdoy++;
      if (datunit)
	 break;

      if (kdoy > ndays) {
	 kdoy = 1;
	 kyr++;
      }

      if ((kyr >= lyr) && (kdoy > ldoy))
	 return(FALSE);
   }
   tyme = 0.0;
   kyear = kyr-65;
   kleap = kyear/4;
   kountwait = 0;

   /* find first physical record */
   if (datunit) {
      while (tyme < tbeg) {
	 cd_read(cd, sizeof(cd_type), 64, datunit);
         index = 63;

         while ((cd[index].org.ef.time == 0.0) && (cd[index].org.ef.mfds == 0))
            index -=1;

         if (index < 0)
            error(4);

         if (cd[index].ef)
	    tyme = cd[index].org.ef.time;
         else
	    tyme = (double) ((float) cd[index].org.cd.millisec / 1000
                 + (float) cd[index].org.cd.second)
                 + (double) ((float)(cd[index].org.cd.hour -24) * HOUR_SEC
                 + (double) (kleap * DAY_SEC))
                 + (double) ((double) kyear * YEAR_DAYS * DAY_SEC);
      }

      /* find logical record */
      tyme = 0.0;
      index = 0;

      while (tyme < tbeg) {
         if (cd[index].ef)
	    tyme = cd[index].org.ef.time;
         else
	    tyme = (double) ((float) cd[index].org.cd.millisec / 1000
                 + (float) cd[index].org.cd.second)
                 + (double) ((float)(cd[index].org.cd.hour -24) * HOUR_SEC
                 + (float) (kleap * DAY_SEC))
                 + (double) (kyear * YEAR_DAYS * DAY_SEC);
         index += 1;
      }
   }
   else
      tyme = tend;

   return(TRUE);
}


/* ----------------------------------------------------------

   cd_read written by R. Baines 2/90

   actually do physical read of raw data or pass reading to
   the compress module to decompress data from a compressed
   data file.
   ---------------------------------------------------------- */

void cd_read(void *ptr, size_t size, size_t n, FILE *stream)
{
   const int len = size * n;
   cd_type buf[64];

   if (Compress)
      LZREAD(&buf, &len);
   else
      fread(buf, size, n, stream);

   memcpy(ptr, buf, len);
}

/* ----------------------------------------------------------

   cd_read_rest written by R. Baines 2/90

   read rest of compressed data file to handle any problems
   that may occur if we leave it unfinished.
   ---------------------------------------------------------- */

void cd_read_rest(size_t size, size_t n)
{
   const int len = size * n;
   cd_type buf[64];

   while(LZREAD(&buf, &len) == len);
}

int doy_to_md (int l_doy)
{
   int md;

   if ((l_doy >= 1) && (l_doy <=31))
      md = 100 + l_doy;
   else if ((l_doy >= 32) && (l_doy <= 59))
      md = 200 + (l_doy - 31);
   else if ((l_doy >= 60) && (l_doy <= 90))
      md = 300 + (l_doy - 59);
   else if ((l_doy >= 91) && (l_doy <= 120))
      md = 400 + (l_doy - 90);
   else if ((l_doy >= 121) && (l_doy <= 151))
      md = 500 + (l_doy - 120);
   else if ((l_doy >= 152) && (l_doy <= 181))
      md = 600 + (l_doy - 151);
   else if ((l_doy >= 182) && (l_doy <= 212))
      md = 700 + (l_doy - 181);
   else if ((l_doy >= 213) && (l_doy <= 243))
      md = 800 + (l_doy - 212);
   else if ((l_doy >= 244) && (l_doy <= 273))
      md = 900 + (l_doy - 243);
   else if ((l_doy >= 274) && (l_doy <= 304))
      md = 1000 + (l_doy - 273);
   else if ((l_doy >= 305) && (l_doy <= 334))
      md = 1100 + (l_doy - 304);
   else if ((l_doy >= 335) && (l_doy <= 365))
      md = 1200 + (l_doy - 334);
   else
      md = 101;

   return (md);
}

int leap_doy_to_md (int l_doy)
{
   int md;

   if ((l_doy >= 1) && (l_doy <=31))
      md = 100 + l_doy;
   else if ((l_doy >= 32) && (l_doy <= 60))
      md = 200 + (l_doy - 31);
   else if ((l_doy >= 61) && (l_doy <= 91))
      md = 300 + (l_doy - 60);
   else if ((l_doy >= 92) && (l_doy <= 121))
      md = 400 + (l_doy - 91);
   else if ((l_doy >= 122) && (l_doy <= 152))
      md = 500 + (l_doy - 121);
   else if ((l_doy >= 153) && (l_doy <= 182))
      md = 600 + (l_doy - 152);
   else if ((l_doy >= 183) && (l_doy <= 213))
      md = 700 + (l_doy - 182);
   else if ((l_doy >= 214) && (l_doy <= 244))
      md = 800 + (l_doy - 213);
   else if ((l_doy >= 245) && (l_doy <= 274))
      md = 900 + (l_doy - 244);
   else if ((l_doy >= 275) && (l_doy <= 305))
      md = 1000 + (l_doy - 274);
   else if ((l_doy >= 306) && (l_doy <= 335))
      md = 1100 + (l_doy - 305);
   else if ((l_doy >= 336) && (l_doy <= 366))
      md = 1200 + (l_doy - 335);
   else
      md = 101;

   return (md);
}