#include <cfortran.h>
#include <xsTypes.h>
#include <XSstreams.h>
#include <XSModel/Model/UniqueEnergy.h>
#include <XSModel/Model/Component/OGIPTable/OGIPTable.h>
#include <XSModel/Parameter/TableModParam.h>
#include <XSUtil/Utils/IosHolder.h>
#include <XSUtil/Utils/XSutility.h>
#include <algorithm>

extern "C" void xsmtbl(float* ear, int ne, float* param, const char* filenm, int ifl, 
                float* photar, float* photer);
extern "C" void xsatbl(float* ear, int ne, float* param, const char* filenm, int ifl, 
                float* photar, float* photer);

FCALLSCSUB7(xsmtbl,XSMTBL,xsmtbl,FLOATV,INT,FLOATV,STRING,INT,FLOATV,FLOATV)
FCALLSCSUB7(xsatbl,XSATBL,xsatbl,FLOATV,INT,FLOATV,STRING,INT,FLOATV,FLOATV)

void xsmtbl(float* ear, int ne, float* param, const char* filenm, int ifl, 
                float* photar, float* photer)
{
   static OGIPTable table("",0);

   string fileName(filenm);
   bool isExp = false;
   try
   {
      if (fileName != table.filename())
      {      
         table.filename(fileName);
         table.read(true,0,false);
      }
      RealArray energyArray(.0,ne+1);
      RealArray fluxArray(.0,ne);
      RealArray fluxErrArray(.0,ne);
      std::copy(&ear[0], &ear[ne+1], &energyArray[0]);

      UniqueEnergy unEng(energyArray);
      // This is probably not necessary.
      unEng.addClient(ifl);
      // In standard xspec usage, table parameters are added to
      // global parameter lists from which their parameters can
      // be adjusted.  These however are only local, so xsmtbl must
      // do all parameter setting.
      const std::vector<Parameter*>& tablePars = table.parameters();
      const size_t nPars = tablePars.size();
      for (size_t i=0; i<nPars; ++i)
      {
         // Just set using default 'v'.  Presumably any adjustments
         // due to limits have already been performed on incoming
         // param array.
         ModParam* tablePar = dynamic_cast<ModParam*>(tablePars[i]);
         tablePar->setValue(param[i]);
      }
      table.energyWeights(&unEng);
      table.interpolate(&unEng, fluxArray, fluxErrArray, false, isExp);
      std::copy(&fluxArray[0], &fluxArray[0]+ne, &photar[0]);
      std::copy(&fluxErrArray[0], &fluxErrArray[0]+ne, &photer[0]);
   }
   catch (...)
   {
      *IosHolder::errHolder() <<"***Error in xsmtbl function."<<std::endl;
   }

} 


void xsatbl(float* ear, int ne, float* param, const char* filenm, int ifl, 
                float* photar, float* photer)
{
   static OGIPTable table("",0);

   string fileName(filenm);
   bool isExp = false;
   try
   {
      if (fileName != table.filename())
      {      
         table.filename(fileName);
         table.read(true,0,true);
      }
      RealArray energyArray(.0,ne+1);
      RealArray fluxArray(.0,ne);
      RealArray fluxErrArray(.0,ne);
      std::copy(&ear[0], &ear[ne+1], &energyArray[0]);

      UniqueEnergy unEng(energyArray);
      // This is probably not necessary.
      unEng.addClient(ifl);
      // In standard xspec usage, table parameters are added to
      // global parameter lists from which their parameters can
      // be adjusted.  These however are only local, so xsatbl must
      // do all parameter setting.
      const std::vector<Parameter*>& tablePars = table.parameters();
      const size_t nPars = tablePars.size();
      for (size_t i=0; i<nPars; ++i)
      {
         // Just set using default 'v'.  Presumably any adjustments
         // due to limits have already been performed on incoming
         // param array.
         ModParam* tablePar = dynamic_cast<ModParam*>(tablePars[i]);
         tablePar->setValue(param[i]);
      }
      table.energyWeights(&unEng);
      table.interpolate(&unEng, fluxArray, fluxErrArray, true, isExp);

      // time dilation. Check for a redshift parameter and apply redshifts to
      // the energy and variance array.

      XSutility::MatchPtrName<Parameter> matchName;

      std::vector<Parameter*>::const_iterator itPar =
              find_if(tablePars.begin(),tablePars.end(),bind2nd(matchName,"z"));

      Real z (0);      
      if ( itPar != tablePars.end() ) z = (*itPar)->value();

      if ( z > 0 ) 
      {
              Real zf (1/(1+z));
              fluxArray *= zf;
              fluxErrArray *= zf*zf;
      }

      std::copy(&fluxArray[0], &fluxArray[0]+ne, &photar[0]);
      std::copy(&fluxErrArray[0], &fluxErrArray[0]+ne, &photer[0]);
   }
   catch (...)
   {
      *IosHolder::errHolder() <<"***Error in xsatbl function."<<std::endl;
   }
}