/* Arguments: input file, output file */

#include        <ctype.h>
#include        <strings.h>
#include        <string.h>
#include        <stdlib.h>
#include        <stdio.h>
#include        <errno.h>
#include        <time.h>
#include        <math.h>
#include        "error.h"

#include        "../../../lib/PPVL-1.10/PPVL.h"
#include        "../../../lib/PPVL-1.10/PPVL_selections.h"

typedef char          *Text;
typedef char unsigned  Byte;
typedef float (*funcptr)( long );

enum Data_order_code
        {
        DATA_ORDER_MSB_FIRST,
        DATA_ORDER_LSB_FIRST
        };

static int Data_order;

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

  PPVL_Parameter *theAggregate, *thisParameter, *theParent, *notesAggregate,
                 *bandbinAggregate, *cubeAggregate, *tempParameter,
		 *isisAggregate, *targetAggregate;
  PPVL_Value     **seq;
  unsigned long  total_scanned;
  static long    Record_bytes, Label_records, File_records, History_record,
                 Image_record, Image_dimensions[3], Pixel_bytes,
		 Suffix_bytes[3], Suffix_multiplier, Data_region[6];
  static Text    Record_type, Data_type_name, Samp_Suffix_type,
                 Band_Suffix_type, Mission_name, Power_state[2],
		 Gain_mode[2], Sampling_mode[3], Back_samp_mode[2];
  Byte           *Cube_data, *Side_plane_data, *Back_plane_data,
                 *Bottom_plane_data;
  static double  Exposure_duration[2], Detector_temp[3], Optics_temp[3];

  Byte		 *trans_short_long(Byte infile[], long dsize);
  void           swap_bytes(Byte indata[], long dsize);
  void           swap_bytes_4(Byte indata[], long dsize);
  static int     Data_Order();
  
  void           fit();

  funcptr        funcs[10];
  char           units[10][64];

  short		 doffset, dataval, testval, line_buf[64];


  int*           new_bg_vals;
  int*           bg_vals;
  int            i, j, temp, rdiff = 0, diff = 0, doswap = 0, Total_lines,
                 Total_bands, line, band, band_suffix, native, badfile,
		 extradatabytes, extradatarecs, head_recs, pad_out,
		 start_band, numinsamp,mwt;
  long           offset, Line_bytes, Samp_suffix_bytes, Band_suffix_bytes,
                 offset_before, offset_after, tmplong, mysum;
  float	         tmpfloat, mymean, mydiff, mydiffsq, myvar, mystd, numstd;
  float          xfd[64],yfd[64],sig[64],siga,sigb,chi2,q,a,b;
  FILE           *infile, *outfile;
  Byte           *line_data, *samp_suffix_data, *band_suffix_data,
                 *out_samp_suffix_data, *out_band_suffix_data;
  char*          version  = "1.0";
  char           *date, *datechomp, *tempc, tmpstr[255], tmptext[64], buff[10],
                 history_line[256];
  time_t         timer;
  int            ivalue,nread;
  double         rvalue;
  const char*    mynull="\0";


  static PPVL_Parameter_Selections
    Parameter_selection_list[] = {
      {"RECORD_TYPE",              PPVL_TYPE_STRING,  1, &Record_type},
      {"RECORD_BYTES",             PPVL_TYPE_INTEGER, 1, &Record_bytes},
      {"FILE_RECORDS",             PPVL_TYPE_INTEGER, 1, &File_records},
      {"LABEL_RECORDS",            PPVL_TYPE_INTEGER, 1, &Label_records},
      {"/^HISTORY",                PPVL_TYPE_INTEGER, 1, &History_record},
      {"/^QUBE",                   PPVL_TYPE_INTEGER, 1, &Image_record},
      {"/QUBE/CORE_ITEMS",         PPVL_TYPE_INTEGER, 3, &Image_dimensions},
      {"/QUBE/CORE_ITEM_TYPE",     PPVL_TYPE_STRING,  1, &Data_type_name},
      {"/QUBE/CORE_ITEM_BYTES",    PPVL_TYPE_INTEGER, 1, &Pixel_bytes},
      {"/QUBE/SUFFIX_ITEMS",       PPVL_TYPE_INTEGER, 3, &Suffix_bytes},
      {"DATA_REGION",              PPVL_TYPE_INTEGER, 6, &Data_region},
      {"/QUBE/SUFFIX_BYTES",       PPVL_TYPE_INTEGER, 1, &Suffix_multiplier},
      {"SAMPLE_SUFFIX_ITEM_TYPE",  PPVL_TYPE_STRING,  1, &Samp_Suffix_type},
      {"BAND_SUFFIX_ITEM_TYPE",    PPVL_TYPE_STRING,  1, &Band_Suffix_type},
      {"MISSION_NAME",             PPVL_TYPE_STRING,  6, &Mission_name},
      {"POWER_STATE_FLAG",         PPVL_TYPE_STRING,  2, &Power_state},
      {"GAIN_MODE_ID",             PPVL_TYPE_STRING,  2, &Gain_mode},
      {"EXPOSURE_DURATION",        PPVL_TYPE_REAL,    2, &Exposure_duration},
      {"DETECTOR_TEMPERATURE",     PPVL_TYPE_REAL,    3, &Detector_temp},
      {"OPTICS_TEMPERATURE",       PPVL_TYPE_REAL,    3, &Optics_temp},
      {"SAMPLING_MODE_ID",         PPVL_TYPE_STRING,  3, &Sampling_mode},


      PPVL_PARAMETER_SELECTIONS_END
      };


  if (argc < 2) {
    printf ("Usage: ir_bg in_cube_filename out_cube_filename\n");
  return;
  }

  if (!(infile = fopen(argv[1], "r"))) {
    printf("ERROR: can't open input file %s\n",argv[1]);
  return;
  }

  /* Read header. */
  theAggregate = PPVL_read_aggregate (infile, 0, &total_scanned);

  /* Get all the header info we need to process the data. */
  i = PPVL_selections(theAggregate,Parameter_selection_list);

  /* Figure out what the background subtraction states are. */
  thisParameter = PPVL_find_parameter(theAggregate,
				      PPVL_SEARCH_FROM_THE_TOP,
				      PPVL_SELECT_NAME,
				      "BACKGROUND_SAMPLING_MODE_ID",
				      &theParent);

  if (PPVL_Type_Is_Sequence(thisParameter->content.value->type)) {
    seq = thisParameter->content.value->data.array;
    for (i=0; i<PPVL_count_all_values(thisParameter->content.value); i++) {
      Back_samp_mode[i] = seq[i]->data.string;
    }
  } else {
    Back_samp_mode[0] = thisParameter->content.value->data.string;
    Back_samp_mode[1] = Back_samp_mode[0];
  }

  /* Open output file. */
  outfile = fopen(argv[2], "wb");

  /* Use PPVL_write_parameter to write the label. */
  PPVL_write_parameter (theAggregate, outfile, 0);

  offset = ftell(outfile);

  /* If the header changed length, adjust for this and update all related
   * keywords.
   * */
  if (offset % Record_bytes == 0) {
    diff = (offset/Record_bytes) - Label_records;
  } else {
    diff = (offset/Record_bytes) + 1 - Label_records;
  }
  if (diff < 0) diff = 0;

  thisParameter = PPVL_find_parameter(theAggregate,
				      PPVL_SEARCH_FROM_THE_TOP,
				      PPVL_SELECT_NAME,
				      "LABEL_RECORDS",
				      &theParent);
  if (thisParameter) thisParameter->content.value->data.integer += diff;
  Label_records += diff;
  thisParameter = PPVL_find_parameter(theAggregate,
				      PPVL_SEARCH_FROM_THE_TOP,
				      PPVL_SELECT_NAME,
				      "FILE_RECORDS",
				      &theParent);
  if (thisParameter) thisParameter->content.value->data.integer += (diff+2);
  thisParameter = PPVL_find_parameter(theAggregate,
				      PPVL_SEARCH_FROM_THE_TOP,
				      PPVL_SELECT_NAME,
				      "^HISTORY",
				      &theParent);
  if (thisParameter) thisParameter->content.value->data.integer += (diff);
  thisParameter = PPVL_find_parameter(theAggregate,
				      PPVL_SEARCH_FROM_THE_TOP,
				      PPVL_SELECT_NAME,
				      "^QUBE",
				      &theParent);
  if (thisParameter) thisParameter->content.value->data.integer += (diff+2);
  Image_record += (diff+2);
  
  rewind(outfile);
  PPVL_write_parameter (theAggregate, outfile, 0);
  offset = ftell(outfile);

  /* pad out to end of header records */
  for (i=0; i < Label_records * Record_bytes - offset; i++) {
    temp = putc(32, outfile);
  }

  offset_before = ftell(outfile);

  /* Read history from input image, if any, write it to the output image. */
  offset = (History_record-1) * Record_bytes;
  fseek(infile, offset, SEEK_SET);

  fgets(history_line, (int)256, infile);
  while (!(strncmp(history_line, "END", 3)  == 0 &&
	   strncmp(history_line, "END_", 4) != 0)) {
    fprintf(outfile, "%s", history_line);
    fgets(history_line, (int)256, infile);
  }

  /* Stick in our history record */
  timer = time(NULL);
  fprintf(outfile, "\015\012");
  fprintf(outfile, "GROUP = IR_BG_ADJUST\015\012");
  fprintf(outfile, "   VERSION = %s\015\012",version);
  date = asctime(localtime(&timer));
  *strchr(date,(char)10) = (char)0;
  fprintf(outfile, "   DATE = \"%s\"\015\012",date);
  fprintf(outfile,
      "   SOFTWARE_DESC = \"IR bg adjust \"\015\012");
  fprintf(outfile, "   INPUT_FILE = \"%s\"\015\012", argv[1]);
  fprintf(outfile, "   OUTPUT_FILE = \"%s\"\015\012", argv[2]);
  fprintf(outfile, "END_GROUP = IR_BG_ADJUST\015\012");
  fprintf(outfile, "END\015\012");
  offset_after = ftell(outfile);
  rdiff = offset_after - offset_before;

  if (rdiff % Record_bytes == 0) {
    head_recs = rdiff/Record_bytes;
    pad_out = 0;
  } else {
    head_recs = rdiff/Record_bytes + 1;
    pad_out = Record_bytes - (rdiff % Record_bytes);
  }

  /* Pad out history record with spaces. */
  if (pad_out > 0) {
    for (i=0; i < pad_out; i++) {
      temp = putc(32, outfile);
    }
  }

  /* If the data doesn't start here, pad out to start of data. */
  offset = ftell(outfile);

  if (offset/Record_bytes+1 < Image_record) {
    rdiff = Image_record - (offset/Record_bytes+1);
    for (i=0; i < Record_bytes*rdiff; i++) {
      temp = putc(32, outfile);
    }
  }

  /* Read data from input image, translate if necessary, write output image. */
  /* swap the bytes if requested */
  /* if there are suffix items that are 2 bytes, transform them to 4 bytes */
  offset = ((Image_record-1) - (diff+2)) * Record_bytes;
  fseek(infile, offset, SEEK_SET);

  Line_bytes = Image_dimensions[0] * Pixel_bytes;
  line_data = (Byte *)calloc(Line_bytes, 1);
  Total_lines = Image_dimensions[2];
  Total_bands = Image_dimensions[1];
  Samp_suffix_bytes = Suffix_bytes[0] * Suffix_multiplier;
  samp_suffix_data = (Byte *)calloc(Samp_suffix_bytes, 1);
  Band_suffix_bytes = Image_dimensions[0] * Suffix_multiplier;
  band_suffix_data = (Byte *)calloc(Band_suffix_bytes, 1);
  bg_vals = (int *)malloc(Total_lines * Total_bands * sizeof(int));
  new_bg_vals = (int *)malloc(Total_lines * Total_bands * sizeof(int));

  /* Go through the cube and pick out all the background values */
  for (line = 0; line < Total_lines; line++) {
    for (band = 0; band < Total_bands; band++) {
      fread(line_data, Line_bytes, 1, infile);

      if (Suffix_bytes[0] > 0) {
        fread(samp_suffix_data, Samp_suffix_bytes, 1, infile);
	if (strcmp(Power_state[1],"ON") != 0) {
	  /* If we are here, the VIS instrument is off, just do IR. */
	  if (strcmp(Samp_Suffix_type, "SUN_INTEGER") == 0) {
	    if (strcmp(Back_samp_mode[0],"ZERO_SUB")!=0 ) {
	      doffset=(short)((samp_suffix_data[2] << 8) + samp_suffix_data[3]);
	    } else {
	      doffset = (short)0;
	    }
	  } else {
	    if (strcmp(Back_samp_mode[0],"ZERO_SUB")!=0 ) {
	      doffset=(short)((samp_suffix_data[1] << 8) + samp_suffix_data[0]);
	    } else {
	      doffset = (short)0;
	    }
	  }
	} else {
	  /* If we are here, the VIS instrument is ON, this code should work
	     whether the IR instrument is on or off.                         */
	  if (strcmp(Samp_Suffix_type, "SUN_INTEGER") == 0) {
	    if (band<96 && line==0 && strcmp(Back_samp_mode[1],"ZERO_SUB")!=0) {
	      doffset=(short)((samp_suffix_data[2] << 8) + samp_suffix_data[3]);
	    } else if (band >= 96 && strcmp(Back_samp_mode[0],"ZERO_SUB")!=0 ) {
	      doffset=(short)((samp_suffix_data[2] << 8) + samp_suffix_data[3]);
	    } else {
	      doffset = (short)0;
	    }
	  } else {
	    if (band<96 && line==0 && strcmp(Back_samp_mode[1],"ZERO_SUB")!=0) {
	      doffset=(short)((samp_suffix_data[1] << 8) + samp_suffix_data[0]);
	    } else if (band >= 96 && strcmp(Back_samp_mode[0],"ZERO_SUB")!=0 ) {
	      doffset=(short)((samp_suffix_data[1] << 8) + samp_suffix_data[0]);
	    } else {
	      doffset = (short)0;
	    }
	  }
	}

        bg_vals[band + line * Total_bands] = doffset;

      } else {
        bg_vals[band + line * Total_bands] = 0;
      }
    }

    /* Skip over backplanes */
    if (Suffix_bytes[1] > 0) {
      for (band_suffix = 0; band_suffix < Suffix_bytes[1]; band_suffix++) {
	nread = fread(band_suffix_data, Band_suffix_bytes, 1, infile);

	/* Back corner if sideplane */
	if (Suffix_bytes[0] > 0) {
	  fread(samp_suffix_data, Samp_suffix_bytes, 1, infile);
	}
      }
    }
  }

  /* Clean up the background values, make averages. */
  numstd = 4.0;  /* Number of standard deviations for outliers */
  if (strcmp(Power_state[1],"ON") != 0) {
    /* Vis instrument is off. */
    start_band = 0;
  } else {
    /* Vis instrument is on. */
    start_band = 96;
  }

  /* For each band, calculate mean and standard deviation of bg points,
     toss out extream points, fit line to remaining points, calculate new
     background values. */
  for (band = start_band; band < Total_bands; band++) {
    mysum = 0;
    myvar = 0.0;
    for (line = 0; line < Total_lines; line++) {
      mysum += bg_vals[band + line * Total_bands];
    }
    mymean = ((float)mysum)/((float)Total_lines);
    for (line = 0; line < Total_lines; line++) {
      mydiff = ((float)bg_vals[band + line * Total_bands]) - mymean;
      mydiffsq = mydiff * mydiff;
      myvar += mydiffsq;
    }
    myvar = myvar/((float)(Total_lines-1));
    mystd = sqrt(myvar);

    numinsamp = -1;
    for (line = 0; line < Total_lines; line++) {
      if (bg_vals[band + line * Total_bands] <= (mymean + numstd * mystd) &&
	  bg_vals[band + line * Total_bands] >= (mymean - numstd * mystd)) {
	numinsamp += 1;
	xfd[numinsamp] = (float)line;
	yfd[numinsamp] = (float)bg_vals[band + line * Total_bands];
      }
    }

    /* Do a linear fit to the data.  Then go back and recalculate new
       values for the data points based on the fit. */
    mwt = 0;
    fit(xfd,yfd,numinsamp,sig,mwt,&a,&b,&siga,&sigb,&chi2,&q);

    /* Set all the new background values in this line to fit values. */
    for (line = 0; line < Total_lines; line++) {
      new_bg_vals[band + line * Total_bands] = (int)((a+b*line)+0.5);
    }
  }

  /* Seek back to the beginning of the data. */
  fseek(infile, offset, SEEK_SET);

  /* Go back through the data, add the background from the sideplane
     back in to the data and then subtract the new background from
     the data.  (IR only) */
  for (line = 0; line < Total_lines; line++) {
    for (band = 0; band < Total_bands; band++) {
      fread(line_data, Line_bytes, 1, infile);

      if (Suffix_bytes[0] > 0) {
        fread(samp_suffix_data, Samp_suffix_bytes, 1, infile);

	if (strcmp(Data_type_name, "SUN_INTEGER") == 0) {
	  for (i=0; i < Line_bytes - 1; i+=2) {
	    dataval =  (short)((line_data[i] << 8) + line_data[i+1]);
	    testval = dataval + bg_vals[band + line * Total_bands];
	    if (band >= start_band &&
		strcmp(Back_samp_mode[0],"ZERO_SUB")!=0 ) {
	      dataval += bg_vals[band + line * Total_bands];
	      dataval -= new_bg_vals[band + line * Total_bands];
	    }
	    line_data[i] = (Byte)(dataval >> 8); 
	    line_data[i+1] = (Byte)(dataval & 0x00FF); 
	  }
	} else {
	  for (i=0; i < Line_bytes - 1; i+=2) {
	    dataval =  (short)((line_data[i+1] << 8) + line_data[i]);
	    testval = dataval + bg_vals[band + line * Total_bands];
	    if (testval >= 4095) {
	      dataval = -32765;
	    } else {
	      if (band >= start_band &&
		  strcmp(Back_samp_mode[0],"ZERO_SUB")!=0 ) {
		dataval += bg_vals[band + line * Total_bands];
		dataval -= new_bg_vals[band + line * Total_bands];
	      }
	    }
	    line_data[i+1] = (Byte)(dataval >> 8); 
	    line_data[i] = (Byte)(dataval & 0x00FF); 
	  }
	}
      
        fwrite(line_data, Line_bytes, 1, outfile);
 
	/* Put the new background value into the side plane. */
	if (band >= start_band && strcmp(Back_samp_mode[0],"ZERO_SUB")!=0 ) {
	  dataval = (short)new_bg_vals[band + line * Total_bands];
	  if (strcmp(Data_type_name, "SUN_INTEGER") == 0) {
	    samp_suffix_data[2] = (Byte)(dataval >> 8); 
	    samp_suffix_data[3] = (Byte)(dataval & 0x00FF); 
	  } else {
	    samp_suffix_data[1] = (Byte)(dataval >> 8); 
	    samp_suffix_data[0] = (Byte)(dataval & 0x00FF); 
	  }
	}
	fwrite(samp_suffix_data, Samp_suffix_bytes, 1, outfile);
      } else {
        fwrite(line_data, Line_bytes, 1, outfile);
      }

    }
    /* backplanes */
    if (Suffix_bytes[1] > 0) {
      for (band_suffix = 0; band_suffix < Suffix_bytes[1]; band_suffix++) {
	nread = fread(band_suffix_data, Band_suffix_bytes, 1, infile);
	fwrite(band_suffix_data, Band_suffix_bytes, 1, outfile);

	if (Suffix_bytes[0] > 0) {
	  fread(samp_suffix_data, Samp_suffix_bytes, 1, infile);
	  fwrite(samp_suffix_data, Samp_suffix_bytes, 1, outfile);
	}
      }
    }
  }

  /* do bottom planes here */
  /* this is currently not implemented */
  
  fclose(infile);

  /* Close output file. */
  fclose(outfile);

}