l1bgen_oci.cpp

Go to the documentation of this file.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <unistd.h>
 
#include <sstream>
#include <fstream>
#include <iomanip>
#include <getopt.h>
#include <libgen.h>
 
#include "nc4utils.h"
#include "l1bgen_oci.h"
 
using namespace std;
using namespace netCDF;
using namespace netCDF::exceptions;
 
#define VERSION "0.045"
 
//    Modification history:
//  Programmer     Organization   Date     Ver   Description of change
//  ----------     ------------   ----     ---   ---------------------
//  Joel Gales     SAIC           04/29/20 0.01  Original development
//  Joel Gales     SAIC           01/13/21 0.02  Add support for SWIR
//  Joel Gales     SAIC           08/11/21 0.03  Add support for handling
//                                               fill values in science data
//  Joel Gales     SAIC           08/12/21 0.04  Add support for hysteresis
//                                               correction
//  Joel Gales     SAIC           09/23/21 0.045 Initialize uninitialized
//                                               variables
 
ofstream tempOut;
 
int main (int argc, char* argv[])
{
 
  cout << "l1bgen_oci " << VERSION << " (" 
       <<  __DATE__ << " " << __TIME__ << ")" << endl;
 
  if ( argc == 1) {
    cout << endl << 
      "l1bgen_oci OCI_L1A_file, cal_LUT_file temp_GEO_LUT_file OCI_L1B_file"
     << endl;
    return 0;
  }
  
  // ********************************* //
  // *** Read calibration LUT file *** //
  // ********************************* //
  NcFile *calLUTfile = new NcFile( argv[2], NcFile::read);
 
  NcGroup gidCommon, gidBlue, gidRed, gidSWIR;
  gidCommon = calLUTfile->getGroup( "common");
  gidBlue = calLUTfile->getGroup( "blue");
  gidRed = calLUTfile->getGroup( "red");
  gidSWIR = calLUTfile->getGroup( "SWIR");
    
  float bwave[NBWAVE];
  float rwave[NRWAVE];
  float swave[NIWAVE];
  float spass[NIWAVE] = {45, 80, 30, 30, 15, 75, 75, 50, 75};
  double K2t[NTIMES];
  float K3T[NTEMPS];
  
  NcVar var;
 
  var = gidCommon.getVar( "blue_wavelength");
  var.getVar( bwave);
 
  var = gidCommon.getVar( "red_wavelength");
  var.getVar( rwave);
 
  var = gidCommon.getVar( "SWIR_wavelength");
  var.getVar( swave);
 
  var = gidCommon.getVar( "K2t");
  var.getVar( K2t);
 
  var = gidCommon.getVar( "K3T");
  var.getVar( K3T);
 
  string tag;
 
  cal_lut_struct blue_lut;
  cal_lut_struct red_lut;
  cal_lut_struct swir_lut;
 
  uint32_t bbanddim, rbanddim, sbanddim;
  
  tag.assign("blue");
  read_oci_cal_lut( calLUTfile, tag, gidBlue, bbanddim, blue_lut);
 
  tag.assign("red");
  read_oci_cal_lut( calLUTfile, tag, gidRed, rbanddim, red_lut);
 
  tag.assign("SWIR");
  read_oci_cal_lut( calLUTfile, tag, gidSWIR, sbanddim, swir_lut);
 
  // Read hysterisis parameters
  float hysttime[9][3];
  float hystamp[9][3];
  
  var = gidSWIR.getVar( "hyst_time_const");
  var.getVar( &hysttime[0][0]);
  var = gidSWIR.getVar( "hyst_amplitude");
  var.getVar( &hystamp[0][0]);
  
  calLUTfile->close();
  
  //  for (size_t i=0; i<NBWAVE; i++) cout << bwave[i] << endl;
 
  
  // ******************************** //
  // *** Read geo LUT (temporary) *** //
  // ******************************** //
  NcFile *geoLUTfile = new NcFile( argv[3], NcFile::read);
  geo_struct geoLUT;
 
  NcGroup gidTime, gidCT, gidRTA_HAM;
 
  gidTime = geoLUTfile->getGroup( "time_params");
  var = gidTime.getVar( "master_clock");
  var.getVar( &geoLUT.master_clock);
  var = gidTime.getVar( "MCE_clock");
  var.getVar( &geoLUT.MCE_clock);
 
  gidCT = geoLUTfile->getGroup( "coord_trans");
  var = gidCT.getVar( "sc_to_tilt");
  var.getVar( &geoLUT.sc_to_tilt);
  var = gidCT.getVar( "tilt_axis");
  var.getVar( &geoLUT.tilt_axis);
  var = gidCT.getVar( "tilt_angles");
  var.getVar( &geoLUT.tilt_angles);
  var = gidCT.getVar( "tilt_to_oci_mech");
  var.getVar( &geoLUT.tilt_to_oci_mech);
  var = gidCT.getVar( "oci_mech_to_oci_opt");
  var.getVar( &geoLUT.oci_mech_to_oci_opt);
 
  gidRTA_HAM = geoLUTfile->getGroup( "RTA_HAM_params");
  var = gidRTA_HAM.getVar( "RTA_axis");
  var.getVar( &geoLUT.rta_axis);
  var = gidRTA_HAM.getVar( "HAM_axis");
  var.getVar( &geoLUT.ham_axis);
  var = gidRTA_HAM.getVar( "HAM_AT_angles");
  var.getVar( &geoLUT.ham_at_angles);
  var = gidRTA_HAM.getVar( "HAM_CT_angles");
  var.getVar( &geoLUT.ham_ct_angles);
  var = gidRTA_HAM.getVar( "RTA_enc_scale");
  var.getVar( &geoLUT.rta_enc_scale);
  var = gidRTA_HAM.getVar( "HAM_enc_scale");
  var.getVar( &geoLUT.ham_enc_scale);
  var = gidRTA_HAM.getVar( "RTA_nadir");
  var.getVar( &geoLUT.rta_nadir);
 
  geoLUTfile->close();
 
  // Open a read data from L1Afile
  NcFile *l1afile = new NcFile( argv[1], NcFile::read);
   
  NcGroup ncGrps[4];
 
  ncGrps[0] = l1afile->getGroup( "scan_line_attributes");
  ncGrps[1] = l1afile->getGroup( "spatial_spectral_modes");
  ncGrps[2] = l1afile->getGroup( "engineering_data");
  ncGrps[3] = l1afile->getGroup( "science_data");
 
  NcGroupAtt att;
 
  // Get date (change this when year and day are added to time field)
  string tstart, tend;
  att = l1afile->getAtt("time_coverage_start");
  att.getValues(tstart);
  cout << tstart << endl;
 
  att = l1afile->getAtt("time_coverage_end");
  att.getValues(tend);
  cout << tend << endl;
 
  uint16_t iyr, imn, idom;
  istringstream iss;
 
  iss.str(tstart.substr(0,4));
  iss >> iyr; iss.clear();
  iss.str(tstart.substr(5,2));
  iss >> imn; iss.clear();
  iss.str(tstart.substr(8,2));
  iss >> idom;
  int32_t jd = jday(iyr, imn, idom);
 
  int32_t iyr32, idy32;
  jdate( jd, &iyr32, &idy32);
 
  // Get numbers of blue and red bands
  NcDim blue_dim = l1afile->getDim("blue_bands");
  uint32_t bbands = blue_dim.getSize();
  NcDim red_dim = l1afile->getDim("red_bands");
  uint32_t rbands = red_dim.getSize();
 
  // Scan time, spin ID and HAM side
  NcDim nscan_dim = l1afile->getDim("number_of_scans");
  uint32_t nscan = nscan_dim.getSize();
  
  double *sstime = new double[nscan];
  var = ncGrps[0].getVar( "scan_start_time");
  var.getVar( sstime);
 
  int32_t *spin = new int32_t[nscan];
  var = ncGrps[0].getVar( "spin_ID");
  var.getVar( spin);
 
  uint8_t *hside = new uint8_t[nscan];
  var = ncGrps[0].getVar( "HAM_side");
  var.getVar( hside);
  
  uint32_t nscan_good=0;
  for (size_t i=0; i<nscan; i++) {
    if ( spin[i] > 0) {
      sstime[nscan_good] = sstime[i];
      hside[nscan_good] = hside[i];
      nscan_good++;
    }
  }
 
  // Check for and fill in missing scan times
 
 
  // ******************************************** //
  // *** Get spatial and spectral aggregation *** //
  // ******************************************** //
  NcDim ntaps_dim = l1afile->getDim("number_of_taps");
  uint32_t ntaps = ntaps_dim.getSize();
  NcDim spatzn_dim = l1afile->getDim("spatial_zones");
  uint32_t spatzn = spatzn_dim.getSize();
 
  int16_t *dtype = new int16_t[spatzn];
  var = ncGrps[1].getVar( "spatial_zone_data_type");
  var.getVar( dtype);
 
  int16_t *lines = new int16_t[spatzn];
  var = ncGrps[1].getVar( "spatial_zone_lines");
  var.getVar( lines);
 
  int16_t *iagg = new int16_t[spatzn];
  var = ncGrps[1].getVar( "spatial_aggregation");
  var.getVar( iagg);
 
  int16_t *bagg = new int16_t[ntaps];
  var = ncGrps[1].getVar( "blue_spectral_mode");
  var.getVar( bagg);
 
  int16_t *ragg = new int16_t[ntaps];
  var = ncGrps[1].getVar( "red_spectral_mode");
  var.getVar( ragg);
 
  
  // ********************************************************************* //
  // *** Get # of EV lines/offset from scan start time to EV mid-time  *** //
  // ********************************************************************* //
  // This will be done by geolocation when integrated into L1B
  int32_t ppr_off;
  double revpsec, secpline;
  int32_t *mspin = new int32_t[nscan];
  int32_t *ot_10us = new int32_t[nscan];
  uint8_t *enc_count = new uint8_t[nscan];
 
  NcDim nenc_dim = l1afile->getDim("encoder_samples");
  uint32_t nenc = nenc_dim.getSize();
 
  float **hamenc = new float *[nscan];
  hamenc[0] = new float[nenc*nscan];
  for (size_t i=1; i<nscan; i++) hamenc[i] = hamenc[i-1] + nenc;
 
  float **rtaenc = new float *[nscan];
  rtaenc[0] = new float[nenc*nscan];
  for (size_t i=1; i<nscan; i++) rtaenc[i] = rtaenc[i-1] + nenc;
 
  read_mce_tlm( l1afile, ncGrps[2], nscan, nenc, ppr_off, revpsec, secpline,
                mspin, ot_10us, enc_count, &hamenc[0], rtaenc);
 
  float clines[32400], slines[4050];
  uint16_t pcdim, psdim;
  int16_t iret;
  double ev_toff, deltc[32400], delts[4050];
  get_ev( secpline, dtype, lines, iagg, pcdim, psdim, ev_toff, clines, slines,
          deltc, delts, iret);
  if ( iret < 0) {
    cout << "No science collect in file: " << argv[1] << endl;
    l1afile->close();
    return 1;
  }
 
  double *evtime = new double[nscan];
  for (size_t i=0; i<nscan_good; i++) evtime[i] = sstime[i] + ev_toff;
 
  size_t ka;
  for (size_t i=0; i<spatzn; i++) {
    if ( dtype[i] != 0 && dtype[i] != 2 && dtype[i] != 10) {
      ka = i;
      break;
    }
  }
 
 
  // *********************************************************** //
  // *** Generate matrices for spectral and gain aggregation *** //
  // *********************************************************** //
  size_t *ia;
  uint32_t iia;
  uint32_t ntb[16];
  ia = new size_t[ntaps];
  
  // Blue bands
  iia=0;
  for (size_t i=0; i<ntaps; i++) {
    if ( bagg[i] > 0) {
      ia[iia] = i;
      iia++;
    }
  }
 
  uint32_t bib=1, bbb=1;
  if ( iia == 0) {
    cout << "All blue taps disabled" << endl;
  } else {
    for (size_t i=0; i<16; i++) ntb[i] = 0;
    for (size_t i=0; i<iia; i++) ntb[ia[i]] = 32 / bagg[ia[i]];
    bib = 0;
    for (size_t i=0; i<16; i++) bib += ntb[i];
  }
  
  // Note: bgmat & rgmat not necessarily contiguous
  float **bamat = new float*[bib];  
  float **bgmat = new float*[512];
 
  if (bib != 1)
    get_agg_mat( ia, iagg[ka], bagg, bib, bbb, ntb, bamat, bgmat);
 
  if (bib != bbands) {
    cout << "Number of blue bands in file: " << argv[1] <<
      " not consistent with spectral aggregation" << endl;
    l1afile->close();
    return 1;
  } else if ( bib < 4) cout << "No blue bands in file: " << argv[1] << endl;
    
  // Red bands
  iia=0;
  for (size_t i=0; i<ntaps; i++) {
    if ( ragg[i] > 0) {
      ia[iia] = i;
      iia++;
    }
  }
 
  uint32_t rib=1, rbb=1;
  if ( iia == 0) {
    cout << "All red taps disabled" << endl;
  } else {
    for (size_t i=0; i<16; i++) ntb[i] = 0;
    for (size_t i=0; i<iia; i++) ntb[ia[i]] = 32 / ragg[ia[i]];
    rib = 0;
    for (size_t i=0; i<16; i++) rib += ntb[i];
  }
  
  float **ramat = new float*[rib];  
  float **rgmat = new float*[512];
 
  if (rib != 1)
    get_agg_mat( ia, iagg[ka], ragg, rib, rbb, ntb, ramat, rgmat);
 
  if (rib != rbands) {
    cout << "Number of red bands in file: " << argv[1] <<
      " not consistent with spectral aggregation" << endl;
    l1afile->close();
    return 1;
  } else if ( rib < 4) cout << "No red bands in file: " << argv[1] << endl;
 
  uint16_t swb = 9;
 
 
  // ********************************* //
  // *** Get dark collect location *** //
  // ********************************* //
  int16_t kd=-1;
  for (size_t i=0; i<spatzn; i++) {
    if ( dtype[i] == 2) {
      kd = (int16_t) i;
    }
  }
  if ( kd == -1) {
    cout << "No dark collect in file: " << argv[1] << endl;
    l1afile->close();
    return 1;
  }
 
  int16_t ldc=0, lds=0;
  for (size_t i=0; i<(size_t) kd; i++) {
    if ( dtype[i] != 0 && dtype [1] != 10) {
      ldc += lines[i] / iagg[i];
      lds += lines[i] / 8;
    }
  }
  int16_t ndc = lines[kd] / iagg[kd];
  int16_t nds = lines[kd] / 8;
 
 
  // *********************************************************************** //
  // *** Generate band gain structs from LUTs, date/time & gain matrices *** //
  // *********************************************************************** //
  gains_struct bgains;
  gains_struct rgains;
  gains_struct sgains;
  /*
  for (size_t i=0; i<64; i++)
    for (size_t j=0; j<512; j++)
      cout << "0: " << i << " " << j << " " << bgmat[j][i] << endl;
  */
  // Note: bgmat & rgmat not necessarily contiguous
  // That's why we pass the location of the array of pointers
  if ( bib >= 4)
    make_oci_gains( bib, bbanddim, iyr, idom, evtime[0], K2t, blue_lut,
                    &bgmat[0], bgains);
 
  if ( rib >= 4)
    make_oci_gains( rib, rbanddim, iyr, idom, evtime[0], K2t, red_lut,
                    &rgmat[0], rgains);
 
  float **sgmat = new float*[swb];
  for (size_t i=0; i<swb; i++) {
    sgmat[i] = new float[swb];
    for (size_t j=0; j<swb; j++) {
      if (i == j) sgmat[i][j] = 1.0; else sgmat[i][j] = 0.0;
    }
  }
  make_oci_gains( swb, swb, iyr, idom, evtime[0], K2t, swir_lut,
                  &sgmat[0], sgains);
 
  // Read selected temperature fields and interpolate to scan times
  float **caltemps = new float *[nscan_good];
  caltemps[0] = new float[30*nscan_good];
  for (size_t i=1; i<nscan_good; i++) caltemps[i] = caltemps[i-1] + 30;
 
  // get_oci_cal_temps JMG not yet defined
  for (size_t i=0; i<nscan_good; i++)
    for (size_t j=1; j<30; j++)
      caltemps[i][j] = 0.0;
  
  // Read dark collects from science data arrays
  vector<size_t> start, count;
  start.push_back(0);
  start.push_back(0);
  start.push_back(ldc);
 
  size_t dims[3];
  dims[0] = nscan_good; dims[1] = bib; dims[2] = ndc;
  uint16_t ***bdark = make3dT<uint16_t>(dims);
  dims[0] = nscan_good; dims[1] = rib; dims[2] = ndc;
  uint16_t ***rdark = make3dT<uint16_t>(dims);
 
  uint16_t bfill;
  uint16_t rfill;
  
  if ( bib > 4) {
    count.push_back(nscan_good);
    count.push_back(bib);
    count.push_back(ndc);
    
    var = ncGrps[3].getVar( "sci_blue");
    var.getVar( start, count, &bdark[0][0][0]);
 
    var.getAtt("_FillValue").getValues(&bfill);
  }
 
  if ( rib > 4) {
    count.clear();
    count.push_back(nscan_good);
    count.push_back(rib);
    count.push_back(ndc);
    
    var = ncGrps[3].getVar( "sci_red");
    var.getVar( start, count, &rdark[0][0][0]);
 
    var.getAtt("_FillValue").getValues(&rfill);
  }
 
  dims[0] = nscan_good; dims[1] = swb; dims[2] = nds;
  uint32_t ***sdark = make3dT<uint32_t>(dims);
 
  uint32_t sfill;
  
  start.clear();
  start.push_back(0);
  start.push_back(0);
  start.push_back(lds);
 
  count.clear();
  count.push_back(nscan_good);
  count.push_back(swb);
  count.push_back(nds);
 
  var = ncGrps[3].getVar( "sci_SWIR");
  var.getVar( start, count, &sdark[0][0][0]);
 
  var.getAtt("_FillValue").getValues(&sfill);
 
  uint32_t fill32;
  
  static l1bFile outfile;
 
  // nscan_good: Number of scans in L1A file
  // pcdim:      Number of CCD band pixels 
  // bbb:        Number of blue bands]
  // rbb:        Number of red bands
  // psdim:      Number of SWIR band pixels
  // swb:        Number of SWIR bands
  outfile.createl1b( (char *) argv[4], nscan_good, pcdim, bbb, rbb, psdim, swb);
 
  // number of scans of dark data to average; will make this an input parameter
  uint16_t ndsc = 1;
  //number of dark pixels to skip; will make this an input parameter
  uint16_t nskp = 0;
  // adjust SWIR skip factor in case of CCD aggregation LT 8
  //  uint32_t nsskp = (nskp * iagg[kd]) / 8;
 
  uint32_t bicount[3] = {1,bib,(uint32_t) pcdim};
  uint32_t ricount[3] = {1,rib,(uint32_t) pcdim};
  uint32_t sicount[3] = {1,swb,(uint32_t) psdim};
  uint32_t bcount[3] = {bbb,1,(uint32_t) pcdim};
  uint32_t rcount[3] = {rbb,1,(uint32_t) pcdim};
  uint32_t scount[3] = {swb,1,(uint32_t) psdim};
 
  // Calibrated data variables
  float **bdn = new float *[bib];
  bdn[0] = new float[pcdim*bib];
  for (size_t i=1; i<bib; i++) bdn[i] = bdn[i-1] + pcdim;
 
  float **rdn = new float *[rib];
  rdn[0] = new float[pcdim*rib];
  for (size_t i=1; i<rib; i++) rdn[i] = rdn[i-1] + pcdim;
 
  float **sdn = new float *[swb];
  sdn[0] = new float[psdim*swb];
  for (size_t i=1; i<swb; i++) sdn[i] = sdn[i-1] + psdim;
 
  float **bcal = new float *[bib];
  bcal[0] = new float[pcdim*bib];
  for (size_t i=1; i<bib; i++) bcal[i] = bcal[i-1] + pcdim;
 
  float **rcal = new float *[rib];
  rcal[0] = new float[pcdim*rib];
  for (size_t i=1; i<rib; i++) rcal[i] = rcal[i-1] + pcdim;
 
  float **scal = new float *[swb];
  scal[0] = new float[psdim*swb];
  for (size_t i=1; i<swb; i++) scal[i] = scal[i-1] + psdim;
 
  double *thetap = new double[pcdim];
  double *thetas = new double[psdim];
 
  float **pview = new float *[pcdim];
  pview[0] = new float[3*pcdim];
  for (size_t i=1; i<pcdim; i++) pview[i] = pview[i-1] + 3;
 
  float **sview = new float *[psdim];
  sview[0] = new float[3*psdim];
  for (size_t i=1; i<psdim; i++) sview[i] = sview[i-1] + 3;
 
  //  ;  pview(3,pdim) R*4   O  View vectors in the instrument frame
  //;  theta(pdim)   R*4     O  Scan angles for science pixels
  
  // Read, calibrate and write science data
  for (size_t iscn=0; iscn<nscan_good; iscn++) {
 
    if ((iscn % 50) == 0) cout << "Calibrating scan: " << iscn << endl;
 
    // Check for valid mirror side
    if ( hside[iscn] == 0 || hside[iscn] == 1) {
 
      // Write evtime
      start.clear();
      start.push_back(iscn);
 
      count.clear();
      count.push_back(1);
 
      var = outfile.ncGrps[1].getVar( "time");
      var.putVar( start, count, &evtime[iscn]);
 
      var = outfile.ncGrps[1].getVar( "HAM_side");
      var.putVar( start, count, &hside[iscn]);
 
      // Get scan angle
      // This will be done by geolocation when integrated into L1B
      get_oci_vecs( nscan, pcdim, geoLUT, ev_toff, spin[iscn], clines, deltc,
                    revpsec, ppr_off, mspin, ot_10us, enc_count, &hamenc[0],
                    &rtaenc[0], pview, thetap, iret);
 
      tempOut.open ("pview.bin", ios::out | ios::trunc | ios::binary);
      tempOut.write((char *) &pview[0][0], sizeof(float)*3*pcdim);
      tempOut.close();
 
      tempOut.open ("thetap.bin", ios::out | ios::trunc | ios::binary);
      tempOut.write((char *) &thetap[0], sizeof(double)*pcdim);
      tempOut.close();
 
      get_oci_vecs( nscan, psdim, geoLUT, ev_toff, spin[iscn], slines, delts,
                    revpsec, ppr_off, mspin, ot_10us, enc_count, &hamenc[0],
                    &rtaenc[0], sview, thetas, iret);
 
      tempOut.open ("sview.bin", ios::out | ios::trunc | ios::binary);
      tempOut.write((char *) &sview[0][0], sizeof(float)*3*psdim);
      tempOut.close();
 
      tempOut.open ("thetas.bin", ios::out | ios::trunc | ios::binary);
      tempOut.write((char *) &thetas[0], sizeof(double)*psdim);
      tempOut.close();
      
      //  Blue bands
      if (bib >= 4) {
 
        start.clear();
        start.push_back(iscn);
        start.push_back(0);
        start.push_back(0);
 
        count.clear();
        count.push_back(bicount[0]);
        count.push_back(bicount[1]);
        count.push_back(bicount[2]);
 
        uint16_t **bsci = new uint16_t *[bib];
        bsci[0] = new uint16_t[pcdim*bib];
        for (size_t i=1; i<bib; i++) bsci[i] = bsci[i-1] + pcdim;
 
        var = ncGrps[3].getVar( "sci_blue");
        var.getVar( start, count, bsci[0]);
 
        tempOut.open ("bsci.bin", ios::out | ios::trunc | ios::binary);
        tempOut.write((char *) bsci[0], sizeof(uint16_t)*pcdim*bib);
        tempOut.close();
      
        // Compute dark offset, correct data, and apply absolute and
        // temporal gain and temperature correction
        float *bdc = new float[bib];
 
        fill32 = bfill;
        int16_t iret;
        get_oci_dark<uint16_t>( iscn, nscan_good, hside, ndsc, nskp,
                                iagg[ka], iagg[kd], ntaps, bagg, fill32,
                                ndc, bdark, bib, bdc, iret);
 
        float *k3 = new float[bib];
        get_oci_temp_corr( bib, bgains, K3T, caltemps[iscn], nscan_good, k3);
        for (size_t j=0; j<bib; j++) {
          for (size_t k=0; k<pcdim; k++) {
 
            // Handle fill value
            if (bsci[j][k] == bfill) {
              bdn[j][k] = -999;
              bcal[j][k] = -999;
              continue;
            }
 
            // Need to save dn for linearity correction
            bdn[j][k] = bsci[j][k] - bdc[j];
            bcal[j][k] = k3[j] * bgains.K1K2[j][hside[iscn]] * bdn[j][k];
          }
        }
 
        delete [] k3;
        delete [] bdc;
        
        // Compute and apply RVS and linearity
        //k4 = fltarr(pdim,nib)
        float **k4 = new float *[bib];
        k4[0] = new float[pcdim*bib];
        for (size_t i=1; i<bib; i++) k4[i] = k4[i-1] + pcdim;
 
        get_oci_rvs_corr( bib, pcdim, hside[iscn], bgains, thetap, k4);
 
        float **k5 = new float *[bib];
        k5[0] = new float[pcdim*bib];
        for (size_t i=1; i<bib; i++) k5[i] = k5[i-1] + pcdim;
        
        get_oci_lin_corr( bib, pcdim, bgains, K3T, caltemps[iscn], bdn, k5);
 
        for (size_t j=0; j<bib; j++) {
          for (size_t k=0; k<pcdim; k++) {
            bcal[j][k] *= k4[j][k] * k5[j][k];
          }
        }
 
        delete [] k4[0];
        delete [] k4;
        delete [] k5[0];
        delete [] k5;
 
        float **bcalb = new float *[bib];
        bcalb[0] = new float[pcdim*bib];
        for (size_t i=1; i<bib; i++) bcalb[i] = bcalb[i-1] + pcdim;
 
        // Aggregate to L1B bands
        //bcalb = transpose(bamat#transpose(bcal))
        for (size_t j=0; j<bib; j++) {
          for (size_t k=0; k<pcdim; k++) {
            float sum = 0.0;
            for (size_t l=0; l<bib; l++)
              if (bcal[l][k] != -999) sum += bamat[j][l]*bcal[l][k];
            bcalb[j][k] = sum;
          }
        }
 
        start.clear();
        start.push_back(0);
        start.push_back(iscn);
        start.push_back(0);
 
        count.clear();
        count.push_back(bcount[0]);
        count.push_back(bcount[1]);
        count.push_back(bcount[2]);
  
        // Output to L1B file
        var = outfile.ncGrps[4].getVar( "Lt_blue");
        var.putVar( start, count, &bcalb[0][0]);
        //        BCALB           FLOAT     = Array[1306, 64]
 
        delete [] bcalb[0];
        delete [] bcalb;
        delete [] bsci[0];
        delete [] bsci;
      } // End blue
 
 
      //  Red bands
      if (rib >= 4) {
 
        start.clear();
        start.push_back(iscn);
        start.push_back(0);
        start.push_back(0);
 
        count.clear();
        count.push_back(ricount[0]);
        count.push_back(ricount[1]);
        count.push_back(ricount[2]);
 
        uint16_t **rsci = new uint16_t *[rib];
        rsci[0] = new uint16_t[pcdim*rib];
        for (size_t i=1; i<rib; i++) rsci[i] = rsci[i-1] + pcdim;
 
        var = ncGrps[3].getVar( "sci_red");
        var.getVar( start, count, rsci[0]);
 
        tempOut.open ("rsci.bin", ios::out | ios::trunc | ios::binary);
        tempOut.write((char *) rsci[0], sizeof(uint16_t)*pcdim*rib);
        tempOut.close();
      
        // Compute dark offset, correct data, and apply absolute and
        // temporal gain and temperature correction
        float *rdc = new float[rib];
        fill32 = rfill;
        int16_t iret;
        get_oci_dark<uint16_t>( iscn, nscan_good, hside, ndsc, nskp,
                                iagg[ka], iagg[kd], ntaps, ragg, fill32,
                                ndc, rdark, rib, rdc, iret);
 
        tempOut.open ("rdc.bin", ios::out | ios::trunc | ios::binary);
        tempOut.write((char *) rdc, sizeof(float)*rib);
        tempOut.close();
        float *k3 = new float[rib];
        get_oci_temp_corr( rib, rgains, K3T, caltemps[iscn], nscan_good, k3);
        tempOut.open ("k3.bin", ios::out | ios::trunc | ios::binary);
        tempOut.write((char *) k3, sizeof(float)*rib);
        tempOut.close();
        for (size_t j=0; j<rib; j++) {
          for (size_t k=0; k<pcdim; k++) {
 
            // Handle fill value
            if (rsci[j][k] == rfill) {
              rdn[j][k] = -999;
              rcal[j][k] = -999;
              continue;
            }
            
            // Need to save dn for linearity correction
            rdn[j][k] = rsci[j][k] - rdc[j];
            rcal[j][k] = k3[j] * rgains.K1K2[j][hside[iscn]] * rdn[j][k];
          }
        }
        tempOut.open ("rdn.bin", ios::out | ios::trunc | ios::binary);
        tempOut.write((char *) rdn[0], sizeof(float)*pcdim*rib);
        tempOut.close();
 
        delete [] k3;
        delete [] rdc;
        
        // Compute and apply RVS and linearity
        //k4 = fltarr(pdim,nib)
        float **k4 = new float *[rib];
        k4[0] = new float[pcdim*rib];
        for (size_t i=1; i<rib; i++) k4[i] = k4[i-1] + pcdim;
        get_oci_rvs_corr( rib, pcdim, hside[iscn], rgains, thetap, k4);
        tempOut.open ("k4.bin", ios::out | ios::trunc | ios::binary);
        tempOut.write((char *) k4[0], sizeof(float)*pcdim*rib);
        tempOut.close();
 
        float **k5 = new float *[rib];
        k5[0] = new float[pcdim*rib];
        for (size_t i=1; i<rib; i++) k5[i] = k5[i-1] + pcdim;
        
        get_oci_lin_corr( rib, pcdim, rgains, K3T, caltemps[iscn], rdn, k5);
        tempOut.open ("k5.bin", ios::out | ios::trunc | ios::binary);
        tempOut.write((char *) k5[0], sizeof(float)*pcdim*rib);
        tempOut.close();
        for (size_t j=0; j<rib; j++) {
          for (size_t k=0; k<pcdim; k++) {
            rcal[j][k] *= k4[j][k] * k5[j][k];
          }
        }
 
        tempOut.open ("rcal.bin", ios::out | ios::trunc | ios::binary);
        tempOut.write((char *) rcal[0], sizeof(float)*pcdim*rib);
        tempOut.close();
        
        delete [] k4[0];
        delete [] k4;
        delete [] k5[0];
        delete [] k5;
 
        float **rcalb = new float *[rib];
        rcalb[0] = new float[pcdim*rib];
        for (size_t i=1; i<rib; i++) rcalb[i] = rcalb[i-1] + pcdim;
 
        // Aggregate to L1B bands
        for (size_t j=0; j<rib; j++) {
          for (size_t k=0; k<pcdim; k++) {
            float sum = 0.0;
            for (size_t l=0; l<rib; l++)
              if (rcal[l][k] != -999) sum += ramat[j][l]*rcal[l][k];
            rcalb[j][k] = sum;
          }
        }
 
        tempOut.open ("rcalb.bin", ios::out | ios::trunc | ios::binary);
        tempOut.write((char *) rcalb[0], sizeof(float)*pcdim*rib);
        tempOut.close();
 
        start.clear();
        start.push_back(0);
        start.push_back(iscn);
        start.push_back(0);
 
        count.clear();
        count.push_back(rcount[0]);
        count.push_back(rcount[1]);
        count.push_back(rcount[2]);
  
        // Output to L1B file
        var = outfile.ncGrps[4].getVar( "Lt_red");
        var.putVar( start, count, &rcalb[0][0]);
 
        delete [] rcalb[0];
        delete [] rcalb;
        delete [] rsci[0];
        delete [] rsci;
      } // End red
 
      //  SWIR bands
      start.clear();
      start.push_back(iscn);
      start.push_back(0);
      start.push_back(0);
 
      count.clear();
      count.push_back(sicount[0]);
      count.push_back(sicount[1]);
      count.push_back(sicount[2]);
 
      uint32_t **ssci = new uint32_t *[swb];
      ssci[0] = new uint32_t[psdim*swb];
      for (size_t i=1; i<swb; i++) ssci[i] = ssci[i-1] + psdim;
 
      var = ncGrps[3].getVar( "sci_SWIR");
      var.getVar( start, count, ssci[0]);
 
      // Compute dark offset, correct data, and apply absolute and
      // temporal gain and temperature correction
      float *sdc = new float[swb];
      //      int16_t sagg = 8;
      int16_t sagg = 1;
      int16_t iret;
      get_oci_dark<uint32_t>( iscn, nscan_good, hside, ndsc, nskp,
                              1, 1, 1, &sagg, sfill, nds, sdark,
                              swb, sdc, iret);
 
      float *k3 = new float[swb];
      if ( iret != -1) {
        get_oci_temp_corr( swb, sgains, K3T, caltemps[iscn], nscan_good, k3);
 
        for (size_t j=0; j<swb; j++) {
          for (size_t k=0; k<psdim; k++) {
 
            // Handle fill value
            if (ssci[j][k] == sfill) {
              sdn[j][k] = -999;
              scal[j][k] = -999;
              continue;
            }
          
            // Need to save dn for linearity correction
            sdn[j][k] = ssci[j][k] - sdc[j];
 
            // Hysteresis correction
            float hc_prev[3];
            float hc[3]={0,0,0};
            float hyst = 0.0;
            for (size_t l=0; l<3; l++) {
              // Compute exponential decay constants
              float e = exp(-1.0/hysttime[j][l]);
 
              if ( k > 0) {
                hc[l] = hc_prev[l]*e + sdn[j][k-1]*hystamp[j][l];
                hyst += hc[l];
              }
              hc_prev[l] = hc[l];
            } // l-loop
 
            scal[j][k] =
              k3[j] * sgains.K1K2[j][hside[iscn]] * (sdn[j][k] - hyst);
          } // k-loop
        }
      } // iret != -1
      
      delete [] k3;
      delete [] sdc;
              
      // Compute and apply RVS and linearity
      float **k4 = new float *[swb];
      k4[0] = new float[psdim*swb];
      for (size_t i=1; i<swb; i++) k4[i] = k4[i-1] + psdim;
      if ( iret != -1)
        get_oci_rvs_corr( swb, psdim, hside[iscn], sgains, thetap, k4);
 
      float **k5 = new float *[swb];
      k5[0] = new float[psdim*swb];
      for (size_t i=1; i<swb; i++) k5[i] = k5[i-1] + psdim;
 
      if ( iret != -1)
        get_oci_lin_corr( swb, psdim, sgains, K3T, caltemps[iscn], sdn, k5);
 
      for (size_t j=0; j<swb; j++) {
        for (size_t k=0; k<psdim; k++) {
          if (scal[j][k] != -999)
            scal[j][k] *= k4[j][k] * k5[j][k];
        }
      }
        
      delete [] k4[0];
      delete [] k4;
      delete [] k5[0];
      delete [] k5;
 
      start.clear();
      start.push_back(0);
      start.push_back(iscn);
      start.push_back(0);
 
      count.clear();
      count.push_back(scount[0]);
      count.push_back(scount[1]);
      count.push_back(scount[2]);
  
      // Output to L1B file
      var = outfile.ncGrps[4].getVar( "Lt_SWIR");
      var.putVar( start, count, &scal[0][0]);
 
      delete [] ssci[0];
      delete [] ssci;
 
    } else {
      cout << "No mirror side index for scan: " << iscn << endl;
    } // Check for valid mirror side
  } // Scan loop
 
  
  // BAMAT           FLOAT     = Array[64, 64]
  // BGMAT           FLOAT     = Array[64, 512]
  // Write spectral band information
  // Calculate band centers for aggregated hyperspectral bands
  // b1bwave = bamat#bgmat#bwave
  if (bib >= 4) {
    // bgmat#bwave
    float *b1bwave = new float[bib];
    float *sum = new float[bib];
    for (size_t i=0; i<bib; i++) {
      sum[i] = 0.0;
      for (size_t j=0; j<512; j++) {
        sum[i] += bwave[j]*bgmat[j][i];
      }
    }
 
    // bamat#sum
    for (size_t i=0; i<bib; i++) {
      b1bwave[i] = 0.0;
      for (size_t j=0; j<bib; j++) {
        b1bwave[i] += bamat[i][j]*sum[j];
      }
    }
    
    start.clear();
    start.push_back(0);
 
    count.clear();
    count.push_back(bbands);
    var = outfile.ncGrps[0].getVar( "blue_wavelength");
    var.putVar( start, count, b1bwave);
 
    delete [] b1bwave;
    delete [] sum;
  }
 
  if (rib >= 4) {
    float *r1bwave = new float[rib];
    float *sum = new float[rib];
    for (size_t i=0; i<rib; i++) {
      sum[i] = 0.0;
      for (size_t j=0; j<512; j++) {
        sum[i] += rwave[j]*rgmat[j][i];
      }
    }
 
    for (size_t i=0; i<rib; i++) {
      r1bwave[i] = 0.0;
      for (size_t j=0; j<rib; j++) {
        r1bwave[i] += ramat[i][j]*sum[j];
      }
    }
    
    start.clear();
    start.push_back(0);
 
    count.clear();
    count.push_back(rbands);
    var = outfile.ncGrps[0].getVar( "red_wavelength");
    var.putVar( start, count, r1bwave);
 
    delete [] r1bwave;
    delete [] sum;
  }
 
  // SWIR wavelengths/bandpass
  start.clear();
  start.push_back(0);
  count.clear();
  count.push_back(NIWAVE);
  
  var = outfile.ncGrps[0].getVar( "SWIR_wavelength");
  var.putVar( start, count, swave);
  var = outfile.ncGrps[0].getVar( "SWIR_bandpass");
  var.putVar( start, count, spass);
 
  
  string l1b_filename;
  l1b_filename.assign(argv[4]);
  outfile.write_granule_metadata( tstart, tend, l1b_filename);
 
  outfile.close();
  
  delete [] sstime;
  delete [] spin;
  delete [] hside;
  delete [] dtype;
  delete [] lines;
  delete [] iagg;
  delete [] bagg;
  delete [] ragg;
  delete [] mspin;
  delete [] ot_10us;
  delete [] enc_count;
  delete [] sgmat;
  delete [] thetap;
  delete [] thetas;
  delete [] ia;
  delete [] evtime;
  
  delete [] hamenc[0];
  delete [] hamenc;
  delete [] rtaenc[0];
  delete [] rtaenc;
  delete [] caltemps[0];
  delete [] caltemps;
  delete [] bdn[0];
  delete [] bdn;
  delete [] rdn[0];
  delete [] rdn;
  delete [] sdn[0];
  delete [] sdn;  
  delete [] bcal[0];
  delete [] bcal;
  delete [] rcal[0];
  delete [] rcal;
  delete [] scal[0];
  delete [] scal;  
  delete [] pview[0];
  delete [] pview;
  delete [] sview[0];
  delete [] sview;  
 
  delete [] blue_lut.K1[0];
  delete [] blue_lut.K1;
  delete [] blue_lut.K5[0];
  delete [] blue_lut.K5;
 
  if (bgains.K1K2 != NULL) delete [] bgains.K1K2[0];
  if (bgains.K1K2 != NULL) delete [] bgains.K1K2;
  if (bgains.K5 != NULL) delete [] bgains.K5[0];
  if (bgains.K5 != NULL) delete [] bgains.K5;
 
  // delete [] arrT[0][0]; delete [] arrT[0]; delete [] arrT;
 
  delete [] red_lut.K1[0];
  delete [] red_lut.K1;
  delete [] red_lut.K5[0];
  delete [] red_lut.K5;
 
  if (rgains.K1K2 != NULL) delete [] rgains.K1K2[0];
  if (rgains.K1K2 != NULL) delete [] rgains.K1K2;
  if (rgains.K5 != NULL) delete [] rgains.K5[0];
  if (rgains.K5 != NULL) delete [] rgains.K5;
  
  delete [] swir_lut.K1[0];
  delete [] swir_lut.K1;
  delete [] swir_lut.K5[0];
  delete [] swir_lut.K5;
 
  delete [] sgains.K1K2[0];
  delete [] sgains.K1K2;
  delete [] sgains.K5[0];
  delete [] sgains.K5;
 
  // Add delete for dark arrays
  
  return 0;
}
 
int read_mce_tlm( NcFile *l1afile, NcGroup egid,
                  uint32_t nscan, uint32_t nenc, int32_t& ppr_off,
                  double& revpsec, double&secpline, int32_t *mspin,
                  int32_t *ot_10us, uint8_t *enc_count,
                  float **hamenc, float **rtaenc) {
 
  NcDim mce_dim = l1afile->getDim("MCE_block");
  uint32_t mce_blk = mce_dim.getSize();
  NcDim ddc_dim = l1afile->getDim("DDC_tlm");
  uint32_t nddc = ddc_dim.getSize();
  NcDim tlm_dim = l1afile->getDim("tlm_packets");
  uint32_t ntlm = tlm_dim.getSize();
 
  uint8_t **mtlm = new uint8_t *[nscan];
  mtlm[0] = new uint8_t[mce_blk*nscan];
  for (size_t i=1; i<nscan; i++) mtlm[i] = mtlm[i-1] + mce_blk;
 
  int16_t **enc = new int16_t *[nscan];
  enc[0] = new int16_t[4*nenc*nscan];
  for (size_t i=1; i<nscan; i++) enc[i] = enc[i-1] + 4*nenc;
  
  uint8_t **ddctlm = new uint8_t *[ntlm];
  ddctlm[0] = new uint8_t[nddc*ntlm];
  for (size_t i=1; i<ntlm; i++) ddctlm[i] = ddctlm[i-1] + nddc;
 
  NcVar var;
  vector<size_t> start, count;
  start.push_back(0);
  start.push_back(0);
  count.push_back(1);
 
  /*
  // Read MCE telemetry byte-by-byte to convert from int8_t to uint8_t
  int8_t *buf = new int8_t[mce_blk];
  var = egid.getVar( "MCE_telemetry");
  count.push_back(mce_blk);
  for (size_t i=0; i<nscan; i++) {
    start[0] = i;
    var.getVar( start, count, buf);
    for (size_t j=0; j<mce_blk; j++) {
      if ( buf[j] & 0x80 == 0)
        mtlm[i][j] = buf[j]; else mtlm[i][j] = 256 + buf[j];
      //      cout << i << " " << j << " " << (int) mtlm[i][j] << endl;
    }
  }
  delete [] buf;
  */
  var = egid.getVar( "MCE_telemetry");
  count.push_back(mce_blk);
  var.getVar( &mtlm[0][0]);
  
  var = egid.getVar( "MCE_encoder_data");
  count.pop_back();
  count.push_back(4*nenc);
  var.getVar( &enc[0][0]);
    
  var = egid.getVar( "MCE_spin_ID");
  var.getVar( mspin);
 
  var = egid.getVar( "DDC_telemetry");
  var.getVar( &ddctlm[0][0]);
 
  int32_t max_enc_cts = 131072; // 2^17
  double clock[2] = {1.36e8,1.0e8};
 
  // Get ref_pulse_divider and compute commanded rotation rate
  uint32_t ui32;
  uint32_t ref_pulse_div[2];
  memcpy( &ui32, &mtlm[0][0], 4);
  ref_pulse_div[0] = SWAP_4( ui32) % 16777216; // 2^24
  memcpy( &ui32, &mtlm[0][4], 4);
  ref_pulse_div[1] = SWAP_4( ui32) % 16777216; // 2^24
 
  int32_t ref_pulse_sel = mtlm[0][9] / 128;
  
  revpsec = clock[ref_pulse_sel] / 2 / max_enc_cts /
    (ref_pulse_div[ref_pulse_sel]/256.0 + 1);
 
  // Get PPR offset and on-time_10us
  memcpy( &ui32, &mtlm[0][8], 4);
  ppr_off = SWAP_4( ui32) % max_enc_cts;
  for (size_t i=0; i<nscan; i++) {
      memcpy( &ui32, &mtlm[i][368], 4);
      ot_10us[i] = SWAP_4( ui32) % 4;
  }
 
  // Get TDI time and compute time increment per line
  uint16_t ui16;
  memcpy( &ui16, &ddctlm[0][346], 2);
  int32_t tditime = SWAP_2( ui16);
  secpline = (tditime+1) / clock[0]; 
 
  // Get valid encoder count, HAM and RTA encoder data
  for (size_t i=0; i<nscan; i++) enc_count[i] = mtlm[i][475];
  float enc_s = 81.0 / 2560;
  for (size_t i=0; i<nscan; i++) {
    for (size_t j=0; j<nenc; j++) {
      hamenc[i][j] = enc[i][4*j+0] * enc_s;
      rtaenc[i][j] = enc[i][4*j+1] * enc_s;
    }
  }
 
  delete[] mtlm[0];
  delete[] mtlm;
 
  delete [] enc[0];
  delete [] enc;
 
  delete [] ddctlm[0];
  delete [] ddctlm;
 
  return 0;
}
 
int get_ev( double secpline, int16_t *dtype, int16_t *lines, int16_t *iagg,
            uint16_t& pcdim, uint16_t& psdim, double& ev_toff,
            float *clines, float *slines, double *deltc, double *delts,
            int16_t &iret) {
 
  // Find end of no-data zone
  int16_t iz=0, line0=0, linen;
  iret = -1;
  while ( dtype[iz] == 0) {
    line0 += lines[iz];
    iz++;
  }
  if (iz == 10) return 0;
  
  // Find number of pixels in Earth views
  pcdim = 0;
  psdim = 0;
  linen = line0;
  for (size_t i=iz; i<9; i++) {
    // Check for Earth, solar, lunar or diagnostic mode
    if ( dtype[i] == 1 || dtype[i] == 3 || dtype[i] == 4 || dtype[i] == 6 ||
         dtype[i] == 7 || dtype[i] == 9) {
      uint16_t np = lines[i] / iagg[i];
      for (size_t j=0; j<np; j++) {
        clines[pcdim+j] = linen + j*iagg[i] + 0.5*(iagg[i]-1);
      }
      pcdim += np;
      uint16_t ns = lines[i] / 8;
      for (size_t j=0; j<ns; j++) {
        slines[psdim+j] = linen + j*8 + 3.5;
      }
      psdim += ns;
      iret = 0;
    }
    linen += lines[i];
  }
 
  // Calculate times
  for (size_t i=0; i<(size_t) pcdim; i++) deltc[i] = secpline * clines[i];
  ev_toff = 0.5 * (deltc[0] + deltc[pcdim-1]);
  for (size_t i=0; i<(size_t) pcdim; i++) deltc[i] -= ev_toff;
    
  for (size_t i=0; i<(size_t) psdim; i++)
    delts[i] = secpline * slines[i] - ev_toff;
 
  return 0;
}
 
 
int get_agg_mat( size_t *ia, int16_t iagg, int16_t jagg[16], uint16_t nib,
                 uint32_t& nbb, uint32_t ntb[16], float **amat, float **gmat) {
 
  // Populate gain aggregation matrix
 
  for (size_t i=0; i<512; i++) {
    gmat[i] = new float[nib];
    for (size_t j=0; j<nib; j++) gmat[i][j] = 0.0;
  }
  
  int16_t ii = 0;
  int16_t itt[16][2];
  for (size_t i=0; i<16; i++) {
    if ( jagg[i] > 0) {
      int16_t iaf = 4;
      if ( iagg*jagg[i] < iaf) iaf = iagg*jagg[i];
      for (size_t k=0; k<ntb[i]; k++) {
        size_t ic = 32*i;
        size_t kj = k*jagg[i];
        for (size_t j=0; j<(size_t) jagg[i]; j++)
          gmat[ic+kj+j][ii+k] = (1.0/jagg[i]) * (4/iaf);
      }
      itt[i][0] = ii;
      ii += ntb[i];
      itt[i][1] = ii - 1;
    }
  }
 
  // Compute number of bands for 8x aggregation with overlapping bands
  nbb = (ntb[0]*3) / 4 + 1;
  //amat(nbb,nib)
  for (size_t i=1; i<16; i++) {
    if (jagg[i] >= jagg[i-1])
      nbb += ntb[i];
    else
      nbb += (ntb[i]*3) / 4 + ntb[i-1]/4;
  }
  
  for (size_t i=0; i<nib; i++) {
    amat[i] = new float[nbb];
    for (size_t j=0; j<nbb; j++) amat[i][j] = 0;
  }
 
  // First tap
  int16_t ib;
  for (size_t k=0; k<=(ntb[ia[0]]*3)/4; k++) {
    amat[k+ntb[ia[0]]/4-1][k] = jagg[ia[0]]/8.0;
  }
  ib = (ntb[ia[0]]*3)/4 + 1;
 
  // Remaining taps
  uint16_t nr;
  for (size_t i=ia[1]; i<16; i++) {
    if (ntb[i] > 0) {
      if (ntb[i] >= ntb[i-1]) {
        // Transition resolution determined by preceding tap
        nr = ntb[i-1]/4 - 1;
        // Remaining bands using preceding tap
        if (nr > 0) {
          for (size_t k=0; k<nr; k++) {
        uint16_t k1 = nr - k - 1;
            uint16_t k2 = ((k+1)*ntb[i]) / ntb[i-1] - 1;
            for (size_t j=0; j<k1; j++)
              amat[itt[i-1][1]-k1+j][ib+k] = jagg[i-1] / 8.0;
            for (size_t j=0; j<k2; j++)
              amat[itt[i][0]+j][ib+k] = jagg[i] / 8.0;
          }
          ib += nr;
        }
      } else {
        // Transition resolution determined using current tap
        nr = ntb[i]/4 - 1;
        // Remaining bands using previous tap
        if (nr > 0) {
          for (size_t k=0; k<nr; k++) {
            uint16_t k1 = ((nr-k)*ntb[i-1]) / ntb[i] - 1;
        uint16_t k2 = k;
            for (size_t j=0; j<k1; j++)
              amat[itt[i-1][1]-k1+j][ib+k] = jagg[i-1] / 8.0;
            for (size_t j=0; j<k2; j++)
              amat[itt[i][0]+j][ib+k] = jagg[i] / 8.0;
          }            
      ib += nr;
        }
      }
      // Remaining bands using this tap
      for (size_t k=0; k<=(ntb[i]*3)/4; k++) {
        for (size_t j=0; j<ntb[i]/4; j++)
        amat[itt[i][0]+k+j][ib+k] = jagg[i] / 8.0;
      }
      ib += (ntb[i]*3) / 4 + 1;
    }
  }
                            
  return 0;
}
 
int read_oci_cal_lut( NcFile *calLUTfile, string tag, NcGroup gidLUT,
                      uint32_t& banddim, cal_lut_struct& cal_lut) {
 
  size_t dims[3];
  NcDim ncDIM;
  uint32_t timedim, tempdim, tcdim, rvsdim, nldim, msdim=2;
  string bandname;
  bandname.assign( tag);
  bandname.append( "_bands");
  
  ncDIM = calLUTfile->getDim(bandname.c_str());
  banddim = ncDIM.getSize();
  
  ncDIM = calLUTfile->getDim("number_of_times");
  timedim = ncDIM.getSize();
  ncDIM = calLUTfile->getDim("number_of_temperatures");
  tempdim = ncDIM.getSize();
  ncDIM = calLUTfile->getDim("number_of_T_coefficients");
  tcdim = ncDIM.getSize();
  ncDIM = calLUTfile->getDim("number_of_RVS_coefficients");
  rvsdim = ncDIM.getSize();
  ncDIM = calLUTfile->getDim("number_of_nonlinearity_coefficients");
  nldim = ncDIM.getSize();
 
  float **K1 = new float *[banddim];
  K1[0] = new float[msdim*banddim];
  for (size_t i=1; i<banddim; i++) K1[i] = K1[i-1] + msdim;
  cal_lut.K1 = K1;
 
  dims[0] = banddim; dims[1] = msdim; dims[2] = timedim;
  float ***K2 = make3dT<float>(dims);
  cal_lut.K2 = K2;
  dims[0] = banddim; dims[1] = tempdim; dims[2] = tcdim;
  float ***K3_coef = make3dT<float>(dims);
  cal_lut.K3_coef = K3_coef;
  dims[0] = banddim; dims[1] = msdim; dims[2] = rvsdim;
  float ***K4_coef = make3dT<float>(dims);
  cal_lut.K4_coef = K4_coef;
 
  double **K5 = new double *[banddim];
  K5[0] = new double[nldim*banddim];
  for (size_t i=1; i<banddim; i++) K5[i] = K5[i-1] + nldim;
  cal_lut.K5 = K5;
  
  cal_lut.ldims[0] = timedim; cal_lut.ldims[1] = tempdim;
  cal_lut.ldims[2] = tcdim; cal_lut.ldims[3] = rvsdim;
  cal_lut.ldims[4] = nldim; cal_lut.ldims[5] = msdim;
 
  NcVar var;
  
  var = gidLUT.getVar( "K1");
  var.getVar( &cal_lut.K1[0][0]);
  var = gidLUT.getVar( "K2");
  var.getVar( &cal_lut.K2[0][0][0]);
  var = gidLUT.getVar( "K3_coef");
  var.getVar( &cal_lut.K3_coef[0][0][0]);
  var = gidLUT.getVar( "K4_coef");
  var.getVar( &cal_lut.K4_coef[0][0][0]);
  var = gidLUT.getVar( "K5");
  var.getVar( &cal_lut.K5[0][0]);
  
  return 0;
}
 
// float  K1K2[nib][msdim]: absolute and temporal gain
// float  K3_coef[nib][tempdim][tcdim]: temperature correction
// float  K4_coef[nib][msdim][rvsdim]: RVS correction
// double K5[nib][nldim]: linearity correction
 
int make_oci_gains( uint32_t nib, uint32_t banddim, uint16_t iyr, uint16_t idom,
                    double stime, double K2t[NTIMES], cal_lut_struct& cal_lut,
                    float **gmat, gains_struct& gains) {
 
  for (size_t i=0; i<6; i++) gains.ldims[i] = cal_lut.ldims[i];
 
  uint16_t timedim = gains.ldims[0];
  uint16_t tempdim = gains.ldims[1];
  uint16_t tcdim = gains.ldims[2];
  uint16_t rvsdim = gains.ldims[3];
  uint16_t nldim = gains.ldims[4];
  uint16_t msdim = gains.ldims[5];
 
  size_t dims[3];
  
  float **K1K2 = new float *[nib];
  K1K2[0] = new float[nib*msdim];
  for (size_t i=1; i<nib; i++) K1K2[i] = K1K2[i-1] + msdim;
  gains.K1K2 = K1K2;
 
  dims[0] = nib; dims[1] = tempdim; dims[2] = tcdim;
  float ***K3_coef = make3dT<float>(dims);
  gains.K3_coef = K3_coef;
  dims[0] = nib; dims[1] = msdim; dims[2] = rvsdim;
  float ***K4_coef = make3dT<float>(dims);
  gains.K4_coef = K4_coef;
 
  double **K5 = new double *[nib];
  K5[0] = new double[nib*nldim];
  for (size_t i=1; i<nib; i++) K5[i] = K5[i-1] + nldim;
  gains.K5 = K5;
  
  // Mirror-side dependent gains
  double *K2 = new double[banddim];
  for (size_t ms=0; ms<msdim; ms++) {
 
    // Get temporal gain and combine with absolute gain
    double d2 = jday(iyr, 1, idom) - 2451545 + stime/864.0;
 
    size_t kd=0;
    for (size_t j=NTIMES-1; j>=0; j--) {
      if ( d2 > K2t[j]) {
        kd = j;
        break;
      }
    }
    if ( kd < (size_t) (timedim-1)) {
      double ff = (d2 - K2t[kd]) / (K2t[kd+1] - K2t[kd]);
      for (size_t j=0; j<banddim; j++)
        K2[j] = cal_lut.K2[j][ms][kd]*(1.0-ff) + cal_lut.K2[j][ms][kd+1]*ff;
    } else {
      for (size_t j=0; j<banddim; j++) K2[j] = cal_lut.K2[j][ms][kd];
    }
    for (size_t j=0; j<banddim; j++) K2[j] *= cal_lut.K1[j][ms];
    for (size_t i=0; i<nib; i++) {
      gains.K1K2[i][ms] = 0;
      for (size_t j=0; j<banddim; j++)
        gains.K1K2[i][ms] += gmat[j][i] * cal_lut.K1[j][ms];
    }
 
    // Generate RVS coefficents
    for (size_t i=0; i<nib; i++) {
      for (size_t k=0; k<rvsdim; k++) {
        gains.K4_coef[i][ms][k] = 0;
        for (size_t j=0; j<banddim; j++)
          gains.K4_coef[i][ms][k] += gmat[j][i] * cal_lut.K4_coef[j][ms][k];
      }
    }
  }
  delete [] K2;
  
  // Generate temperature coefficients
  for (size_t i=0; i<nib; i++) {
    for (size_t k=0; k<tcdim; k++) {
      for (size_t l=0; l<tempdim; l++) {
        gains.K3_coef[i][l][k] = 0;
        for (size_t j=0; j<banddim; j++)
          gains.K3_coef[i][l][k] += gmat[j][i] * cal_lut.K3_coef[j][l][k];
      }
    }
  }
 
  // Generate linearity coefficients
  for (size_t i=0; i<nib; i++) {
    for (size_t k=0; k<nldim; k++) {
      gains.K5[i][k] = 0;
      for (size_t j=0; j<banddim; j++)
        gains.K5[i][k] += gmat[j][i] * cal_lut.K5[j][k];
    }
  }
  /*
  for (size_t i=0; i<nib; i++) {
    for (size_t j=0; j<banddim; j++)
      cout << "2: " << i << " " << j << " " << gmat[j][i] << endl;
  }
  */
  return 0;
}
 
 
int get_oci_vecs( uint32_t nscan, uint16_t pdim, geo_struct& geoLUT,
                  double ev_toff, int32_t spin, float *slines, double *delt,
                  double revpsec, int32_t ppr_off, int32_t *mspin,
                  int32_t *ot_10us, uint8_t *enc_count, float **hamenc,
                  float **rtaenc, float **pview, double *theta, int16_t& iret) {
 
  // This program generates the OCI Earth view vectors for one spin.  
  // It uses MCE telemetry and encoder data.  Further refinements will be made
  // as the instrument optics model and test results become available.  
  // Reference: "Use of OCI Telemetry to Determine Pixel Line-of-Sight",
  // F. Patt, 2020-05-18
 
  int32_t max_enc_cts = 131072; // 2^17
  double dtenc = 0.001;
  constexpr double pi = 3.14159265358979323846;
  
  double rad2asec = (180/pi) * 3600;
 
  // Compute scan angle corresponding to PPR
  float pprang = 2 * pi * (ppr_off - geoLUT.rta_nadir) / max_enc_cts;
  if (pprang > pi) pprang -= 2*pi;
  
  // Compute ideal scan angles for science pixels
  double *toff = new double[pdim];
  for (size_t i=0; i<pdim; i++) toff[i] = delt[i] + ev_toff;
 
  for (size_t i=0; i<pdim; i++)
    theta[i] = pprang + 2 * pi * revpsec * toff[i];
  // Interpolate encoder data to pixel times and add to scan angles
  // RTA only for now, include HAM when we know how.
 
  double *thetacor = new double[pdim];
 
  int isp = -1;
  for (size_t i=0; i<nscan; i++) {
    if ( mspin[i] == spin) {
      isp = (int) i;
      break;
    }
  }
  if ( isp == -1) {
    cout << "No MCE encoder data for spin: " << spin << endl;
    iret = 1;
  }
  size_t ip = 0, ke;
  double tenc_ke;
  while ( ip < pdim) {
    // Encoder sample times at 1 KHz
    for (size_t j=0; j<enc_count[isp]; j++) {
      double tenc = j*dtenc;
      if ( tenc < toff[ip] && (tenc+dtenc) > toff[ip]) {
        ke = j;
        tenc_ke = tenc;
        break;
      }
    }
    size_t njp=0;
    for (size_t i=0; i<pdim; i++) {
      if ( toff[i] >= tenc_ke && toff[i] < tenc_ke+dtenc) {
        double ft = (toff[i] - tenc_ke) / dtenc;
        thetacor[i] = (1-ft)*rtaenc[isp][ke] + ft*rtaenc[isp][ke+1];
        njp++;
      }
    }
    ip += njp;
  }
 
  // Simple, planar view vector model to start
  // Update when optical model is available.
  for (size_t i=0; i<pdim; i++) {
    theta[i] = theta[i] + thetacor[i] / rad2asec;
    pview[i][1] = sin(theta[i]);
    pview[i][2] = cos(theta[i]);
  }
 
  delete [] thetacor;
  delete [] toff;
 
  return 0;
}
 
 
template <typename T>
int get_oci_dark( size_t iscn, uint32_t nscan, uint8_t *hside, uint16_t ndsc,
                  uint16_t nskp, int16_t iags, int16_t iagd, uint32_t ntaps,
                  int16_t *jagg, uint32_t dfill, int16_t ndc, T ***dark,
                  uint32_t nib, float *dc, int16_t& iret) {
 
  // Program to generate dark corrections for OCI data by averaging the
  // dark collect data and correcting for bit shift/truncation if necessary
 
  // Determine number of bands per tap for hyperspectral data
  int16_t *nbndt = new int16_t[ntaps];
 
  if (ntaps == 16) {
    // hyperspectral bands
    for (size_t i=0; i<ntaps; i++)
      if ( jagg[i] > 0) nbndt[i] = 32 / jagg[i]; else nbndt[i] = 0;
  } else {
    for (size_t i=0; i<ntaps; i++) nbndt[i] = 9;
  }
  int16_t nbnd = 0;
  for (size_t i=0; i<ntaps; i++) nbnd += nbndt[i];
 
  // Select data for HAM side and determine scan indices
  int32_t *kh = new int32_t[nscan];
  int32_t nkh=0;
  for (size_t i=0; i<nscan; i++) {
    if ( hside[i] == hside[iscn]) {
      kh[i] = (int32_t) i;
      nkh++;
    } else {
      kh[i] = -1;
    }
  }
 
  int32_t js=0;
  for (size_t i=0; i<nscan; i++) {
    if ( kh[i] == (int32_t) iscn) {
      js = (int32_t) i;
      break;
    }
  }
 
  // Check for valid dark collect data within specified range
  uint16_t ndscl = ndsc;
  bool valid_dark_found = false;
  int32_t is1=js, is2=js;
 
  while (!valid_dark_found && ndscl <= nkh) {
    if ( ndsc > 1) {
      is1 = js - ndsc/2;
      is2 = js + ndsc/2;
      // Check for start or end of granule
      if (is1 < 0) {
        is1 = 0;
        is2 = ndsc - 1;
      }
      if (is2 >= nkh) {
        is1 = nkh - ndsc;
        is2 = nkh - 1;
      }
    }
 
    // If no valid dark data, expand scan range
    for (size_t i=is1; i<=(size_t) is2; i++) {
      for (size_t j=nskp; j<(size_t) ndc; j++) {
        if ( dark[kh[i]][0][j] != dfill) {
          valid_dark_found = true;
          break;
        }
      }
    }
    if ( !valid_dark_found) ndscl += 2;
  }
 
  if ( !valid_dark_found) {
    iret =-1;
    return 0;
  }
  
  // Loop through taps and compute dark correction
  int16_t ibnd=0;
  for (size_t i=0; i<ntaps; i++) {
    if ( jagg[i] > 0) {
      float ddiv = 1.0;
      float doff = 0.0;
      if ( iags*jagg[i] > 4) {
        ddiv = iagd * jagg[i] / 4.0;
        doff = (ddiv-1) / (2*ddiv);
      }
 
      for (size_t j=0; j<(size_t) nbndt[i]; j++) {
        float sum = 0.0;
        int nv = 0;
        for (size_t k=kh[is1]; k<=(size_t) kh[is2]; k++) {
          for (size_t l=nskp; l<(size_t) ndc; l++) {
            if (dark[k][ibnd+j][l] != dfill) {
              sum += dark[k][ibnd+j][l];
              nv++;
            }
          }
        }
        //        dc[ibnd+j] = ((sum/(ndc-nskp))/ddiv) - doff;
        dc[ibnd+j] = ((sum/(nv))/ddiv) - doff;
      } // j loop
    } // if ( 
    ibnd += nbndt[i];
  } // i loop
 
  delete [] kh;
  delete [] nbndt;
 
  iret = 0;
  if (ndscl > ndsc) iret = 1;
  
  return 0;
}
 
 
int get_oci_temp_corr( uint32_t nib, gains_struct gains, float K3T[NTEMPS],
                       float *caltemps, uint32_t nscan, float *k3) {
 
  uint16_t tempdim = gains.ldims[1];
  uint16_t tcdim = gains.ldims[2];
 
  for (size_t i=0; i<nib; i++) k3[i] = 1.0;
 
  for (size_t i=0; i<tempdim; i++) {
    float td = caltemps[i] - K3T[i];
      for (size_t j=0; j<tcdim; j++) {
        for (size_t k=0; k<nib; k++) {
          k3[k] -= gains.K3_coef[k][i][j] * powf(td, j+1);
        }
      }
  }
 
  return 0;
}
 
int get_oci_rvs_corr( uint32_t nib, uint16_t pdim, uint8_t hside,
                      gains_struct gains, double *theta, float **k4) {
 
  // Program to compute RVS correction from coefficients and scan angle
 
  uint16_t rvsdim = gains.ldims[3];
 
  for (size_t i=0; i<nib; i++)
    for (size_t j=0; j<pdim; j++)
      k4[i][j] = 1.0;
 
  for (size_t i=0; i<rvsdim; i++)
    for (size_t j=0; j<nib; j++)
      for (size_t k=0; k<pdim; k++)
        k4[j][k] += gains.K4_coef[j][hside][i] * powf(theta[k], j+1);
  
  return 0;
}
 
int get_oci_lin_corr( uint32_t nib, uint16_t pdim, gains_struct gains,
                      float K3T[NTEMPS], float *caltemps, float **dn,
                      float **k5) {
 
  // Program to compute OCI linearity correction from coefficients, dn and
  // temperatures 3rd-order polynomial of dn with 2nd-order temperature effects
  // for three temperatures.
  // This is currently a SWAG at the functional form and will be revised when
  // better known
 
  // Index for which temperatures to use for linearity; revisit this when known
  int16_t ibtmp[3] = {1,2,3};
  float td[3];
  for (size_t i=0; i<3; i++) td[i] = caltemps[ibtmp[i]] - K3T[ibtmp[i]];
 
  for (size_t i=0; i<pdim; i++) {
    // Zeroth-order correction
    for (size_t j=0; j<nib; j++) k5[j][i] = gains.K5[j][0];
    // First and second-order corrections are quadratic functions of
    // temperatures
    for (size_t k=1; k<=2; k++) {
      for (size_t l=0; l<nib; l++) {
        k5[l][i] += (gains.K5[l][k] + gains.K5[l][k+3]*powf(td[0], k) +
                     gains.K5[l][k+5]*powf(td[1], k) +
                     gains.K5[l][k+7]*powf(td[2], k))*powf(dn[l][i], k);
      }
      // Third-order correction
      for (size_t l=0; l<nib; l++)
        k5[l][i] += gains.K5[l][3]*powf(dn[l][i], 3);
    }
  }
 
  return 0;
}
 
 
/*------------------------------------------------------------------ */
/* Create an Generic NETCDF4 level1B file                            */
/* ----------------------------------------------------------------- */
int l1bFile::createl1b( char* l1b_filename, uint16_t nscan_good,
                        uint16_t pcdim, uint16_t bbb, uint16_t rbb,
                        uint16_t psdim, uint16_t swb) {
 
  // nscan_good: Number of scans in L1A file
  // pcdim:      Number of CCD band pixels 
  // bbb:        Number of blue bands]
  // rbb:        Number of red bands
  // psdim:      Number of SWIR band pixels
  // swb:        Number of SWIR bands
  
  try {
    l1bfile = new NcFile( l1b_filename, NcFile::replace);
  }
  catch ( NcException& e) {
    e.what();
    cerr << "Failure creating OCI L1B file: " << l1b_filename << endl;
    exit(1);
  }
 
  fileName.assign( l1b_filename);
   
  ifstream oci_l1b_data_structure;
  string line;
  string dataStructureFile;
  dataStructureFile.assign("$OCDATAROOT/oci/OCI_Level-1B_Data_Structure.cdl");
  expandEnvVar( &dataStructureFile);
   
  oci_l1b_data_structure.open( dataStructureFile.c_str(), ifstream::in);
  if ( oci_l1b_data_structure.fail() == true) {
    cout << "\"" << dataStructureFile.c_str() << "\" not found" << endl;
    exit(1);
  }
 
  // Find "dimensions" section of CDL file
  while(1) {
    getline( oci_l1b_data_structure, line);
    size_t pos = line.find("dimensions:");
    if ( pos == 0) break;
  }
 
  // Define dimensions from "dimensions" section of CDL file
  ndims = 0;
 
  while(1) {
     getline( oci_l1b_data_structure, line);
     if (line.substr(0,2) == "//") continue;
     
     size_t pos = line.find(" = ");
     size_t semi = line.find(" ;");
     if ( pos == string::npos) break;
 
     uint32_t dimSize;
     istringstream iss;
     //     istringstream iss(line.substr(pos+2, string::npos));
     string dimString = line.substr(pos+2, semi-(pos+2));
     if (dimString.find("UNLIMITED") == string::npos) {
       iss.str(dimString);
       iss >> dimSize;
     } else {
       dimSize = NC_UNLIMITED;
     }
     
     iss.clear(); 
     iss.str( line);
     iss >> skipws >> line;
 
     if (line.compare("ccd_pixels") == 0) {
       dimSize = pcdim;
     }
 
     if (line.compare("SWIR_pixels") == 0) {
       dimSize = psdim;
     }
 
     if (line.compare("blue_bands") == 0) {
       dimSize = bbb;
     }
 
     if (line.compare("red_bands") == 0) {
       dimSize = rbb;
     }
 
     if (line.compare("SWIR_bands") == 0) {
       dimSize = swb;
     }
 
     try {
       ncDims[ndims++] = l1bfile->addDim( line, dimSize);
     }
     catch ( NcException& e) {
       e.what();
       cerr << "Failure creating dimension: " << line.c_str() << endl;
       exit(1);
     }
 
     cout << "Dimension Name: " << line.c_str() << " Dimension Size: "
          << dimSize << endl;
 
   } // while loop
 
   // Read global attributes (string attributes only) 
   while(1) {
     getline( oci_l1b_data_structure, line);
     size_t pos = line.find("// global attributes");
     if ( pos == 0) break;
   }
 
   while(1) {
     getline( oci_l1b_data_structure, line);
     size_t pos = line.find(" = ");
     if ( pos == string::npos) break;
 
     string attValue;
 
     // Remove leading and trailing quotes
     attValue.assign(line.substr(pos+4));
     size_t posQuote = attValue.find("\"");
     attValue.assign(attValue.substr(0, posQuote));
 
     istringstream iss(line.substr(pos+2));
     iss.clear(); 
     iss.str( line);
     iss >> skipws >> line;
 
     // Skip commented out attributes
     if (line.substr(0,2) == "//") continue;
 
     string attName;
     attName.assign(line.substr(1).c_str());
 
     // Skip non-string attributes
     if (attName.compare("orbit_number") == 0) continue;
     if (attName.compare("history") == 0) continue;
     if (attName.compare("format_version") == 0) continue;
     if (attName.compare("instrument_number") == 0) continue;
     if (attName.compare("pixel_offset") == 0) continue;
     if (attName.compare("number_of_filled_scans") == 0) continue;
 
     // cout << attName.c_str() << " " << attValue.c_str() << endl;
     try {
       l1bfile->putAtt(attName, attValue);
     }
     catch ( NcException& e) {
       e.what();
       cerr << "Failure creating attribute: " + attName << endl;
       exit(1);
     }
     
   } // while(1)
 
 
   ngrps = 0;
   // Loop through groups
   while(1) {
     getline( oci_l1b_data_structure, line);
 
     // Check if end of CDL file
     // If so then close CDL file and return
     if (line.substr(0,1).compare("}") == 0) {
       oci_l1b_data_structure.close();
       return 0;
     }
 
     // Check for beginning of new group
     size_t pos = line.find("group:");
 
     // If found then create new group and variables
     if ( pos == 0) {
 
       // Parse group name
       istringstream iss(line.substr(6, string::npos));
       iss >> skipws >> line;
 
       ncGrps[ngrps++] = l1bfile->addGroup(line);
 
       int numDims=0;
       string sname;
       string lname;
       string standard_name;
       string units;
       string flag_values;
       string flag_meanings;
       double valid_min=0.0;
       double valid_max=0.0;
       double fill_value=0.0;
 
       vector<NcDim> varVec;
 
       int ntype=0;
       NcType ncType;
       
       // Loop through datasets in group
       getline( oci_l1b_data_structure, line); // skip "variables:"
       while(1) {
         getline( oci_l1b_data_structure, line);
 
         if (line.substr(0,2) == "//") continue;
         if (line.length() == 0) continue;
         if (line.substr(0,1).compare("\r") == 0) continue;
         if (line.substr(0,1).compare("\n") == 0) continue;
 
         size_t pos = line.find(":");
 
         // No ":" found, new dataset or empty line or end-of-group
         if ( pos == string::npos) {
           
           if ( numDims > 0) {
             // Create previous dataset
 
             createField( ncGrps[ngrps-1], sname.c_str(), lname.c_str(),
                          standard_name.c_str(), units.c_str(),
                          (void *) &fill_value, 
                          flag_values.c_str(), flag_meanings.c_str(),
                          valid_min, valid_max, ntype, varVec);
 
             flag_values.assign("");
             flag_meanings.assign("");
             units.assign("");
             varVec.clear();
           }
 
           valid_min=0.0;
           valid_max=0.0;
           fill_value=0.0;
 
           if (line.substr(0,10).compare("} // group") == 0) break;
 
           // Parse variable type
           string varType;
           istringstream iss(line);
           iss >> skipws >> varType;
 
           // Get corresponding NC variable type
           if ( varType.compare("char") == 0) ntype = NC_CHAR;
           else if ( varType.compare("byte") == 0) ntype = NC_BYTE;
           else if ( varType.compare("short") == 0) ntype = NC_SHORT;
           else if ( varType.compare("int") == 0) ntype = NC_INT;
           else if ( varType.compare("long") == 0) ntype = NC_INT;
           else if ( varType.compare("float") == 0) ntype = NC_FLOAT;
           else if ( varType.compare("real") == 0) ntype = NC_FLOAT;
           else if ( varType.compare("double") == 0) ntype = NC_DOUBLE;
           else if ( varType.compare("ubyte") == 0) ntype = NC_UBYTE;
           else if ( varType.compare("ushort") == 0) ntype = NC_USHORT;
           else if ( varType.compare("uint") == 0) ntype = NC_UINT;
           else if ( varType.compare("ulong") == 0) ntype = NC_UINT;
           else if ( varType.compare("int64") == 0) ntype = NC_INT64;
           else if ( varType.compare("uint64") == 0) ntype = NC_UINT64;
 
           // Parse short name (sname)
           pos = line.find("(");
           size_t posSname = line.substr(0, pos).rfind(" ");
           sname.assign(line.substr(posSname+1, pos-posSname-1));
           cout << "sname: " << sname.c_str() << endl;
 
           // Parse variable dimension info
           this->parseDims( line.substr(pos+1, string::npos), varVec);
           numDims = varVec.size();
 
         } else {
           // Parse variable attributes
           size_t posEql = line.find("=");
           size_t pos1qte = line.find("\"");
           size_t pos2qte = line.substr(pos1qte+1, string::npos).find("\"");
       // cout << line.substr(pos+1, posEql-pos-2).c_str() << endl;
 
           string attrName = line.substr(pos+1, posEql-pos-2);
 
           // Get long_name
           if ( attrName.compare("long_name") == 0) {
             lname.assign(line.substr(pos1qte+1, pos2qte));
             //             cout << "lname: " << lname.c_str() << endl;
           }
 
           // Get units
           else if ( attrName.compare("units") == 0) {
             units.assign(line.substr(pos1qte+1, pos2qte));
             //             cout << "units: " << units.c_str() << endl;
           }
 
           // Get _FillValue
           else if ( attrName.compare("_FillValue") == 0) {
             iss.clear(); 
             iss.str( line.substr(posEql+1, string::npos));
             iss >> fill_value;
           }
 
           // Get flag_values
           else if ( attrName.compare("flag_values") == 0) {
             flag_values.assign(line.substr(pos1qte+1, pos2qte));
           }
 
           // Get flag_meanings
           else if ( attrName.compare("flag_meanings") == 0) {
             flag_meanings.assign(line.substr(pos1qte+1, pos2qte));
           }
 
           // Get valid_min
           else if ( attrName.compare("valid_min") == 0) {
             iss.clear(); 
             iss.str( line.substr(posEql+1, string::npos));
             iss >> valid_min;
             //             cout << "valid_min: " << valid_min << endl;
           }
 
           // Get valid_max
           else if ( attrName.compare("valid_max") == 0) {
             iss.clear(); 
             iss.str( line.substr(posEql+1, string::npos));
             iss >> valid_max;
             //             cout << "valid_max: " << valid_max << endl;
           }
 
         } // if ( pos == string::npos)
       } // datasets in group loop
     } // New Group loop
   } // Main Group loop
 
   return 0;
}
 
 
int createField( NcGroup &ncGrp, const char *sname, const char *lname, 
                 const char *standard_name, const char *units,
                 void *fill_value, const char *flag_values,
                 const char *flag_meanings,
                 double low, double high, int nt, vector<NcDim>& varVec) {
 
  size_t dimlength;
 
 
  /* Create the NCDF dataset */
  NcVar ncVar;
  try {
    ncVar = ncGrp.addVar(sname, nt, varVec);   
  }
  catch ( NcException& e) {
    cout << e.what() << endl;
    cerr << "Failure creating variable: " << sname << endl;
    exit(1);
  }
 
  // Set fill value
  double fill_value_dbl;
  memcpy( &fill_value_dbl, fill_value, sizeof(double));
 
  int8_t i8;
  uint8_t ui8;
  int16_t i16;
  uint16_t ui16;
  int32_t i32;
  uint32_t ui32;
  float f32;
 
  if ( low != fill_value_dbl) {
    if ( nt == NC_BYTE) {
      i8 = fill_value_dbl;
      ncVar.setFill(true, (void *) &i8);
    } else if ( nt == NC_UBYTE) {
      ui8 = fill_value_dbl;
      ncVar.setFill(true, (void *) &ui8);
    } else if ( nt == NC_SHORT) {
      i16 = fill_value_dbl;
      ncVar.setFill(true, (void *) &i16);
    } else if ( nt == NC_USHORT) {
      ui16 = fill_value_dbl;
      ncVar.setFill(true, (void *) &ui16);
    } else if ( nt == NC_INT) {
      i32 = fill_value_dbl;
      ncVar.setFill(true, (void *) &i32);
    } else if ( nt == NC_UINT) {
      ui32 = fill_value_dbl;
      ncVar.setFill(true, (void *) &ui32);
    } else if ( nt == NC_FLOAT) {
      f32 = fill_value_dbl;
      ncVar.setFill(true, (void *) &f32);
    } else {
      ncVar.setFill(true, (void *) &fill_value_dbl);
    }
  }
 
  /* vary chunck size based on dimensions */ 
  int do_deflate = 0;
  vector<size_t> chunkVec;
  if ( varVec.size() == 3 && (strncmp(sname, "EV_", 3) == 0)) {
    dimlength = varVec[2].getSize();
 
    chunkVec.push_back(1);
    chunkVec.push_back(16);
    chunkVec.push_back(dimlength/10);
 
    do_deflate = 1;
  }
 
  /* Set compression */
  if ( do_deflate) {
    /* First set chunking */
    try {
      ncVar.setChunking(ncVar.nc_CHUNKED, chunkVec);
    }
    catch ( NcException& e) {
      e.what();
      cerr << "Failure setting chunking: " << sname << endl;
      exit(1);
    }
 
    try {
      ncVar.setCompression(true, true, 5);
    }
    catch ( NcException& e) {
      e.what();
      cerr << "Failure setting compression: " << sname << endl;
      exit(1);
    }
  }
 
 
  /* Add a "long_name" attribute */
  try {
    ncVar.putAtt("long_name", lname);
  }
  catch ( NcException& e) {
    e.what();
    cerr << "Failure creating 'long_name' attribute: " << lname << endl;
    exit(1);
  }
 
  if ( strcmp( flag_values, "") != 0) {
 
    size_t curPos=0;
 
    string fv;
    fv.assign( flag_values);
    size_t pos = fv.find("=", curPos);
    fv = fv.substr(pos+1);
 
    size_t semicln = fv.find(";");
    pos = 0;
 
    int8_t vec[1024];
    int n = 0;
    while(pos != semicln) {
      pos = fv.find(",", curPos);
      if ( pos == string::npos) 
        pos = semicln;
 
      string flag_value;
      istringstream iss(fv.substr(curPos, pos-curPos));
      iss >> skipws >> flag_value;
      vec[n++] = atoi( flag_value.c_str());
      curPos = pos + 1;
    }
 
    try {
      ncVar.putAtt("flag_values", NC_BYTE, n, vec);
    }
    catch ( NcException& e) {
      e.what();
      cerr << "Failure creating 'flag_values' attribute: " << lname << endl;
      exit(1);
    }
  }
 
  /* Add a "flag_meanings" attribute if specified*/
  if ( strcmp( flag_meanings, "") != 0) {
 
    try {
      ncVar.putAtt("flag_meanings", flag_meanings);
    }
    catch ( NcException& e) {
      e.what();
      cerr << "Failure creating 'flag_meanings' attribute: "
           << flag_meanings << endl;
      exit(1);
    }
  }
 
  /* Add "valid_min/max" attributes if specified */
  if (low < high) {
    switch(nt) {              /* Use the appropriate number type */
    case NC_BYTE:
      {
    uint8_t vr[2];
    vr[0] = (uint8_t)low;
    vr[1] = (uint8_t)high;
 
        try {
          ncVar.putAtt("valid_min", NC_BYTE, 1, &vr[0]);
        }
        catch ( NcException& e) {
          e.what();
          cerr << "Failure creating 'valid_min' attribute: " << vr[0] << endl;
          exit(1);
        }
 
        try {
          ncVar.putAtt("valid_max", NC_BYTE, 1, &vr[1]);
        }
        catch ( NcException& e) {
          e.what();
          cerr << "Failure creating 'valid_max' attribute: " << vr[1] << endl;
          exit(1);
        }
      }
      break;
    case NC_UBYTE:
      {
    uint8_t vr[2];
    vr[0] = (uint8_t)low;
    vr[1] = (uint8_t)high;
 
        try {
          ncVar.putAtt("valid_min", NC_UBYTE, 1, &vr[0]);
        }
        catch ( NcException& e) {
          e.what();
          cerr << "Failure creating 'valid_min' attribute: " << vr[0] << endl;
          exit(1);
        }
 
        try {
          ncVar.putAtt("valid_max", NC_UBYTE, 1, &vr[1]);
        }
        catch ( NcException& e) {
          e.what();
          cerr << "Failure creating 'valid_max' attribute: " << vr[1] << endl;
          exit(1);
        }
      }
      break;
    case NC_SHORT:
      {
    int16_t vr[2];
    vr[0] = (int16_t)low;
    vr[1] = (int16_t)high;
 
        try {
          ncVar.putAtt("valid_min", NC_SHORT, 1, &vr[0]);
        }
        catch ( NcException& e) {
          e.what();
          cerr << "Failure creating 'valid_min' attribute: " << vr[0] << endl;
          exit(1);
        }
 
        try {
          ncVar.putAtt("valid_max", NC_SHORT, 1, &vr[1]);
        }
        catch ( NcException& e) {
          e.what();
          cerr << "Failure creating 'valid_max' attribute: " << vr[1] << endl;
          exit(1);
        }
      }
      break;
    case NC_USHORT:
      {
    uint16_t vr[2];
    vr[0] = (uint16_t)low;
    vr[1] = (uint16_t)high;
 
        try {
          ncVar.putAtt("valid_min", NC_USHORT, 1, &vr[0]);
        }
        catch ( NcException& e) {
          e.what();
          cerr << "Failure creating 'valid_min' attribute: " << vr[0] << endl;
          exit(1);
        }
 
        try {
          ncVar.putAtt("valid_max", NC_USHORT, 1, &vr[1]);
        }
        catch ( NcException& e) {
          e.what();
          cerr << "Failure creating 'valid_max' attribute: " << vr[1] << endl;
          exit(1);
        }
              }
      break;
    case NC_INT:
      {
    int32_t vr[2];
    vr[0] = (int32_t)low;
    vr[1] = (int32_t)high;
 
        try {
          ncVar.putAtt("valid_min", NC_INT, 1, &vr[0]);
        }
        catch ( NcException& e) {
          e.what();
          cerr << "Failure creating 'valid_min' attribute: " << vr[0] << endl;
          exit(1);
        }
 
        try {
          ncVar.putAtt("valid_max", NC_INT, 1, &vr[1]);
        }
        catch ( NcException& e) {
          e.what();
          cerr << "Failure creating 'valid_max' attribute: " << vr[1] << endl;
          exit(1);
        }
 
      }
      break;
    case NC_UINT:
      {
    uint32_t vr[2];
    vr[0] = (uint32_t)low;
    vr[1] = (uint32_t)high;
 
        try {
          ncVar.putAtt("valid_min", NC_UINT, 1, &vr[0]);
        }
        catch ( NcException& e) {
          e.what();
          cerr << "Failure creating 'valid_min' attribute: " << vr[0] << endl;
          exit(1);
        }
 
        try {
          ncVar.putAtt("valid_max", NC_UINT, 1, &vr[1]);
        }
        catch ( NcException& e) {
          e.what();
          cerr << "Failure creating 'valid_max' attribute: " << vr[1] << endl;
          exit(1);
        }
        
      }
      break;
    case NC_FLOAT:
      {
    float vr[2];
    vr[0] = (float)low;
    vr[1] = (float)high;
 
        try {
          ncVar.putAtt("valid_min", NC_FLOAT, 1, &vr[0]);
        }
        catch ( NcException& e) {
          e.what();
          cerr << "Failure creating 'valid_min' attribute: " << vr[0] << endl;
          exit(1);
        }
 
        try {
          ncVar.putAtt("valid_max", NC_FLOAT, 1, &vr[1]);
        }
        catch ( NcException& e) {
          e.what();
          cerr << "Failure creating 'valid_max' attribute: " << vr[1] << endl;
          exit(1);
        }        
      }
      break;
    case NC_DOUBLE:
      {
    double vr[2];
    vr[0] = low;
    vr[1] = high;
 
        try {
          ncVar.putAtt("valid_min", NC_DOUBLE, 1, &vr[0]);
        }
        catch ( NcException& e) {
          e.what();
          cerr << "Failure creating 'valid_min' attribute: " << vr[0] << endl;
          exit(1);
        }
 
        try {
          ncVar.putAtt("valid_max", NC_DOUBLE, 1, &vr[1]);
        }
        catch ( NcException& e) {
          e.what();
          cerr << "Failure creating 'valid_max' attribute: " << vr[1] << endl;
          exit(1);
        }
      }
      break;
    default:
      fprintf(stderr,"-E- %s line %d: ",__FILE__,__LINE__);
      fprintf(stderr,"Got unsupported number type (%d) ",nt);
      fprintf(stderr,"while trying to create NCDF variable, \"%s\", ",sname);
      return(1);
    }
  }           
    
  /* Add a "units" attribute if one is specified */
  if(units != NULL && *units != 0) {
 
    try {
      ncVar.putAtt("units", units);
    }
    catch ( NcException& e) {
      e.what();
      cerr << "Failure creating 'units' attribute: " << units << endl;
      exit(1);
    }
  }
 
  /* Add a "standard_name" attribute if one is specified */
  if(standard_name != NULL && *standard_name != 0) {
    try {
      ncVar.putAtt("standard_name", units);
    }
    catch ( NcException& e) {
      e.what();
      cerr << "Failure creating 'standard_name' attribute: "
           << standard_name << endl;
      exit(1);
    }
  }
  
  return 0;
}
 
 
int l1bFile::parseDims( string dimString, vector<NcDim>& varDims) {
 
  size_t curPos=0;
  //  char dimName[NC_MAX_NAME+1];
  string dimName;
 
  while(1) {
    size_t pos = dimString.find(",", curPos);
    if ( pos == string::npos) 
      pos = dimString.find(")");
 
    string varDimName;
    istringstream iss(dimString.substr(curPos, pos-curPos));
    iss >> skipws >> varDimName;
 
    for (int i=0; i<ndims; i++) {
      
      try {      
        dimName = ncDims[i].getName();
      }
      catch ( NcException& e) {
        e.what();
        cerr << "Failure accessing dimension: " + dimName << endl;
        exit(1);
      }
      
      if ( varDimName.compare(dimName) == 0) {
        cout << "     " << dimName << " " << ncDims[i].getSize() << endl;
        varDims.push_back(ncDims[i]);
        break;
      }
    }
    if ( dimString.substr(pos, 1).compare(")") == 0) break;
 
    curPos = pos + 1;
  }
 
  return 0;
}
 
 
int l1bFile::write_granule_metadata( std::string tstart, std::string tend,
                                     std::string l1b_name) {
 
  // Date created
  char buf[32];
  strcpy( buf, unix2isodate(now(), 'G'));
  l1bfile->putAtt("date_created", buf);
 
  l1bfile->putAtt("time_coverage_start", tstart.c_str());
 
  l1bfile->putAtt("time_coverage_end", tend.c_str());
 
  // Write product file name
  l1bfile->putAtt("product_name", l1b_name.c_str());
 
  return 0;
}
 
 
int l1bFile::close() {
 
  try { 
    l1bfile->close();
  }
  catch ( NcException& e) {
    cout << e.what() << endl;
    cerr << "Failure closing: " + fileName << endl;
    exit(1);
  }
  
  return 0;
}
 
 
template <typename T>
T*** make3dT( size_t dims[3]) {
 
  T ***arr3d = new T **[dims[0]];
 
  arr3d[0] = new T*[dims[0]*dims[1]];
  arr3d[0][0] = new T[dims[0]*dims[1]*dims[2]];
 
  for (size_t i=1; i<dims[0]; i++) arr3d[i] = arr3d[i-1] + dims[1];
 
  for (size_t i=0; i<dims[0]; i++) {
    if ( i > 0) arr3d[i][0] = arr3d[i-1][0] + dims[1]*dims[2];
    for (size_t j=1; j<dims[1]; j++)
      arr3d[i][j] = arr3d[i][j-1] + dims[2];
  }
 
  return arr3d;
}