/* STEREO IMPACT SEP SEPT 
** 
** L1 Data Processing
**
** (C) Copyright 2005-2006 Stephan I. Böttcher, CAU Kiel - IEAP
**
** Stephan I. Böttcher     <boettcher@physik.uni-kiel.de>
** Reinhold Müller-Mellin  <mueller-mellin@physik.uni-kiel.de>
**
** $Id: sept_l1.c,v 1.8 2010/05/31 11:49:35 bottcher Exp $
** $Log: sept_l1.c,v $
** Revision 1.8  2010/05/31 11:49:35  bottcher
** add rounding to dec14, test case
**
** Revision 1.7  2010/05/10 08:45:51  bottcher
** Fix bug in float 14 decompression for normalized numbers, version 0x040000
**
** Revision 1.6  2006/12/22 07:22:54  bottcher
** software version 0x030000: double time
**
** Revision 1.5  2006/12/21 08:54:33  bottcher
** Fix bit-lossage in specta unpacking
**
** Revision 1.4  2006/10/13 06:36:50  bottcher
** turn line ends to Unix
**
** Revision 1.3  2006/03/22 15:46:39  bottcher
** init_called for automatic initialization
** verbosity added
**
** Revision 1.2  2006/03/22 12:58:55  bottcher
** Usage documentation added
**
** Revision 1.1  2005/12/27 22:35:38  bottcher
** Initial commit, feature complete.
**
*/

static const unsigned int version = 0x040000;

/*
** version & 0xff000000 : errors
** version & 0x00ff0000 : software version
** version & 0x0000ff00 : calibration version
** version & 0x000000ff : calibration subversion
*/

#include "sept_l1.h"

static struct sept_calibration
{
  int version;
  int temp_offset;
  double bin_width[4][32];
  double geo_factor[4];
} calib[2][2];


/* 
**  Decode the SEPT count rates format:
**
**   4 bit exponent, 10 bit mantissa
*/

static unsigned int dec14(unsigned int c)
{
  unsigned int e, m;
  c &= 0x3fff;
  if (c<0x800)
	  return c;
  e = c >> 10; 
  m = (c & 0x3ff) | 0x400;
  return (m << (e-1)) | (1<<(e-2));
}

/*
** Decompress a spectrum
**  14 bits per counter, 
**  7 bytes for 4 counters,
**  56 bytes for 32 counters.
*/

static void decompress(const unsigned char *ii, sept_science_count_t *oo)
{
  int o = 0;
  int i;
  for (i=0; i<56; i+=7)
    {
      oo[o++] = dec14( ii[i+0]      <<6               | ii[i+1]>>2);
      oo[o++] = dec14((ii[i+1]&0x03)<<12 | ii[i+2]<<4 | ii[i+3]>>4); 
      oo[o++] = dec14((ii[i+3]&0x0f)<<10 | ii[i+4]<<2 | ii[i+5]>>6);
      oo[o++] = dec14((ii[i+5]&0x3f)<<8  | ii[i+6]);
    }
}

/*
** SEPT Science data packet structure
** multi-byte data is Big Endian 
*/

struct septpacket
{
  unsigned char spec[4][56];
  unsigned char single[3];
  unsigned char pdfe[3][4];
  unsigned char intreg[2];
  unsigned char timer[2][3];
  unsigned char status;
  unsigned char state;
  unsigned char sequence;
  unsigned char pulser;
  unsigned char filter;
  unsigned char heater;
  unsigned char temp[2];
  unsigned char spare[5];
  unsigned char checksum;
};


/*
**  SEPT Level 1 data processing
**
**  returns 0 for valid packets
**  returns 1 for wrong APID
*/

static int init_called = 0;

int sept_level1(const struct sept_raw_packet *raw,
		struct sept_science_level1 *l1)
{
  int i,j;

  struct septpacket *l0 = (struct septpacket *) (raw->data);
  struct sept_calibration * cal = &calib[raw->s_c&1][raw->apid&1];

  if ((raw->apid &~1) != 600)
    {
      if (sept_level1_verbose >= 3)
	fprintf(stderr, 
		"sept_level1() called with APID %d, expecting 600 or 601\n",
		raw->apid);
      return 1;
    }
  if (!init_called)
    sept_level1_init();

  l1->version = version | cal->version;

  l1->s_c = raw->s_c & 1;
  l1->inst = raw->apid & 1;
  l1->time = raw->time;

  for (i=0; i<4; i++)
    decompress(l0->spec[i], l1->spec[i]);

  for (i=0; i<2; i++)
    {
      l1->timer[i] 
	= 256*l0->timer[i][0] 
	+     l0->timer[i][1] 
	+     l0->timer[i][2]/256.0
	;
      if (l1->timer[i]<=0)
	l1->timer[i]=0.001;
    }

  for (i=0; i<4; i++)
    for (j=0; j<32; j++)
      l1->spec[i][j] 
	/= cal->geo_factor[i]
	* l1->timer[i>>1] 
	* cal->bin_width[i][j]
	;
	
  /*
  **  Detectors 0,1,4,5 are telescope A 
  **  Detectors 2,3,6,7 are telescope B
  **  Detectors 0,1,2,3 are main channels
  **  Detectors 4,5,6,7 are guard channels
  **  Detectors 0,2,4,6 are electron counters
  **  Detectors 1,3,4,7 are proton counters
  */
  l1->single_counter_id = (l0->status>>4)&7;
  l1->single_counter_rate =
    ( 256*256*l0->single[0]
      +   256*l0->single[1]
      +       l0->single[2] ) 
    / l1->timer[(l1->single_counter_id>>1) & 1];
  

  l1->interrupts = l0->intreg[0]*256 + l0->intreg[1];
  
  for (i=0; i<3; i++)
    for (j=0; j<4; j++)
      l1->pdfe[i][j] = l0->pdfe[i][j];

  l1->pulser = l0->pulser;
  l1->filter = l0->filter;

  /* extract the bits that tell which telescopes are alive */ 
  for (i=0; i<2; i++)
    {
      l1->alive[i] = l0->status >> i & 1;
    }

  /*
  ** Determine the mode from the actual configuration, 
  ** not from state, which is telling the mode for the next packet.
  */
  for (i=0; i<2; i++)
    {
      unsigned int pdfe0 = l0->pdfe[0][2*i+0] >> 5;
      unsigned int pdfe1 = l0->pdfe[0][2*i+1] >> 5;
      unsigned int filter = (l0->filter >> (4*i)) & 0x0f;
      unsigned int pulser = l0->pulser;

      if (! l1->alive[i])
	l1->alive[i] = SEPT_MODE_OFF;
      else if (pdfe0==4 && pdfe1==4 && filter==0x0a && pulser==0)
	l1->alive[i] = SEPT_MODE_SCIENCE;
      else if (pdfe0==5 && pdfe1==5 && filter==0x0f && pulser==0)
	l1->alive[i] = SEPT_MODE_MIPS;
      else if (pulser!=0)
	l1->alive[i] = SEPT_MODE_ITPG;
      else
	l1->alive[i] = SEPT_MODE_UNKNOWN;
    }

  l1->sequence = l0->sequence;
  l1->state = l0->state;
  l1->status = l0->status;
  
  for (i=0; i<2; i++)
    l1->temperature[i] = l0->temp[i] - cal->temp_offset;

  l1->heater = l0->heater;

  if (sept_level1_verbose >= 4)
    {
      fprintf(stderr,
	      "sept_level1: S/C=%c apid=%d time=%.0f\n",
	      "AB"[l1->s_c], l1->inst+600, (double)l1->time);
      if (sept_level1_verbose >= 9)
	sept_level1_print(stderr, l1);
    }
  return 0;
}

/* source SEPT FPGA data sheet, 2.5.4, Table 5 */
static const unsigned int sept_bin_widths[] =
  {
    3,  1,  1,  1,  1,  1,  2,  2,
    2,  2,  3,  3,  3,  4,  4,  4,
    5,  6,  6,  7,  8,  9,  9, 11,
    13, 14, 15, 18, 19, 22, 24, 33 
  };


/* 
**  Read a calibration data file, if available.
**
**  The filename is
**   1. provided in the global variable sept_level1_calibration_file
**   2. the value of the environment variable SEPT_LEVEL1_CALIB
**   3. "sept_level1_calibration.txt"
**  The first of these which is defined is opened for reading.
**  If this fails, a default configuration is used.
**
**  returns 0 for success
**  returns 1 for missing calibration file
**  returns 2 for syntax errors
**  returns 4 for incomplete data
*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>

int sept_level1_verbose = 2;
const char *sept_level1_calibration_file = NULL;

int sept_level1_init()
{
  
  /* source: SEPT FPGA data sheet, 2.5.4, Table 5 */
  static const unsigned int bin_widths[] =
    {
      3,  1,  1,  1,  1,  1,  2,  2,
      2,  2,  3,  3,  3,  4,  4,  4,
      5,  6,  6,  7,  8,  9,  9, 11,
     13, 14, 15, 18, 19, 22, 24, 33 
    };
  const double null_bin_width_MeV = 0.01;
  static const double null_geo_fac[] = { 0.22, 0.26, 0.22, 0.26 };
  const int null_version = 0;
  const int null_temp_offset = 166;

  FILE *f;
  const char *filename = sept_level1_calibration_file;

  if (!filename)
    filename = getenv("SEPT_LEVEL1_CALIB");
  if (!filename)
    filename = "sept_level1_calibration.txt";

  int i,j, sc, inst;
  char line[81];

  int err;
  int gerr = 0;
  
  f = fopen(filename, "r");
  
  if (!f)
    {
      if (sept_level1_verbose >= 1)
	fprintf(stderr, 
		"SEPT Level1: cannot open calibration file %s\n"
		"SEPT Level1: using default calibration version %d\n",
		filename,
		null_version);
      
      for (sc = 0; sc<2; sc++)
	for (inst = 0; inst<2; inst++)
	  {
	    struct sept_calibration *cal = &calib[sc][inst];
	    cal->version = null_version;
	    cal->temp_offset = null_temp_offset;
	    for (i=0; i<4; i++)
	      {
		cal->geo_factor[i] = null_geo_fac[i];
		for (j=0; j<32; j++)
		  cal->bin_width[i][j] = null_bin_width_MeV*bin_widths[j];
	      }
	  }
  
      return 1;
    }
  
  if (sept_level1_verbose >= 2)
    fprintf(stderr, 
	    "SEPT Level1: reading calibration from file %s\n",
	    filename);

  /* make all numbers invalid */
  for (sc = 0; sc<2; sc++)
    for (inst = 0; inst<2; inst++)
      {
	struct sept_calibration *cal = &calib[sc][inst];
	cal->version = 0xfffffff;
	cal->temp_offset = 0xffffffff;
	for (i=0; i<4; i++)
	  {
	    cal->geo_factor[i] = 0.0;
	    for (j=0; j<32; j++)
	      cal->bin_width[i][j] = 0.0;
	  }
      }

  sc = 0;
  inst = 0;
  i = 0;
  j = 32;

  while (fgets(line, sizeof(line), f))
    {
      char *k, *v;
      int lk;

      int iv;
      double dv;
      char cv;
      int viv, vdv;

      k=line;
      while (*k && isspace(*k))
	k++;

      if (!*k || *k=='#' || *k=='\n')
	continue;

      v = k+1;

      while (*v && *v!=':')
	v++;

      if (!*v && sept_level1_verbose>=0)
	{
	  gerr |= 2;
	  fprintf(stderr, 
		  "SEPT Level1: ERROR parsing calibration file: no colon\n"
		  "SEPT Level1: %s\n",
		  line );
	  continue;
	}

      *v = '\0';
      v++;
      
      while (*v && isspace(*v))
	v++;
	
      lk = strlen(k);

      vdv = 1==sscanf(v, "%lf", &dv);
      viv = 1==sscanf(v, "%i", &iv);
      cv = toupper(*v);
      err = 0;

      if (!strcmp(k, "S/C"))
	{
	  sc = cv=='B';
	  err = !sc && cv!='A';
	}
      else if (lk>=4 && !strncmp(k, "Instrument", lk))
	{
	  inst = cv=='E';
	  err = !inst && cv!='N';
	}
      else if (!strcmp(k, "Version"))
	{
	  err = !viv;
	  if (!err)
	    calib[sc][inst].version = iv;
	}
      else if (lk>=4 && !strncmp(k, "TemperatureOffset", lk))
	{
	  err = !viv;
	  if (!err)
	    calib[sc][inst].temp_offset = iv;
	}
      else if (lk>=4 && !strncmp(k, "Telescope", lk))
	{
	  err = !viv || iv<0 || iv>=4;
	  if (!err)
	    i = iv;
	  j=0;
	}
      else if (lk>=3 && !strncmp(k, "BinWidth", lk))
	{
	  err = !vdv || j>=32 || dv<=0;
	  if (!err)
	    calib[sc][inst].bin_width[i][j++] = dv;
	}
      else if (lk>=3 && !strncmp(k, "EnergyPerBin", lk))
	{
	  err = !vdv || j!=0 || dv<=0;
	  for (j=0; j<32; j++)
	    calib[sc][inst].bin_width[i][j] = dv*bin_widths[j];
	}
      else if (lk>=3 && !strncmp(k, "GeometryFactor", lk))
	{
	  err = !vdv || dv<=0;
	  if (!err)
	    calib[sc][inst].geo_factor[i] = dv;
	}
      else
	{
	  gerr |= 2;
	  if (sept_level1_verbose>=0)
	    fprintf(stderr, 
		    "SEPT Level1: ERROR parsing calibration file: illegal key\n"
		    "SEPT Level1: %s\n",
		    line );
	}

      if (err)
	{
	  gerr |= 2;
	  if (sept_level1_verbose>=0)
	    fprintf(stderr, 
		    "SEPT Level1: ERROR parsing calibration file: illegal value\n"
		    "SEPT Level1: %s\n",
		    line );
	}
      
    }

  fclose(f);

  /* check all numbers for invalid */
  for (sc = 0; sc<2; sc++)
    for (inst = 0; inst<2; inst++)
      {
	err = 0;
	struct sept_calibration *cal = &calib[sc][inst];
	if (cal->version == 0xfffffff)
	  {
	    err++;
	    if (sept_level1_verbose >= 1)
	      fprintf(stderr, 
		      "SEPT Level1: %s %s: Version missing\n", 
		      sc ? "B":"A", 
		      inst ? "E":"N/S");
	    cal->version = 0x0000ffff;
	  }
	if (cal->temp_offset == 0xffffffff)
	  {
	    err++;
	    if (sept_level1_verbose >= 1)
	      fprintf(stderr, 
		      "SEPT Level1: %s %s: TemperatureOffset missing\n", 
		      sc ? "B":"A", 
		      inst ? "E":"N/S");
	    cal->temp_offset = null_temp_offset;
	  }
	for (i=0; i<4; i++)
	  {
	    if (cal->geo_factor[i] == 0.0)
	      {
		err++;
		if (sept_level1_verbose >= 1)
		  fprintf(stderr, 
			  "SEPT Level1: %s %s %d GeometryFactor missing\n", 
			  sc ? "B":"A", 
			  inst ? "E":"N/S",
			  i);
		cal->geo_factor[i] = null_geo_fac[i];
	      }
	    for (j=0; j<32; j++)
	      if (cal->bin_width[i][j] == 0.0)
		{
		  err++;
		  if (sept_level1_verbose >= 1)
		    fprintf(stderr, 
			    "SEPT Level1: %s %s %d BinWidth %d missing\n", 
			    sc ? "B":"A", 
			    inst ? "E":"N/S",
			    i, j);
		  cal->bin_width[i][j] = null_bin_width_MeV*bin_widths[j];
		}
	  }
	if (err)
	  {
	    gerr |= 4;
	    cal->version |= 0x8100000;
	    if (sept_level1_verbose >= 0)
	      fprintf(stderr, 
		      "SEPT Level1: ERROR Incomplete Calibration Data, S/C %s, Inst %s\n", 
		      sc ? "B":"A", 
		      inst ? "E":"N/S");
	  }
      }

  init_called ++;

  /*
  ** gerr & 1 : no calibration file
  ** gerr & 2 : syntax error
  ** gerr & 4 : calibration incomplete  
  */
  return gerr;
}

#include <math.h>

int sept_level1_print(FILE *f, const struct sept_science_level1 *l1)
{
  int j,e,s;

  int r=0;
  int rr;
  rr=fprintf(f, 
	     "SEPT science %s %s,  version %x.%04x,  time %.0f\n"
	     "Mode:\t%d\t%d\t"       "Single %d: %.3g\n"
	     "Timer:\t%.3g\t%.3g\t"  "Status: %02x %02x %02x %02x %04x %d\n"
	     "Temp:\t%d\t%d\t"       "PDFE: %02x %02x %02x  %02x %02x %02x  %02x %02x %02x  %02x %02x %02x\n",
	     l1->s_c ? "Behind" : "Ahead", 
	     l1->inst ? "Ecliptic" : "North/South",
	     l1->version >> 16, l1->version & 0xffff,
	     (double)l1->time,
	     l1->alive[0], l1->alive[1],
	     l1->single_counter_id, l1->single_counter_rate,
	     l1->timer[0], l1->timer[1],
	     l1->status, l1->state, l1->pulser, l1->filter, l1->interrupts, l1->sequence,
	     l1->temperature[0], l1->temperature[1],
	     l1->pdfe[0][0], l1->pdfe[1][0], l1->pdfe[2][0], 
	     l1->pdfe[0][1], l1->pdfe[1][1], l1->pdfe[2][1], 
	     l1->pdfe[0][2], l1->pdfe[1][2], l1->pdfe[2][2], 
	     l1->pdfe[0][3], l1->pdfe[1][3], l1->pdfe[2][3] );
  if (rr>=0)
    r+= rr;
  else
    return rr;
  for (e=7; e>0; e--)
    {
      for (s=0; s<2; s++)
	{
	  if (s==0)
	    rr=fprintf(f, e==7 ? " e-|" : "\n   |");
	  else
	    rr=fprintf(f, e==7 ? "    p|" : "     |");
	  if (rr>=0)
	    r+= rr;
	  else
	    return rr;
 
	  for (j=0; j<32; j++)
	    {
	      double c = l1->spec[s][j] + l1->spec[s+2][j];
	      int l = c < 1 ? 0 : 1+(int)floor(log10(c));
	      rr=fprintf(f, l>=e ? "X" : " ");
	      if (rr>=0)
		r+= rr;
	      else
		return rr;
	    }
	}
    }
  rr=fprintf(f, 
	     "\n   +--------------------------------"
	     "     +--------------------------------\n");
  if (rr>=0)
    r+= rr;
  else 
    return rr;
  return r;
}