/* vg1pws_sa_list.c
 *
 * Example 'C' code to extract Voyager 1 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 vg1pws_sa_list VG1PWS_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 verify1.txt.
 * Typical usage would be:
 *   vg1pws_sa_list < T790305.DAT
 */

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


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

int vg1pwssd (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.
 *
 * The calibration file VG1PWSCL.TAB is read from the current directory.
 *
 * 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 */

  /* 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[0][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 */
  }

  /* 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;

} /* vg1pwssd - 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 */
  const char *calfile = "VG1PWSCL.TAB"; /* Calibration assumed in local dir */
  FILE* outfp;

  if (!(outfp = fopen("verify1.txt", "wb"))){
    fprintf (stderr, "ERROR opening verify1.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 (vg1pwssd (calfile, rec, sd)) {
      fprintf (stderr, "ABORTING vg1pws_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 - vg1pws_sa_list example program */