#include "table.h"
#include "SPio.h"
#include "SPutils.h"

//-------------------------------------------------------------------------------
// Class table

// default constructor

table::table()
  : m_Parameters(),
    m_Spectra(),
    m_ModelName(BlankStr),
    m_ModelUnits(BlankStr),
    m_NumIntParams(0),
    m_NumAddParams(0),
    m_isError(false),
    m_isRedshift(false),
    m_isAdditive(true),
    m_isEscale(false),
    m_Energies(),
    m_EnergyUnits("keV"),
    m_LowEnergyLimit(0.0),
    m_HighEnergyLimit(0.0),
    m_Filename(BlankStr)
{
}

// Copy constructor

table::table(const table& a)
{
  copy(a);
}

// Read constructor

table::table(const string infilename, bool loadAll)
{
  read(infilename, loadAll);
}

// load constructor

table::table(const string modName, const string modUnits, const Integer nInt,
	     const Integer nAdd, const bool isz, const bool isAdd,
	     const string eUnits, const Real lowElim, const Real highElim,
	     const string filename, const vector<Real>& energies,
	     const vector<tableParameter>& paramObjects,
	     const vector<tableSpectrum>& spectrumObjects, const bool isEscale,
	     const bool isErr)
{
  load(modName, modUnits, nInt, nAdd, isz, isAdd, eUnits, lowElim, highElim,
       filename, energies, paramObjects, spectrumObjects, isEscale, isErr);
}

// destructor

table::~table()
{
  // clear vectors with guaranteed reallocation
  vector<tableParameter>().swap(m_Parameters);
  vector<tableSpectrum>().swap(m_Spectra);
  vector<Real>().swap(m_Energies);
}

// load method

void table::load(const string modName, const string modUnits, const Integer nInt,
		 const Integer nAdd, const bool isz, const bool isAdd,
		 const string eUnits, const Real lowElim, const Real highElim,
		 const string filename, const vector<Real>& energies,
		 const vector<tableParameter>& paramObjects,
		 const vector<tableSpectrum>& spectrumObjects, const bool isEscale,
		 const bool isErr)
{
  m_ModelName = modName;
  m_ModelUnits = modUnits;
  m_NumIntParams = nInt;
  m_NumAddParams = nAdd;
  m_isError = isErr;
  m_isRedshift = isz;
  m_isAdditive = isAdd;
  m_isEscale = isEscale;
  m_EnergyUnits = eUnits;
  m_LowEnergyLimit = lowElim;
  m_HighEnergyLimit = highElim;
  m_Filename = filename;
  m_Energies.resize(energies.size());
  for (size_t i=0; i<energies.size(); i++) m_Energies[i] = energies[i];
  for (size_t i=0; i<paramObjects.size(); i++) pushParameter(paramObjects[i]);
  for (size_t i=0; i<spectrumObjects.size(); i++) pushSpectrum(spectrumObjects[i]);
}

// Read the table from a FITS file. if loadAll is false then don't actualy read
// in the Spectra but set up the objects

Integer table::read(const string infilename, const bool loadAll)
{

  string hduName;
  string DefString;
  bool DefBool;
  Real DefReal;
  Integer DefInteger;

  // open the FITS file at the primary extension

  unique_ptr<FITS> pInfile((FITS*)0);

  try {
    pInfile.reset(new FITS(infilename, Read, false));
  } catch(...) {
    string msg = "Failed to read "+infilename;
    SPreportError(NoSuchFile, msg);
    return(NoSuchFile);
  }

  m_Filename = infilename;

  PHDU& primary = pInfile->pHDU();

  // read the keywords we need from this extension

  m_ModelName = SPreadKey(primary, "MODLNAME", BlankStr);
  DefString = "ph/cm^2/s";
  m_ModelUnits = SPreadKey(primary, "MODLUNIT", DefString);
  DefBool = true;
  m_isAdditive = SPreadKey(primary, "ADDMODEL", DefBool);
  DefBool = false;
  m_isRedshift = SPreadKey(primary, "REDSHIFT", DefBool);
  DefBool = false;
  m_isEscale = SPreadKey(primary, "ESCALE", DefBool);
  DefReal = -1.0;
  m_LowEnergyLimit = SPreadKey(primary, "LOELIMIT", DefReal);
  DefReal = -1.0;
  m_HighEnergyLimit = SPreadKey(primary, "HIELIMIT", DefReal);
  DefInteger = 0;
  Integer numFiltExps = SPreadKey(primary, "NXFLTEXP", DefInteger);
  if ( numFiltExps > 0 ) {
    m_FiltExps.resize(numFiltExps);
    for (size_t i=1; i<=m_FiltExps.size(); i++) {
      ostringstream keyStream;
      keyStream << "XFXP" << setfill('0') << right << setw(4) << i;
      m_FiltExps[i-1] = SPreadKey(primary, keyStream.str(), BlankStr);
    }    
  }

  // at the moment no errors associated with model spectra

  m_isError = false;

  // Move to the parameters extension

  hduName = "PARAMETERS";
  ExtHDU& param = pInfile->extension(hduName);

  // set number of parameters

  m_NumIntParams = SPreadKey(param, "NINTPARM", (Integer)0);
  m_NumAddParams = SPreadKey(param, "NADDPARM", (Integer)0);

  size_t Nparams = m_NumIntParams + m_NumAddParams;

  // read the columns

  if ( Nparams > 0 ) {

    vector<string> name, units;
    vector<Integer> method, numbvals;
    vector<Real> initial, delta, minimum, bottom, top, maximum;
    vector<vector<Real> > value;

    SPreadCol(param, "NAME", name);
    SPreadCol(param, "UNITS", units);
    if ( units.size() == 0 ) {
      units.resize(name.size());
      for (size_t i=0; i<units.size(); i++) units[i] = BlankStr;
    }
    SPreadCol(param, "METHOD", method);
    SPreadCol(param, "INITIAL", initial);
    SPreadCol(param, "DELTA", delta);
    SPreadCol(param, "MINIMUM", minimum);
    SPreadCol(param, "BOTTOM", bottom);
    SPreadCol(param, "TOP", top);
    SPreadCol(param, "MAXIMUM", maximum);
    SPreadCol(param, "NUMBVALS", numbvals);
    SPreadVectorCol(param, "VALUE", value);

    // loop round parameters setting up the objects and loading them
    // note for the name we replace any spaces by underbars to avoid
    // subsequent problems in xspec

    for (size_t iPar=0; iPar<Nparams; iPar++) {
      tableParameter in;
      std::replace(name[iPar].begin(),name[iPar].end(),' ','_');
      in.setName(name[iPar]);
      in.setUnits(units[iPar]);
      in.setInterpolationMethod(method[iPar]);
      in.setInitialValue(initial[iPar]);
      in.setDelta(delta[iPar]);
      in.setMinimum(minimum[iPar]);
      in.setBottom(bottom[iPar]);
      in.setTop(top[iPar]);
      in.setMaximum(maximum[iPar]);
      vector<Real> thisParVals(numbvals[iPar]);
      for (size_t i=0; i<thisParVals.size(); ++i) thisParVals[i]=value[iPar][i];
      in.setTabulatedValues(thisParVals);

      // Force InterpolationMethod=-1 if there are no tabulated parameters
      if ( in.NumberTabulatedValues() == 0 ) in.setInterpolationMethod(-1);

      m_Parameters.push_back(in);
    }

  }

  // Now read the energies

  hduName = "ENERGIES";
  ExtHDU& ener = pInfile->extension(hduName);

  vector<Real> eLow, eHigh;

  SPreadCol(ener, "ENERG_LO", eLow);
  SPreadCol(ener, "ENERG_HI", eHigh);

  size_t Nbins = eLow.size();

  m_Energies.resize(Nbins+1);
  for (size_t i=0; i<Nbins; i++) {
    m_Energies[i] = eLow[i];
  }
  m_Energies[Nbins] = eHigh[Nbins-1];
  m_EnergyUnits = "keV";

  // move to the SPECTRA extension

  hduName = "SPECTRA";
  ExtHDU& spec = pInfile->extension(hduName);

  // get the number of rows

  size_t Nspec = (size_t)SPreadKey(spec, "NAXIS2", (Integer)0);

  // construct the names of any columns with additive parameter spectra

  vector<string> addColName(m_NumAddParams);
  for (size_t iaF=0; iaF<(size_t)m_NumAddParams; iaF++) {
    ostringstream sname;
    sname << "ADDSP" << setfill('0') << setw(3) << iaF+1;
    addColName[iaF] = sname.str();
  }

  // If loading the spectra then Loop round taking each row in the extension at a time

  if ( loadAll ) {

    for (size_t iSpec=0; iSpec<Nspec; iSpec++) {

      tableSpectrum in;

      vector<Real> paramval;
      SPreadVectorColRow(spec, "PARAMVAL", iSpec+1, paramval);

      in.setParameterValues(paramval);

      vector<Real> intpspec;
      SPreadVectorColRow(spec, "INTPSPEC", iSpec+1, intpspec);
      in.setFlux(intpspec);

      // not reading flux errors yet
      intpspec.resize(0);
      in.setFluxError(intpspec);

      // Need to resize addFlux before calling setaddFluxElement
      vector<vector<Real> > tmpaddFlux(m_NumAddParams);
      in.setaddFlux(tmpaddFlux);
      for (size_t iaF=0; iaF<(size_t)m_NumAddParams; iaF++) {

	vector<Real> addpspec;
	SPreadVectorColRow(spec, addColName[iaF], iSpec+1, addpspec);
	in.setaddFluxElement((Integer)iaF, addpspec);
	// not reading flux errors yet
	addpspec.resize(0);
	in.setaddFluxErrorElement((Integer)iaF, addpspec);

      }

      m_Spectra.push_back(in);

    }

  } else {

    // we are not loading the spectra so just set up dummy objects

    for (size_t iSpec=0; iSpec<Nspec; iSpec++) {

      tableSpectrum in;
      m_Spectra.push_back(in);
    
    }

  }

  return(OK);
}

// read the listed Spectra

template <class T> Integer table::readSpectra(T& spectrumList)
{

  // first check whether we have not already loaded the required spectra

  bool done = true;
  for (size_t iList=0; iList<spectrumList.size(); iList++) {
    if ( m_Spectra[spectrumList[iList]].NumberFluxes() == 0 ) {
      done = false;
      break;
    }
  }
  if ( done ) return(OK);

  // open the FITS file using the name which has been recorded in the table object 

  unique_ptr<FITS> pInfile((FITS*)0);

  try {
    pInfile.reset(new FITS(m_Filename, Read, false));
  } catch(...) {
    string msg = "Failed to read "+m_Filename;
    SPreportError(NoSuchFile, msg);
    return(NoSuchFile);
  }

  // move to the SPECTRA extension

  string hduName = "SPECTRA";
  ExtHDU& spec = pInfile->extension(hduName);

  // construct the names of any columns with additive parameter spectra

  vector<string> addColName(m_NumAddParams);
  for (size_t iaF=0; iaF<(size_t)m_NumAddParams; iaF++) {
    ostringstream sname;
    sname << "ADDSP" << setfill('0') << setw(3) << iaF+1;
    addColName[iaF] = sname.str();
  }

  // Loop round reading each row required

  for (size_t iList=0; iList<spectrumList.size(); iList++) {

    size_t iSpec = spectrumList[iList];

    tableSpectrum& in = m_Spectra[iSpec];

    if (in.NumberFluxes() == 0 ) {

      vector<Real> paramval;
      SPreadVectorColRow(spec, "PARAMVAL", iSpec+1, paramval);
      in.setParameterValues(paramval);

      vector<Real> intpspec;
      SPreadVectorColRow(spec, "INTPSPEC", iSpec+1, intpspec);
      in.setFlux(intpspec);

      // not reading flux errors yet
      intpspec.resize(0);
      in.setFluxError(intpspec);

      for (size_t iaF=0; iaF<(size_t)m_NumAddParams; iaF++) {

	vector<Real> addpspec;
	SPreadVectorColRow(spec, addColName[iaF], iSpec+1, addpspec);
	in.setaddFluxElement((Integer)iaF, addpspec);
	// not reading flux errors yet
	addpspec.resize(0);
	in.setaddFluxErrorElement((Integer)iaF, addpspec);

      }

    }

  }

  return(OK);
}

// required for linker instantiation
template Integer table::readSpectra(IntegerArray& spectrumList);
template Integer table::readSpectra(valarray<Integer>& spectrumList);
template Integer table::readSpectra(valarray<size_t>& spectrumList);
template Integer table::readSpectra(vector<size_t>& spectrumList);

// Push table Parameter object

void table::pushParameter(const tableParameter& paramObject)
{
  m_Parameters.push_back(paramObject);
  return;
}

// Push spectrum Parameter object

void table::pushSpectrum(const tableSpectrum& spectrumObject)
{
  m_Spectra.push_back(spectrumObject);
  return;
}

// Get table Parameter object (counts from zero)

tableParameter table::getParameter(const Integer number) const
{
  tableParameter paramObject;
  if ( number >= 0 && number < (Integer)m_Parameters.size() ) paramObject = m_Parameters[number];
  return paramObject;
}

// Get table Spectrum object (counts from zero)

tableSpectrum table::getSpectrum(const Integer number) const
{
  tableSpectrum spectrumObject;
  if ( number >= 0 && number < (Integer)m_Spectra.size() ) spectrumObject = m_Spectra[number];
  return spectrumObject;
}

// deep copy

table& table::copy(const table& a)
{
  setParameters(a.getParameters());
  setSpectra(a.getSpectra());
  setModelName(a.getModelName());
  setModelUnits(a.getModelUnits());
  setNumIntParams(a.getNumIntParams());
  setNumAddParams(a.getNumAddParams());
  setisError(a.getisError());
  setisRedshift(a.getisRedshift());
  setisAdditive(a.getisAdditive());
  setisEscale(a.getisEscale());
  setEnergies(a.getEnergies());
  setEnergyUnits(a.getEnergyUnits());
  setLowEnergyLimit(a.getLowEnergyLimit());
  setHighEnergyLimit(a.getHighEnergyLimit());
  setFilename(a.getFilename());
  setFiltExps(a.getFiltExps());
  return *this;
}

table& table::operator= (const table& a)
{
  return copy(a);
}

// display information about the table - return as a string

string table::disp(const bool headerOnly) const
{
  ostringstream outstr;

  outstr << "Table information : " <<endl;
  outstr << "Read from                          = " << m_Filename << endl;
  outstr << "Model name                         = " << m_ModelName << endl;
  outstr << "Model units                        = " << m_ModelUnits<< endl;
  outstr << "Number of interpolation parameters = " << m_NumIntParams << endl;
  outstr << "Number of additional parameters    = " << m_NumAddParams << endl;
  if ( m_isError ) outstr << "Model contains errors" << endl;
  if ( m_isRedshift ) outstr << "Model includes redshift" << endl;
  if ( m_isEscale ) outstr << "Model includes energy scaling" << endl;
  if ( m_isAdditive ) {
    outstr << "Model is additive" << endl;
  } else {
    outstr << "Model is multiplicative" << endl;
  }
  outstr << "Number of model energies           = " << m_Energies.size() << endl;
  outstr << "Energy units                       = " << m_EnergyUnits << endl;
  if ( m_LowEnergyLimit >= 0.0 )
    outstr << "Model below tabulated energies     = " << m_LowEnergyLimit << endl;
  if ( m_HighEnergyLimit >= 0.0 )
    outstr << "Model above tabulated energies     = " << m_HighEnergyLimit << endl;
  if ( m_FiltExps.size() > 0 ) {
    for (size_t i=0; i<m_FiltExps.size(); i++)
      outstr << "Filter expression " << i+1 << "    = " << m_FiltExps[i] << endl;
  }
  outstr << " " << endl;

  if ( headerOnly ) return outstr.str();

  for (size_t i=0; i<m_Parameters.size(); i++) {
    outstr << "Parameter " << i+1 << " : " << endl;
    outstr << m_Parameters[i].disp() << endl;
  }
  for (size_t i=0; i<m_Spectra.size(); i++) {
    outstr << "Spectrum " << i+1 << " : " << endl;
    outstr << m_Spectra[i].disp() << endl;
  }

  return outstr.str();
}

// clear contents of table object (mainly useful for Python)

void table::clear()
{
  m_Parameters.clear();
  m_Spectra.clear();
  m_ModelName = BlankStr;
  m_ModelUnits = BlankStr;
  m_NumIntParams = 0;
  m_NumAddParams = 0;
  m_isError = false;
  m_isRedshift = false;
  m_isAdditive = false;
  m_isEscale = false;
  m_Energies.clear();
  m_EnergyUnits = BlankStr;
  m_LowEnergyLimit = -1.0;
  m_HighEnergyLimit = -1.0;
  m_Filename = BlankStr;
  return;
}

string table::check() const
{
  ostringstream outstr;

  // check that energies are in increasing order

  bool isIncreasing(true);
  for (size_t i=0; i<m_Energies.size()-1; i++) {
    if ( m_Energies[i] > m_Energies[i+1] ) isIncreasing = false;
  }
  if ( !isIncreasing ) {
    outstr << "The Energies array is not in increasing order" << endl;
  }

  // check consistency of energy arrays. note that Flux may validly have size
  // zero if we have not loaded the complete Spectra objects

  for (size_t i=0; i<m_Spectra.size(); i++) {
    if ( (Integer)m_Energies.size()-1 != m_Spectra[i].NumberFluxes() && m_Spectra[i].NumberFluxes() != 0 ) {
      outstr << "The size of the Energies array (" << m_Energies.size() 
	   << ") is not one more than that for the Flux array (" 
	   << m_Spectra[i].NumberFluxes() << ") for spectrum " << i << endl;
    }
    for (size_t j=0; j<m_Spectra[i].getaddFlux().size(); j++ ) {
      if ( m_Spectra[i].NumberAdditiveFluxes(j) != m_Spectra[i].NumberFluxes() ) {
	outstr << "The size of the addFlux array (" << m_Spectra[i].NumberAdditiveFluxes(j)
	     << ") for additional parameter " << j+1 
	     << " is not equal to that for the Flux array (" 
	     << m_Spectra[i].NumberFluxes() << ") for spectrum in row " << i+1 << endl;
      }
    }
  }

  // check consistency of parameter arrays

  if ( m_NumIntParams + m_NumAddParams != (Integer)m_Parameters.size() ) {
    outstr << "The number of Parameters objects (" << m_Parameters.size() 
	 << ") does not match the sum of interpolated and additional parameters ("
	 << m_NumIntParams+m_NumAddParams << ")" << endl;
  }

  // check that there are the correct number of spectrum objects

  Integer Nspec(1);
  for (size_t i=0; i<(size_t)(m_NumIntParams+m_NumAddParams); i++) {
    if ( m_Parameters[i].getInterpolationMethod() >= 0 && m_Parameters[i].NumberTabulatedValues() > 0 ) {
      Nspec *= m_Parameters[i].NumberTabulatedValues();
    }
  }
  if ( m_FiltExps.size() > 1 ) Nspec *= m_FiltExps.size();
  if ( Nspec != (Integer)m_Spectra.size() ) {
    outstr << "The number of Spectra objects (" << m_Spectra.size()
	 << ") does not match that expected from the Parameters objects ("
	 << Nspec << ")" << endl;
  }

  // check that the spectra all have the same size. Some spectra may validly have
  // zero size if they have not been read in yet.

  Integer Nbins(0);
  for (size_t iSpec=0; iSpec<(size_t)Nspec; iSpec++) {
    if ( m_Spectra[iSpec].NumberFluxes() != 0 ) {
      Nbins = m_Spectra[iSpec].NumberFluxes();
      break;
    }
  }

  for (size_t iSpec=0; iSpec<(size_t)Nspec; iSpec++) {
    if ( m_Spectra[iSpec].NumberFluxes() != Nbins && m_Spectra[iSpec].NumberFluxes() != 0 ) {
      outstr << "The spectrum for row " << iSpec+1 << " is not the same size as the first" << endl;
    }
    for (size_t iaF=0; iaF<m_Spectra[iSpec].getaddFlux().size(); iaF++) {
      if ( m_Spectra[iSpec].NumberAdditiveFluxes(iaF) != Nbins && m_Spectra[iSpec].NumberAdditiveFluxes(iaF) != 0 ) {
	outstr << "The " << iaF+1 << "th additional spectrum for row " 
	       << iSpec+1 << " is not the same size as the first spectrum" << endl;
      }
    }
  }


  return outstr.str();
}

// convert to standard units (keV and ph/cm^2/s). Note that this does not attempt to
// convert LowEnergyLimit or HighEnergyLimit because this cannot be done in all cases.

Integer table::convertUnits()
{
  Real xfactor(1.0);
  bool xwave(false);

  // set up energy/wave conversion factors and check for valid units

  Integer status(OK);

  status = calcXfactor(m_EnergyUnits, xwave, xfactor);
  if ( status != OK ) return(status);

  if ( xfactor != 1.0 ) {
    if ( xwave ) {
      for (size_t i=0; i<m_Energies.size(); i++) {
	m_Energies[i] = xfactor/m_Energies[i];
      }
    } else {
      for (size_t i=0; i<m_Energies.size(); i++) {
	m_Energies[i] *= xfactor;
      }
    }
  }

  m_EnergyUnits = "keV";

  // if the model is not multiplicative then we do not want to do any model
  // flux conversions because the units are just multiplicative factors

  if ( !m_isAdditive ) return(OK);

  bool energy(true);
  bool perwave(false);
  bool perenergy(false);
  Real yfactor(1.0);

  // Find the model units set and check for validity

  status = calcYfactor(m_ModelUnits, energy, perwave, perenergy, yfactor);
  if ( status != OK ) return(status);

  // if nothing more to be done leave now

  if ( yfactor == 1.0 && !energy && !perwave && !perenergy ) return(OK);

  // set up arrays of mean energies/wavelengths and bin sizes
  // both are in keV.

  vector<Real> xmean(m_Energies.size()-1);
  vector<Real> xgmult(m_Energies.size()-1);
  vector<Real> xbinsize(m_Energies.size()-1);

  for (size_t i=0; i<xmean.size(); i++) {
    xmean[i] = (m_Energies[i]+m_Energies[i+1])/2.0;
    xgmult[i] = (m_Energies[i] * m_Energies[i+1]);
    xbinsize[i] = abs(m_Energies[i+1]-m_Energies[i]);
  }

  // Now do the conversions of Spectra.Flux and Spectra.addFlux. Note the six 
  // different possibilities based on whether flux is in photons or energy, 
  // is per wavelength or not, is per energy or not (both per wave and per 
  // energy can be true at the same time).

  for (size_t ispec=0; ispec<m_Spectra.size(); ispec++) {

    vector<Real> thisFlux = m_Spectra[ispec].getFlux();
    for (size_t i=0; i<thisFlux.size(); i++) {
      if ( !energy && !perwave && !perenergy ) {
	// this is the one we have already done - flux must be ph/cm^2/s
      } else if ( !energy && !perwave && perenergy ) {
	thisFlux[i] *= yfactor*xbinsize[i];
      } else if ( !energy && perwave && !perenergy ) {
	thisFlux[i] *= yfactor*xbinsize[i]/xgmult[i];
      } else if ( energy && !perwave && !perenergy ) {
	thisFlux[i] *= yfactor/xmean[i];
      } else if ( energy && !perwave && perenergy ) {
	thisFlux[i] *= yfactor*xbinsize[i]/xmean[i];
      } else if ( energy && perwave && !perenergy ) {
	thisFlux[i] *= yfactor*xbinsize[i]/xmean[i]/xgmult[i];
      }
    }
    m_Spectra[ispec].setFlux(thisFlux);

    for (Integer iadd=0; iadd<m_Spectra[ispec].NumberAdditiveParameters(); iadd++) {
      vector<Real> thisaddFlux = m_Spectra[ispec].getaddFluxElement(iadd);
      for (size_t i=0; i<thisaddFlux.size(); i++) {
	if ( !energy && !perwave && !perenergy ) {
	  // this is the one we have already done - flux must be ph/cm^2/s
	} else if ( !energy && !perwave && perenergy ) {
	  thisaddFlux[i] *= yfactor*xbinsize[i];
	} else if ( !energy && perwave && !perenergy ) {
	  thisaddFlux[i] *= yfactor*xbinsize[i]/xgmult[i];
	} else if ( energy && !perwave && !perenergy ) {
	  thisaddFlux[i] *= yfactor/xmean[i];
	} else if ( energy && !perwave && perenergy ) {
	  thisaddFlux[i] *= yfactor*xbinsize[i]/xmean[i];
	} else if ( energy && perwave && !perenergy ) {
	  thisaddFlux[i] *= yfactor*xbinsize[i]/xmean[i]/xgmult[i];
	}
      }
      m_Spectra[ispec].setaddFluxElement(iadd, thisaddFlux);
    }
  }

  return(OK);
}

// reverse the rows. useful in the case that energies are not in increasing order

void table::reverseRows()
{

  // reverse the energies

  size_t Ne(m_Energies.size());
  vector<Real> TempE(m_Energies);
  for (size_t i=0; i<Ne; i++) m_Energies[i] = TempE[Ne-i-1];

  // loop round the tableSpectrum objects

  for (size_t iSpec=0; iSpec<m_Spectra.size(); iSpec++) {

    // reverse the Flux array

    vector<Real> thisFlux = m_Spectra[iSpec].getFlux();
    vector<Real> TempF(thisFlux.size());
    size_t N(TempF.size());
    for (size_t i=0; i<N; i++) TempF[N-i-1] = thisFlux[i];
    m_Spectra[iSpec].setFlux(TempF);

    // loop over any addFlux vectors

    for (Integer iaF=0; iaF<m_Spectra[iSpec].NumberAdditiveParameters(); iaF++) {

      // reverse this addFlux array
      vector<Real> thisaddFlux = m_Spectra[iSpec].getaddFluxElement(iaF);
      vector<Real> TempaF(thisaddFlux.size());
      for (size_t i=0; i<N; i++) TempaF[N-i-1] = thisaddFlux[i];
      m_Spectra[iSpec].setaddFluxElement(iaF,TempaF);
      
    }

  }

  return;
}

// write to a FITS file

Integer table::write(string outfilename) const
{

  vector<string> ttype;
  vector<string> tform;
  vector<string> tunit;

  // Create a new FITS file instance

  std::unique_ptr<FITS> pFits((FITS*)0);

  try {                
    pFits.reset( new FITS(outfilename,Write) );
  } catch (FITS::CantCreate) {
    string msg = "Failed to create "+outfilename+" for table model file";
    SPreportError(CannotCreate, msg);
    return(CannotCreate);       
  }

  // Write the keywords to the primary header

  PHDU& primary = pFits->pHDU();
  
  SPwriteKey(primary,"HDUCLASS", (string)"OGIP", BlankStr);
  SPwriteKey(primary,"HDUCLAS1", (string)"XSPEC TABLE MODEL", BlankStr);
  SPwriteKey(primary,"HDUVERS", (string)"1.1.0"," ");
  SPwriteKey(primary,"MODLNAME", m_ModelName,"Table model name");
  SPwriteKey(primary,"MODLUNIT", m_ModelUnits,"Table model units");
  SPwriteKey(primary,"REDSHIFT", m_isRedshift,"Add redshift parameter?");
  if ( m_isEscale ) SPwriteKey(primary,"ESCALE", m_isEscale,
			       "Add energy scaling parameter?");
  SPwriteKey(primary,"ADDMODEL", m_isAdditive,"Is model additive?");
  if ( m_LowEnergyLimit >= 0.0 )
    SPwriteKey(primary,"LOELIMIT", m_LowEnergyLimit,
			 "Value of model below tabulated energies");
  if ( m_HighEnergyLimit >= 0.0 )
    SPwriteKey(primary,"HIELIMIT", m_HighEnergyLimit,
			 "Value of model above tabulated energies");
  if ( m_FiltExps.size() > 0 ) {
    SPwriteKey(primary,"NXFLTEXP", m_FiltExps.size(),"Number of XFLT expressions");
    for (size_t i=1; i<=m_FiltExps.size(); i++) {
      ostringstream keyStream;
      keyStream << "XFXP" << setfill('0') << right << setw(4) << i;
      SPwriteKey(primary,keyStream.str(), m_FiltExps[i-1],"XFLT expression");
    }    
  }

  // Set up and create the PARAMETERS extension

  ttype.resize(11);
  tform.resize(11);
  tunit.resize(11);

  ttype[0] = "NAME";
  ttype[1] = "UNITS";
  ttype[2] = "METHOD";
  ttype[3] = "INITIAL";
  ttype[4] = "DELTA";
  ttype[5] = "MINIMUM";
  ttype[6] = "BOTTOM";
  ttype[7] = "TOP";
  ttype[8] = "MAXIMUM";
  ttype[9] = "NUMBVALS";
  ttype[10] = "VALUE";

  tform[0] = "12A";
  tform[1] = "12A";
  tform[2] = "J";
  for (size_t i=3; i<9; i++) tform[i] = "E";
  tform[9] = "J";
  tform[10] = "PE";

  for (size_t i=0; i<11; i++) tunit[i] = BlankStr;

  Table* pparam = pFits->addTable("PARAMETERS",m_Parameters.size(),ttype,tform,tunit);
  Table& param = *pparam;

  SPwriteKey(param,"HDUCLASS", (string)"OGIP", BlankStr);
  SPwriteKey(param,"HDUCLAS1", (string)"XSPEC TABLE MODEL", BlankStr);
  SPwriteKey(param,"HDUCLAS2", (string)"PARAMETERS", BlankStr);
  SPwriteKey(param,"HDUVERS", (string)"1.0.0", BlankStr);
  SPwriteKey(param,"NINTPARM", m_NumIntParams,"Number of interpolation parameters ");
  SPwriteKey(param,"NADDPARM", m_NumAddParams,"Number of additional parameters ");

  // write the parameter info. Note use of valarray because CCfits doesn't
  // support Column::write(vector<T>, Integer).

  size_t Nparams(m_Parameters.size());

  vector<string> names;
  for (size_t ipar=0; ipar<Nparams; ipar++) names.push_back(m_Parameters[ipar].getName());
  param.column("NAME").write(names,1);

  vector<string> units;
  for (size_t ipar=0; ipar<Nparams; ipar++) units.push_back(m_Parameters[ipar].getUnits());
  param.column("UNITS").write(units,1);

  valarray<Integer> ivalues(Nparams);
  for (size_t ipar=0; ipar<Nparams; ipar++) ivalues[ipar] = m_Parameters[ipar].getInterpolationMethod();
  param.column("METHOD").write(ivalues,1);

  valarray<Real> rvalues(Nparams);
  for (size_t ipar=0; ipar<Nparams; ipar++) rvalues[ipar] = m_Parameters[ipar].getInitialValue();
  param.column("INITIAL").write(rvalues,1);

  for (size_t ipar=0; ipar<Nparams; ipar++) rvalues[ipar] = m_Parameters[ipar].getDelta();
  param.column("DELTA").write(rvalues,1);

  for (size_t ipar=0; ipar<Nparams; ipar++) rvalues[ipar] = m_Parameters[ipar].getMinimum();
  param.column("MINIMUM").write(rvalues,1);

  for (size_t ipar=0; ipar<Nparams; ipar++) rvalues[ipar] = m_Parameters[ipar].getBottom();
  param.column("BOTTOM").write(rvalues,1);

  for (size_t ipar=0; ipar<Nparams; ipar++) rvalues[ipar] = m_Parameters[ipar].getTop();
  param.column("TOP").write(rvalues,1);

  for (size_t ipar=0; ipar<Nparams; ipar++) rvalues[ipar] = m_Parameters[ipar].getMaximum();
  param.column("MAXIMUM").write(rvalues,1);

  for (size_t ipar=0; ipar<Nparams; ipar++) ivalues[ipar] = m_Parameters[ipar].NumberTabulatedValues();
  param.column("NUMBVALS").write(ivalues,1);

  // workaround required here because Column::writeArrays will not take
  // vector<vector<T> > only vector<valarray<T> >.

  vector<valarray<Real> > pvalues(Nparams);
  for (size_t ipar=0; ipar<Nparams; ipar++) {
    vector<Real> values = m_Parameters[ipar].getTabulatedValues();
    pvalues[ipar].resize(values.size());
    for (size_t i=0; i<values.size(); i++) pvalues[ipar][i] = values[i];
  }
  param.column("VALUE").writeArrays(pvalues,1);

  // Write/update checksum keywords

  param.writeChecksum();

  // Create the ENERGIES extension

  ttype.resize(2);
  tform.resize(2);
  tunit.resize(2);

  ttype[0] = "ENERG_LO";
  tform[0] = "D";
  tunit[0] = BlankStr;

  ttype[1] = "ENERG_HI";
  tform[1] = "D";
  tunit[1] = BlankStr;

  Table* penergies = pFits->addTable("ENERGIES",m_Energies.size()-1,ttype,tform,tunit);
  Table& energies = *penergies;

  SPwriteKey(energies,"HDUCLASS", (string)"OGIP", BlankStr);
  SPwriteKey(energies,"HDUCLAS1", (string)"XSPEC TABLE MODEL", BlankStr);
  SPwriteKey(energies,"HDUCLAS2", (string)"ENERGIES", BlankStr);
  SPwriteKey(energies,"HDUVERS", (string)"1.0.0", BlankStr);

  // write the energies

  size_t Nenergies(m_Energies.size()-1);
  rvalues.resize(Nenergies);

  for (size_t ien=0; ien<Nenergies; ien++) rvalues[ien] = m_Energies[ien];
  energies.column("ENERG_LO").write(rvalues,1);

  for (size_t ien=0; ien<Nenergies; ien++) rvalues[ien] = m_Energies[ien+1];
  energies.column("ENERG_HI").write(rvalues,1);

  // Write/update checksum keywords

  energies.writeChecksum();
  
  // Create the SPECTRA extension
  // not writing flux errors yet

  ttype.resize(2+m_NumAddParams);
  tform.resize(2+m_NumAddParams);
  tunit.resize(2+m_NumAddParams);

  stringstream RepeatStream;
  RepeatStream << m_NumIntParams;
  string Repeat(RepeatStream.str());

  ttype[0] = "PARAMVAL";
  tform[0] = Repeat+"E";
  tunit[0] = BlankStr;

  RepeatStream.str("");
  RepeatStream << Nenergies;
  Repeat = RepeatStream.str();

  ttype[1] = "INTPSPEC";
  tform[1] = Repeat+"E";
  tunit[1] = BlankStr;

  for (size_t iadd=1; iadd<=(size_t)m_NumAddParams; iadd++) {
    RepeatStream.str("");
    RepeatStream << iadd;
    Repeat = RepeatStream.str();
    if ( iadd < 10 ) {
      ttype[1+iadd] = "ADDSP00"+Repeat;
    } else if ( iadd < 100 ) {
      ttype[1+iadd] = "ADDSP0"+Repeat;
    } else if ( iadd < 1000 ) {
      ttype[1+iadd] = "ADDSP"+Repeat;
    }
    RepeatStream.str("");
    RepeatStream << Nenergies;
    Repeat = RepeatStream.str();
    tform[1+iadd] = Repeat+"E";
    tunit[1+iadd] = BlankStr;
  }

  Table* pspectra = pFits->addTable("SPECTRA",m_Spectra.size(),ttype,tform,tunit);
  Table& spectra = *pspectra;

  SPwriteKey(spectra,"HDUCLASS", (string)"OGIP", BlankStr);
  SPwriteKey(spectra,"HDUCLAS1", (string)"XSPEC TABLE MODEL", BlankStr);
  SPwriteKey(spectra,"HDUCLAS2", (string)"MODEL SPECTRA", BlankStr);
  SPwriteKey(spectra,"HDUVERS", (string)"1.0.0", BlankStr);

  // Write the spectra

  size_t Nspectra(m_Spectra.size());
  vector<valarray<Real> > rarray(Nspectra);

  if ( m_NumIntParams > 1 ) {
    for (size_t isp=0; isp<Nspectra; isp++) {
      vector<Real> values = m_Spectra[isp].getParameterValues();
      rarray[isp].resize(values.size());
      for (size_t i=0; i<values.size(); i++) rarray[isp][i] = values[i];
    }
    spectra.column("PARAMVAL").writeArrays(rarray,1);
  } else if ( m_NumIntParams > 0 ) {
    rvalues.resize(Nspectra);
    for (size_t isp=0; isp<Nspectra; isp++) rvalues[isp] = m_Spectra[isp].getParameterValuesElement(0);
    spectra.column("PARAMVAL").write(rvalues,1);
  }

  if ( Nenergies > 1 ) {
    for (size_t isp=0; isp<Nspectra; isp++) {
      vector<Real> thisFlux = m_Spectra[isp].getFlux();
      rarray[isp].resize(thisFlux.size());
      for (size_t i=0; i<thisFlux.size(); i++) rarray[isp][i] = thisFlux[i];
    }
    spectra.column("INTPSPEC").writeArrays(rarray,1);

    for (size_t iadd=1; iadd<=(size_t)m_NumAddParams; iadd++) {
      for (size_t isp=0; isp<Nspectra; isp++) {
	vector<Real> thisaddFlux = m_Spectra[isp].getaddFluxElement(iadd-1);
	rarray[isp].resize(thisaddFlux.size());
	for (size_t i=0; i<thisaddFlux.size(); i++) rarray[isp][i] = thisaddFlux[i];
      }
      spectra.column(ttype[iadd+1]).writeArrays(rarray,1);
    }
  } else {
    rvalues.resize(Nspectra);
    for (size_t isp=0; isp<Nspectra; isp++) rvalues[isp] = m_Spectra[isp].getFluxElement(0);
    spectra.column("INTPSPEC").write(rvalues,1);

    for (size_t iadd=1; iadd<=(size_t)m_NumAddParams; iadd++) {
      for (size_t isp=0; isp<Nspectra; isp++) rvalues[isp] = m_Spectra[isp].getaddFluxElementElement(iadd-1,0);
      spectra.column(ttype[iadd+1]).write(rvalues,1);
    }

  }  

  // Write/update checksum keywords

  spectra.writeChecksum();

  return(OK);

}

// get values from table for input parameters using interpolation.

template <class T> Integer table::getValues(const T& parameterValues, const Real minEnergy, 
					    const Real maxEnergy, T& tableEnergyBins, 
					    T& tableValues, T& tableErrors)
{
  const std::map<string,Real> xfltInfo;
  return getValues(parameterValues, xfltInfo, minEnergy, maxEnergy, tableEnergyBins,
		   tableValues, tableErrors);
}

// now the version with xflt info

template <class T> Integer table::getValues(const T& parameterValues, 
					    const std::map<string,Real>& xfltInfo,
					    const Real minEnergy, 
					    const Real maxEnergy, T& tableEnergyBins, 
					    T& tableValues, T& tableErrors)
{
  // if redshift is included then it will be the last entry in parameterValues so
  // ignore that for the moment but store the redshift factor

  size_t Nparin(parameterValues.size());
  Real zfact(1.0);
  if ( m_isRedshift) {
    zfact = 1.0 + parameterValues[Nparin-1];
    Nparin--;
  }

  // if escale is included then it will now be the last entry in parameterValues so
  // ignore that for the moment but store the escale.

  Real escale(1.0);
  if ( m_isEscale ) {
    escale = parameterValues[Nparin-1];
    Nparin--;
  }

  if ( Nparin != (size_t)(m_NumIntParams+m_NumAddParams) ) {
    return(InconsistentNumTableParams);
  }

  // test for a filter specification. at present the only filter
  // specification allowed is a keyword name. setFilterNum is then
  // set to the value of that keyword in xfltInfo to specify which
  // model spectrum is required for a given parameter vector.

  size_t numFilters = m_FiltExps.size();
  size_t selFilterNum(0);
  if ( numFilters > 0 && xfltInfo.size() == 0 ) return(InconsistentTableFilter);

  // at the moment we assume that all the m_FiltExps use the same keyword
  // so just need to get it from the first one
  // extract the keyword and value for this filter
  if (numFilters)
  {
     size_t colPos = m_FiltExps[0].find(":");
     if ( colPos == string::npos ) return(InconsistentTableFilter);
     string key = m_FiltExps[0].substr(0,colPos);
     // search for the same keyword in the input xfltInfo. if it is found then
     // set selFilterNum to its value. if it isn't then we have a problem so
     // return an error
     if ( xfltInfo.find(key) == xfltInfo.end() ) return(InconsistentTableFilter);
     selFilterNum = (size_t)xfltInfo.find(key)->second;
  }
  // set up a temporary version of m_Energies with energy scaling

  vector<Real> scaledEnergies(m_Energies);
  size_t nEnergies = scaledEnergies.size();
  if ( m_isEscale ) {
    for (size_t i=0; i< nEnergies; i++) scaledEnergies[i] *= escale;
  }

  // test for the case of no overlap between minEnergy to maxEnergy and the
  // tabulated energies. In this case return the minimum and maximum
  // table energies in tableEnergy and either m_LowEnergyLimit or 
  // m_HighEnergyLimit in tableValues as appropriate.

  if ( maxEnergy < scaledEnergies[0] || minEnergy > scaledEnergies[nEnergies-1] ) {
    tableEnergyBins.resize(2);
    tableEnergyBins[0] = scaledEnergies[0];
    tableEnergyBins[1] = scaledEnergies[nEnergies-1];
    tableValues.resize(1);
    if ( maxEnergy < scaledEnergies[0] ) {
      tableValues[0] = m_LowEnergyLimit;
    } else {
      tableValues[0] = m_HighEnergyLimit;
    }
    tableErrors.resize(0);
    return(0);
  }

  // find the range of tabulated energies between minEnergy and maxEnergy.
  // minEindex should be the last entry in Energies below minEnergy and
  // maxEindex should be the first entry in Energies above maxEnergy.
  // minEnergy and maxEnergy are assumed to be in observed frame while Energies
  // are in source frame so need to take into account (1+z) factor.
  // if maxEnergy < minEnergy then use all the tabulated energies

  size_t minEindex = 0;
  size_t maxEindex = nEnergies-1;
/*  if ( minEnergy <= maxEnergy ) {
    while ( scaledEnergies[minEindex] < minEnergy*zfact ) minEindex++;
    if ( minEindex > 0 ) minEindex--;
    while ( scaledEnergies[maxEindex] > maxEnergy*zfact ) maxEindex--;
    if ( maxEindex < scaledEnergies.size()-1 ) maxEindex++;
  }
*/
  const size_t NE(maxEindex-minEindex+1);

  // copy the table energies into the output array if the redshift is non-zero
  // then the output tableEnergyBins will be in the observed frame

  tableEnergyBins.resize(NE);
  if ( zfact == 1.0 ) {
    for (size_t ie=0; ie<NE; ie++) tableEnergyBins[ie] = scaledEnergies[ie+minEindex];
  } else {
    for (size_t ie=0; ie<NE; ie++) tableEnergyBins[ie] = scaledEnergies[ie+minEindex]/zfact;
  }

  // now split out interpolation from additional parameters
  vector<Real> interParamValues, addParamValues;
  vector<Integer> interParamIndex, addParamIndex;
  for (size_t i=0; i<Nparin; i++) {
    if ( m_Parameters[i].getInterpolationMethod() != -1 ) {
      interParamValues.push_back(parameterValues[i]);
      interParamIndex.push_back(i);
    } else {
      addParamValues.push_back(parameterValues[i]);
      addParamIndex.push_back(i);
    }
  }

  size_t Ninter = interParamValues.size();
  if ( Ninter != (size_t)m_NumIntParams) {
    return(InconsistentNumTableParams);
  }
  const size_t NA(m_NumAddParams);
  
  // special case of no interpolation parameters in which case there is only
  // one tableSpectrum

  if ( Ninter == 0 ) {

    const size_t NF = NE-1;
    tableValues.resize(NF);
    tableSpectrum& tabSpec = m_Spectra[0];
    vector<Real> thisFlux = tabSpec.getFlux();
    for (size_t ie=0; ie<NF; ie++) tableValues[ie] = thisFlux[ie+minEindex];
    if ( m_isError ) {
      tableErrors.resize(NF);
      vector<Real> thisFluxError = tabSpec.getFluxError();
      for (size_t ie=0; ie<NF; ie++) {
	tableErrors[ie] = thisFluxError[ie+minEindex]*thisFluxError[ie+minEindex];
      }
    }

    // accumulate additional parameter spectra
    for (size_t k=0; k<NA; k++) {
      Real parValue = addParamValues[k];
      vector<Real> thisaddFlux = tabSpec.getaddFluxElement(k);
      for (size_t ie=0; ie<NF; ie++) tableValues[ie] += thisaddFlux[ie+minEindex] * parValue;
      if ( m_isError ) {
	vector<Real> thisaddFluxError = tabSpec.getaddFluxErrorElement(k);
	for (size_t ie=0; ie<NF; ie++) tableErrors[ie] += pow(2.0,thisaddFluxError[ie+minEindex] * parValue);
      }
    }

    // include time dilation factor if required and set tableErrors to be sigmas
    if ( zfact != 1.0 && getisAdditive() ) {
      for (size_t ie=0; ie<NF; ie++) tableValues[ie] /= zfact;
    }
    if ( m_isError ) {
      if ( zfact == 1.0 || !getisAdditive() ) {
	for (size_t ie=0; ie<NF; ie++) tableErrors[ie] = sqrt(tableErrors[ie]);
      } else {
	for (size_t ie=0; ie<NF; ie++) tableErrors[ie] = sqrt(tableErrors[ie])/zfact;
      }
    }

    return(OK);
    
  }

  // Ninter is non-zero so now, test for limits.
  size_t totalBlockSize(1);
  vector<IntegerArray> prepRecordNumbers(Ninter);
  vector<bool> exactMatch(Ninter);
  for (size_t j= 0; j < Ninter; ++j) {
    Real value = interParamValues[j];
    tableParameter& tabParam = m_Parameters[interParamIndex[j]];
    vector<Real> tabValue = tabParam.getTabulatedValues();
    size_t N = tabValue.size();
    // accumulate product of numbers of parameter values,
    // will be used to calculate table offsets below.
    totalBlockSize *= N;
    size_t kExact (0);
    if ( N > 1 ) {
      // Condition for out-of-bounds test must be consistent with
      // the test for exactness below.  Note that tabValue's original
      // input comes from floats in a FITS file, not doubles.
      const Real fuzz = (value == 0.0) ? 0.0 : FUZZY;
      const Real magnitude = (value == 0.0) ? 1.0 : std::abs(value);
      if ((tabValue[0] - value)/magnitude > fuzz || 
	  (value - tabValue[N-1])/magnitude > fuzz ) {
	ostringstream outstr;
	outstr << "Requested parameter value " << value
	       << " outside tabulated range (" << tabValue[0] << ", "
	       << tabValue[N-1] << ")";
	SPreportError(TableParamValueOutsideRange, outstr.str());
 	return(TableParamValueOutsideRange);
      } else {
	for (; kExact < N-1; ++kExact) {
	  if ( (exactMatch[j] 
		= (std::abs((value - tabValue[kExact])/magnitude) <= fuzz)) )
	    break;
	}       
      }
    } else {
      exactMatch[j] = true;
    }

    if (exactMatch[j]) {
      // just store this here for now, will resize
      // correctly later.
      prepRecordNumbers[j].resize(1,kExact);
    } else {
      size_t n (1);
      // 2^(j+1) records needed.
      for (size_t l = 0; l <= j ; ++l ) if (!exactMatch[l]) n *= 2; 
      prepRecordNumbers[j].resize(n,0);
    }
  }
  
  IntegerArray blockOffset(Ninter,1);
  blockOffset[0] = totalBlockSize;

  size_t jj(0);
  do {
    totalBlockSize /= m_Parameters[interParamIndex[jj]].NumberTabulatedValues(); 
    blockOffset[jj] = totalBlockSize;    
    ++jj;  
  } while ( jj < Ninter - 1);
  
  IntegerArray bracket(Ninter);
  for (size_t j = 0; j < Ninter; ++ j) {
    Real value = interParamValues[j];
    tableParameter& tabParam = m_Parameters[interParamIndex[j]];
    vector<Real> tabValue = tabParam.getTabulatedValues();
    if ( !exactMatch[j] ) {
      // bracket[j] is the arraypoint in TabulatedValues below the target
      // value, which is therefore straddled by (bracket[j],bracket[j]+1).
      // if the range is exceed, perform constant extrapolation.
      SPfind(tabValue,value,bracket[j]);
      if ( bracket[j] >= static_cast<int>(tabValue.size() - 1)) {
	exactMatch[j] = true;
	prepRecordNumbers[j][0] = tabValue.size() - 1;   
	for ( size_t l  = j ;  l < Ninter; ++l) {
	  prepRecordNumbers[l].resize(prepRecordNumbers[l].size()/2,0);   
	}
      } else if ( bracket[j] < 0) {
	exactMatch[j] = true;
	prepRecordNumbers[j][0] = 0;   
	for ( size_t l  = j ;  l < Ninter; ++l) {
	  prepRecordNumbers[l].resize(prepRecordNumbers[l].size()/2,0);   
	}      
      } else {
	if ( j == 0 ) {
	  prepRecordNumbers[0][0] = blockOffset[0]*(bracket[0]);   
	  prepRecordNumbers[0][1] = blockOffset[0]*(bracket[0] + 1);  
	} else {
	  IntegerArray& previous = prepRecordNumbers[j-1];
	  size_t MP = previous.size();
	  for (size_t k = 0; k < MP; ++k) {
	    prepRecordNumbers[j][2*k] = previous[k]  + blockOffset[j]*(bracket[j]);
	    prepRecordNumbers[j][2*k + 1] 
	      = previous[k]  + blockOffset[j]*(bracket[j] + 1);
	  }
	}
	
      }
    }  
    
    if (exactMatch[j]) {
      Integer kExact = prepRecordNumbers[j][0];
      if ( j == 0 ) {
	// this will evaluate correctly to the #of blocks
	// before the one containing the record of interest
	// because kExact is 0 based. 
	prepRecordNumbers[j][0] = blockOffset[0]*kExact;
      } else {
	IntegerArray& previous = prepRecordNumbers[j-1];
	prepRecordNumbers[j].resize(previous.size());
	size_t MP = previous.size();
	for (size_t k = 0; k < MP; ++k) {
	  prepRecordNumbers[j][k] = previous[k] + 
	    blockOffset[j]*kExact;
	}
      }
    }     
  }

  // now for each interpolation parameter i the input values lies between the
  // bracket[i] and bracket[i+1] entries in TabulatedValues. recordNumbers points
  // to the entries in Spectra which will be required. However this assumes so
  // far that there is one spectrum per parameter vector. If filters are in use
  // there will be m_FiltExps.size() for each and we need the one specified by
  // selFilterNum.

  IntegerArray recordNumbers(prepRecordNumbers[Ninter-1]); 
  size_t numFiltExps(m_FiltExps.size());
  if ( numFiltExps > 1 ) {
    for (size_t i=0; i<recordNumbers.size(); i++) {
      recordNumbers[i] *= numFiltExps;
      recordNumbers[i] += selFilterNum;
    }
  }

  // calculate the fractions for each interpolation
  RealArray fraction(0.0, interParamValues.size());
  for (size_t j=0; j<fraction.size(); ++j) {
    if (!exactMatch[j]) {
      Real parValue = interParamValues[j];
      tableParameter& tabParam = m_Parameters[interParamIndex[j]];
      vector<Real> tParValues = tabParam.getTabulatedValues();
      Real x1 = tParValues[bracket[j]];
      Real x2 = tParValues[bracket[j]+1];
      if (tabParam.getInterpolationMethod() == 0) {
	fraction[j] = (parValue - x1)/(x2 - x1);
      } else {
	// we know x1 < parVal < x2.
	// now, we ought to check that x1,x2 > 0 earlier
	// than this point!
	fraction[j] = log(parValue/x1)/log(x2/x1);           
      }
    }
  }
  
  // now set up to do the interpolation

  const size_t NP(interParamValues.size());
  const size_t NR(recordNumbers.size());

  // make sure that we have read in the Spectra objects which we will be using
  readSpectra(recordNumbers);

  // load all the spectra and variances into work arrays
  vector<RealArray> spectrumEntries(NR);
  vector<RealArray> varianceEntries(NR);

  const size_t NF = NE-1;
  for (size_t i=0; i<NR; i++) {
    tableSpectrum& tabSpec = m_Spectra[recordNumbers[i]];
    spectrumEntries[i].resize(NF);
    vector<Real> thisFlux = tabSpec.getFlux();
    for (size_t ie=0; ie<NF; ie++) spectrumEntries[i][ie] = thisFlux[ie+minEindex];
    if ( m_isError ) {
      varianceEntries[i].resize(NF);
      vector<Real> thisFluxError = tabSpec.getFluxError();
      for (size_t ie=0; ie<NF; ie++) {
	varianceEntries[i][ie] = thisFluxError[ie+minEindex]*thisFluxError[ie+minEindex];
      }
    }
  }

  // accumulate additional parameter spectra for each recordNumber
  for (size_t j=0; j<NR; j++) {
    tableSpectrum& tabSpec = m_Spectra[recordNumbers[j]];
    for (size_t k=0; k<NA; k++) {
      Real parValue = addParamValues[k];
      vector<Real> thisaddFlux = tabSpec.getaddFluxElement(k);
      for (size_t ie=0; ie<NF; ie++) spectrumEntries[j][ie] += thisaddFlux[ie+minEindex] * parValue;
      if ( m_isError ) {
	vector<Real> thisaddFluxError = tabSpec.getaddFluxErrorElement(k);
	for (size_t ie=0; ie<NF; ie++) varianceEntries[j][ie] += pow(2.0,thisaddFluxError[ie+minEindex] * parValue);
      }
    }
  } 

  // and the interpolation

  if ( NR > 1 ) {

    size_t nd = NR;
    vector<RealArray> reducedSp(nd);
    vector<RealArray> reducedVar(nd);
    for (size_t i=0; i<nd; i++) reducedSp[i].resize(NF);
    if ( m_isError ) {
      for (size_t i=0; i<nd; i++) reducedVar[i].resize(NF);
    }
    for (Integer j=NP-1; j>=0; --j ) {
      Real factor = fraction[j];
      Real complfactor = 1.0-fraction[j];
      if ( !exactMatch[j] ) {
	nd /= 2;
	for (size_t k=0; k<nd; ++k) {
	  const size_t k2 = 2*k;
	  reducedSp[k] = complfactor*spectrumEntries[k2] + factor*spectrumEntries[k2+1];
	  if ( m_isError ) {
	    reducedVar[k] = complfactor*varianceEntries[k2]
		   + factor*varianceEntries[k2+1];       
	  }  
	}
      } else {
	for ( size_t k = 0; k < nd; ++k ) reducedSp[k] = spectrumEntries[k];
	if ( m_isError ) {
	  for (size_t k=0; k<nd; ++k) {
	    reducedVar[k] = varianceEntries[k];
	  }
	}
      }
      for (size_t i=0; i<nd; ++i) spectrumEntries[i] = reducedSp[i];
      if ( m_isError ) {
	for (size_t i=0; i<nd; ++i) varianceEntries[i] = reducedVar[i];
      }

    }

    if (nd != 1 ) {
      // we should never end up here
    }

  }

  // set the output arrays including time dilation factor if redshift is non-zero

  tableValues.resize(NF);
  if ( zfact == 1.0 || !getisAdditive() ) {
    for (size_t ie=0; ie<NF; ie++) tableValues[ie] = spectrumEntries[0][ie];
  } else {
    for (size_t ie=0; ie<NF; ie++) tableValues[ie] = spectrumEntries[0][ie]/zfact;
  }
  if ( m_isError ) {
    tableErrors.resize(NF);
    if ( zfact == 1.0 || !getisAdditive() ) {
      for (size_t ie=0; ie<NF; ie++) tableErrors[ie] = sqrt(varianceEntries[0][ie]);
    } else {
      for (size_t ie=0; ie<NF; ie++) tableErrors[ie] = sqrt(varianceEntries[0][ie])/zfact;
    }
  }

  return(OK);
}

// required for linker instantiation
template Integer table::getValues(const RealArray&, const Real, const Real,
				  RealArray&, RealArray&, RealArray&);
template Integer table::getValues(const vector<Real>&, const Real, const Real,
				  vector<Real>&, vector<Real>&, vector<Real>&);
template Integer table::getValues(const RealArray&, const std::map<string,Real>&,
				  const Real, const Real,
				  RealArray&, RealArray&, RealArray&);
template Integer table::getValues(const vector<Real>&, const std::map<string,Real>&,
				  const Real, const Real,
				  vector<Real>&, vector<Real>&, vector<Real>&);


//-------------------------------------------------------------------------------
// Class tableParameter

// default constructor

tableParameter::tableParameter()
  : m_Name(BlankStr),
    m_Units(BlankStr),
    m_InterpolationMethod(0),
    m_InitialValue(0.0),
    m_Delta(0.01),
    m_Minimum(0.0),
    m_Bottom(0.0),
    m_Top(1.0e6),
    m_Maximum(1.0e6),
    m_TabulatedValues()
{
}

// load constructor

tableParameter::tableParameter(const string name, string units, 
			       const Integer interp,
			       const Real initial, const Real delta,
			       const Real min, const Real bot, const Real top,
			       const Real max, const vector<Real>& values)
{
  load(name, units, interp, initial, delta, min, bot, top, max, values);
}

// backwards-compatible load constructor without the units

tableParameter::tableParameter(const string name, const Integer interp,
			       const Real initial, const Real delta,
			       const Real min, const Real bot, const Real top,
			       const Real max, const vector<Real>& values)
{
  string units = BlankStr;
  load(name, units, interp, initial, delta, min, bot, top, max, values);
}

// destructor

tableParameter::~tableParameter()
{
  // clear vector with guaranteed reallocation
  vector<Real>().swap(m_TabulatedValues);
}

// load method

void tableParameter::load(const string name, string units, 
			  const Integer interp, const Real initial,
			  const Real delta, const Real min, const Real bot, 
			  const Real top, const Real max, 
			  const vector<Real>& values)
{
  m_Name = name;
  m_Units = units;
  m_InterpolationMethod = interp;
  m_InitialValue = initial;
  m_Delta = delta;
  m_Minimum = min;
  m_Bottom = bot;
  m_Top = top;
  m_Maximum = max;
  m_TabulatedValues.resize(values.size());
  for (size_t i=0; i<values.size(); i++) m_TabulatedValues[i] = values[i];
}

// backwards-compatible load method without the units

void tableParameter::load(const string name,  
			  const Integer interp, const Real initial,
			  const Real delta, const Real min, const Real bot, 
			  const Real top, const Real max, 
			  const vector<Real>& values)
{
  string units = BlankStr;
  load(name, units, interp, initial, delta, min, bot, top, max, values);
}

// display information about the table parameter - return as a string

string tableParameter::disp() const
{
  ostringstream outstr;

  outstr << "Parameter information : " << endl;
  outstr << "Parameter name       = " << m_Name << endl;
  outstr << "Parameter units      = " << m_Units << endl;
  outstr << "Interpolation method = ";
  if ( m_InterpolationMethod == 0 ) {
    outstr << "Linear interpolation" << endl;
  } else if ( m_InterpolationMethod == 1 ) {
    outstr << "Logarithmic interpolation" << endl;
  } else if ( m_InterpolationMethod == -1 ) {
    outstr << "Additional (non-interpolated)" << endl;
  } else {
    outstr << "Unrecognized interpolation method" << endl;
  }

  outstr << "Initial value        = " << m_InitialValue << endl;
  outstr << "Delta                = " << m_Delta << endl;
  outstr << "Minimum              = " << m_Minimum << endl;
  outstr << "Bottom               = " << m_Bottom << endl;
  outstr << "Top                  = " << m_Top << endl;
  outstr << "Maximum              = " << m_Maximum << endl;
  if ( m_InterpolationMethod != -1 ) {
    outstr << "Tabulated values     = ";
    for (size_t i=0; i<m_TabulatedValues.size(); i++) outstr << m_TabulatedValues[i] << "  ";
  }
  outstr << endl;

  return outstr.str();
}

// clear contents of the table parameter (mainly useful for Python)

void tableParameter::clear()
{
  m_Name = BlankStr;
  m_Units = BlankStr;
  m_InterpolationMethod = 0;
  m_InitialValue = 0.0;
  m_Delta = 0.0;
  m_Minimum = 0.0;
  m_Bottom = 0.0;
  m_Top = 0.0;
  m_Maximum = 0.0;
  m_TabulatedValues.clear();
  return;
}

//-------------------------------------------------------------------------------
// Class tableSpectrum

// default constructor

tableSpectrum::tableSpectrum()
  : m_Flux(),
    m_FluxError(),
    m_ParameterValues(),
    m_addFlux(),
    m_addFluxError()
{
}

// load constructor

tableSpectrum::tableSpectrum(const vector<Real>& parVals, const vector<Real>& flux,
			     const vector<vector<Real> >& addFlux,
			     const vector<Real>& fluxErr,
			     const vector<vector<Real> >& addFluxErr)
{
  load(parVals, flux, addFlux, fluxErr, addFluxErr);
}

// destructor

tableSpectrum::~tableSpectrum()
{
  // clear vectors with guaranteed reallocation
  vector<Real>().swap(m_Flux);
  vector<Real>().swap(m_ParameterValues);
  for (size_t i=0; i<getaddFlux().size(); i++) vector<Real>().swap(m_addFlux[i]);
  vector<vector<Real> >().swap(m_addFlux);

}

// load method

void tableSpectrum::load(const vector<Real>& parVals, const vector<Real>& flux,
			 const vector<vector<Real> >& addFlux,
			 const vector<Real>& fluxErr,
			 const vector<vector<Real> >& addFluxErr)
{
  m_ParameterValues.resize(parVals.size());
  for (size_t i=0; i<parVals.size(); i++) m_ParameterValues[i] = parVals[i];
  m_Flux.resize(flux.size());
  for (size_t i=0; i<flux.size(); i++) m_Flux[i] = flux[i];
  m_addFlux.resize(addFlux.size());
  for (size_t i=0; i<addFlux.size(); i++) {
    m_addFlux[i].resize(addFlux[i].size());
    for (size_t j=0; j<addFlux[i].size(); j++) m_addFlux[i][j] = addFlux[i][j];
  }
  m_FluxError.resize(fluxErr.size());
  for (size_t i=0; i<fluxErr.size(); i++) m_FluxError[i] = fluxErr[i];
  m_addFluxError.resize(addFluxErr.size());
  for (size_t i=0; i<addFluxErr.size(); i++) {
    m_addFluxError[i].resize(addFluxErr[i].size());
    for (size_t j=0; j<addFluxErr[i].size(); j++) m_addFluxError[i][j] = addFluxErr[i][j];
  }
}
  
// push an additional parameter spectrum

void tableSpectrum::pushaddFlux(vector<Real> input)
{
  m_addFlux.push_back(input);
  return;
}

// Display information about the table spectrum - return as a string

string tableSpectrum::disp() const
{
  ostringstream outstr;

  outstr << "Spectrum information : " << endl;
  outstr << "Number of model flux bins        = " << m_Flux.size() << endl;
  outstr << "Parameter values                 = ";
  for (size_t i=0; i<m_ParameterValues.size(); i++) outstr << m_ParameterValues[i] << "  ";
  outstr << endl;
  outstr << "Number of additional flux arrays = " << m_addFlux.size();
  return outstr.str();
}

// clear contents of the table parameter (mainly useful for Python)

void tableSpectrum::clear()
{
  m_Flux.clear();
  m_ParameterValues.clear();
  m_addFlux.clear();
  return;
}