/* vg2pws_sa_list.c
 *
 * Example 'C' code to extract Voyager 2 PWS full-resolution spectrum
 * analyzer samples and list them as calibrated spectral density values
 * along with spacecraft event time and telemetry mode information.
 * This code constitutes a complete application and should work on any
 * system with an ANSI C compiler regardless of machine byte order.
 *
 * On UNIX systems this example can be compiled with a command like:
 *   cc -O -o vg2pws_sa_list VG2PWS_SA_LIST.C
 *
 * The calibration file VG1PWSCL.TAB is assumed to be in the current
 * directory.
 *
 * Input is read from stdin and output is to a file named verify2.txt.
 * Typical usage would be:
 *   vg2pws_sa_list < T790305.DAT
 *
 * Adjustments for decreased sensitivity in the upper 8 SA channels are
 * applied for data on or after September 24, 1977.
 */

#include <stdlib.h>
#include <stdio.h>
#include "VGPWS_SA_LIST.H"

/* Helper function to determine data correction coefficent set 
 *
 * Switch failure occured on 1977-09-24 JPL didn't stop trying different 
 * fixes until about Dec. 1st 1977.  A spontaneous state change in
 * the channels was observed in 2006 and seemed to revert back to
 * a state that was observered in 1977 
 */
static const double tonl[] = {2.0, 1.0, -1.0, -2.0, -3.0, 1.0, 2.0, 1.0};
static const double cross_pt[]    = {76 - 4, 76 + 10};
static const double cor_fac_low[] = {-530.4, -650.8 };
static const double cor_fac_hi[]  = {20.113, 6.253  };

int datcoridx(const unsigned char* buf){
  
  int nYear = PWSA_y1900(buf) + 1900;
  int nHoy  = PWSA_hoy(buf);
  int nSoh  = PWSA_soh(buf); 
  double rDoy = ( nSoh + (nHoy * 3600. )) / 86400.0;

  /* Use correction 0 from: 
   *   1977-267T00:47 to 1977-283T16:00
   *   1978-010T20:04 to 2006-324T20:50 
	*/
  if( (nYear == 1977) && 
      (rDoy >= 267.032638889) && (rDoy < 283.666666667)) return 0;
  if( (nYear == 1978) && (rDoy >= 10.835995)) return 0;
  if( (nYear > 1978) && (nYear < 2006) ) return 0;
  if( (nYear == 2006) && (rDoy < 324.868055555)) return 0;
  
  /* Use correction 1 from:
   *   1977-283T16:00 to 1977-312T20:12
   *   1977-335T21:54 to 1978-010T20:04    
   *   2006-324T20:50 to present (at least 2016-245) 
	*/
  if(nYear == 1977){
     if((rDoy >= 283.666666667)&&(rDoy < 312.841666667)) return 1;
	  if(rDoy >= 335.912222222) return 1;
  }
  if((nYear == 1978) && (rDoy < 10.835995)) return 1;
  
  if((nYear == 2006) && (rDoy >= 324.868055555)) return 1;
  if(nYear > 2006) return 1;
  
  /* Use no correction otherwise */
  return -1;
}



/* The following calibration function could be saved out in a separate
 * source file for more convenient use in other applications.
 */

int vg2pwssd (const char* calfile, const unsigned char *buf, double *sd) {

/* Given a pointer to a 48-byte buffer containing a Voyager 1 PWS
 * full-resolution uncalibrated data record and a pointer to a
 * 16-element double precision array, fill the array with the
 * corresponding calibrated spectral density (V**2/m**2/Hz) values.
 *
 * Return 0 on success and non-zero on failure.
 */

  static int loadtable = 1;     /* flag to load calibration table */
  FILE *calfp;                  /* calibration file pointer */
  static double cal[256][16];   /* calibration table */
  double value;                 /* calibrated value */
  int dn[16];                   /* array of uncalibrated samples */
  int i, j, idata;              /* loop indices and temp variables */
  double data, y0, y1;          /* used in linear interpolation */

  double ddn;

  /* Load calibration table once */
  if (loadtable) {

    if (!(calfp = fopen(calfile, "r"))) {
      fprintf (stderr, "ERROR opening %s\n", calfile);
      return -1;
    }
    for (j = 0; j < 256; j++) {
      if (!fscanf (calfp, " %3d", &idata) || (idata != j)) {
        fprintf (stderr, "ERROR loading %s\n", calfile);
        return -1;
      }
      for (i = 0; i < 16; i++) {
        if (!fscanf (calfp, ",%8lE", &value)) {
          fprintf (stderr, "ERROR reading %s\n", calfile);
          return -1;
        }
        /* Convert volts to spectral density (V**2/m**2/Hz) and store */
        cal[j][i] = (value * value) / (PWS_elength * PWS_elength) /
          PWSA_bw[1][i];
      } /* for all 16 channels */
    } /* for all 256 uncalibrated values */

    fclose (calfp);
    loadtable = 0;

  } /* if table needs to be loaded */

  /* Extract raw uncalibrated samples (data numbers, dn) */
  for (i = 0; i < 16; i++) {
    dn[i] = PWSA_dn(buf,i);
    if (dn[i] & 0x8000) dn[i] = 0;      /* zero flagged interference */
  }
  
  /* Determine set of correcton coeffients to use */
  int iCor = datcoridx(buf);
	 
  if(iCor >= 0){
		 
    switch (PWSA_mode(buf)){

      /* 8-bit modes */
      case CR2: case CR3: case CR4: case CR5: case CR6: case CR1:
      case GS10A: case GS3: case GS7: case GS6: case OC2: case OC1:
      case GS10: case GS8:
        for (i = 8; i < 16; i++) {
          if (dn[i] > 0) {

            if (dn[i] < 64) dn[i] = 64;
            if (dn[i] <= cross_pt[iCor])
              dn[i] = (int)(tonl[i - 8] + cor_fac_low[iCor] + 8.6 * dn[i]);
            else
              dn[i] = (int)(tonl[i - 8] + cor_fac_hi[iCor] + 0.99 * dn[i]);
				  
            if(dn[i] > 255) dn[i] = 255;
          }
        }
        break;

      /* 10-bit sums of 4 instrument samples */
      case CR5A: case UV5A:
         for (i = 8; i < 16; i++) {
           if (dn[i] > 0) {
             ddn = dn[i] / 4.0;

             if (ddn < 64.0) ddn = 64.0;
             if (ddn <= cross_pt[iCor])
                ddn = tonl[i - 8] + cor_fac_low[iCor] + 8.6 * ddn;
             else
                ddn = tonl[i - 8] + cor_fac_hi[iCor] + 0.99 * ddn;

             if(ddn > 255) ddn = 255;
             dn[i] = (int)(ddn * 4.0);
           }
         }
         break;
      default:
        fprintf (stderr, "ERROR unexpected telemetry mode %d\n",
                 PWSA_mode(buf));
        return -1;

    } /* switch on telemetry mode */

  } /* if adjust upper 8 ADC counts */

  /* Extract telemetry mode information and handle 10-bit sum cases */
  switch (PWSA_mode(buf)) {

    /* 8-bit modes */
    case CR2: case CR3: case CR4: case CR5: case CR6: case CR1:
    case GS10A: case GS3: case GS7: case GS6: case OC2: case OC1:
    case GS10: case GS8:
      for (i = 0; i < 16; i++) sd[i] = cal[dn[i]][i];
      break;
    /* 10-bit sums of 4 instrument samples */
    case CR5A: case UV5A:
      for (i = 0; i < 16; i++) {
        data = (double)dn[i] / 4.0;
        idata = (int)data;
        if (idata == 255) sd[i] = cal[255][i];
        else {
          y0 = cal[idata][i];
          y1 = cal[idata+1][i];
          sd[i] = y0 + (y1 - y0) * (data - (double)idata);
        }
      }
      break;
    default:
      fprintf
        (stderr, "ERROR unexpected telemetry mode %d\n", PWSA_mode(buf));
      return -1;

  } /* switch on telemetry mode */

  return 0;

} /* vg2pwssd - extract calibrated spectral density */


/* The main program . . .
 */

int main (int argc, char *argv[]) {

  unsigned char rec[PWSA_RECSIZE];      /* fixed-length records */
  unsigned long int count = 0;          /* number of records processed */
  int i;                                /* loop index */
  double sd[16];                        /* 16 channels of spectral density */
  char *calfile = "VG2PWSCL.TAB";       /* Calibration assumed in local dir */
  FILE* outfp;

  if (!(outfp = fopen("verify2.txt", "wb"))){
    fprintf (stderr, "ERROR opening verify2.txt for writing\n");
     return 3;
  }

  /* Loop through all available records from standard input */
  while (fread (rec, PWSA_RECSIZE, 1, stdin)) {
    count++;

    fprintf (outfp, "\r\nRECORD %5ld ", count);

    /* Use convenience macros to extract SCET and SCLK */
    fprintf (outfp, " SCET %04d %03d %02d:%02d:%02d.%03d ", PWSA_year(rec),
      PWSA_day(rec), PWSA_hour(rec), PWSA_minute(rec), PWSA_sec(rec),
      PWSA_msec(rec));
    fprintf (outfp, " SCLK %05d.%02d.%03d ", PWSA_mod16(rec), PWSA_mod60(rec),
      PWSA_line(rec));

    /* Telemetry mode and spacecraft number */
    fprintf (outfp, " MODE %4s ", TLM_mode_name[PWSA_mode(rec)]);
    fprintf (outfp, " VG%d\r\n", PWSA_sc(rec));

    /* Calibrated spectral density values */
    if (vg2pwssd (calfile, rec, sd)) {
      fprintf (stderr, "ABORTING vg2pws_sa_list\n");
      return 7;
    }

    for (i = 0; i < 16; i++) {
      fprintf (outfp, " %8.2E", sd[i]);
    }
    fprintf (outfp, "\r\n");

  } /* while data to read */

  fclose(outfp);
  return 0;

} /* main - vg2pws_sa_list example program */