#include <XSUtil/FunctionUtils/FunctionUtility.h>
#include <XSUtil/FunctionUtils/xsFortran.h>
#include <sstream>

extern "C" void smdem0_(int& itype, float* ear, int& ne, float* nabun, 
                        float& dens, float& z, int& ninputt, float* inputt,
                        float* dem, int& ifl, float* photar, int& status);

extern "C" void smdem1_(int& itype, float* ear, int& ne, float* nabun, 
                        float& dens, float& z, int& ninputt, float* inputt,
                        float* dem, int& ifl, float* photar, int& status);

extern "C" void smdem2_(int& itype, float* ear, int& ne, float* nabun, 
                        float& dens, float& z, int& ninputt, float* inputt,
                        float* dem, int& ifl, int& qtherm, float& velocity, 
                        float* photar, int& status);

// This strange looking interface is for the benefit of Fortran functions which
// call it.  Too many args for cfortran.
extern "C" void sumdem_(int* pItype, int* pSwtch, float* ear, int* pNe, float* abun,
                        float* pDens, float* pZ, int* pNinputt, float* inputt, float* dem,
                        int* pIfl, int* pQtherm, float* pVelocity, float* photar, int* status)
{

// Subroutine to calculate the summed emission from plasma with a
// DEM given by the input array.
// Arguments :
//      itype   I        i: type of plasma emission file
//                            1 = R-S
//                            2 = Mekal
//                            3 = Meka
//                            4 = APEC
//      swtch  I        i: 0==calculate, 1==interpolate, 2=apec interpolate
//      ear     R        i: model energy ranges
//      ne      I        i: number of model energies
//      abun    R        i: abundances
//      dens    R        i: density (cm^-3)
//      z       R        i: redshift
//      ninputt I        i: number of temperatures
//      inputt  R        i: temperatures
//      dem     R        i: emission measures for input temperatures
//      ifl     I        i: dataset number (unused at present)
//      qtherm  R        i: apply thermal broadening (APEC models only)
//      velocityR        i: gaussian velocity broadening (APEC models only)
//      photar  R        r: spectrum
//      status  I        r: 0==OK
// 
// Translated from sumdem.f (CG 10/08)

   int itype = *pItype;
   int swtch = *pSwtch;
   int ne = *pNe;
   float dens = *pDens;
   float z = *pZ;
   int ninputt = *pNinputt;
   int ifl = *pIfl;
   int qtherm = *pQtherm;
   float velocity = *pVelocity;

   const int TOTEL=30, MKEL=14, RSEL=12, APEL=13;
   const char *celt[TOTEL] = {"H", "He", "Li", "Be", "B", "C", "N", "O", "F", 
                          "Ne", "Na", "Mg", "Al", "Si", "P", "S", "Cl", "Ar", 
                          "K", "Ca", "Sc", "Ti", "V", "Cr", "Mn", "Fe", "Co", 
                          "Ni", "Cu", "Zn"};   
   const float andgrv[] = {1.00e+00, 9.77e-02, 1.45e-11, 1.41e-11, 3.98e-10, 
                3.63e-04, 1.12e-04, 8.51e-04, 3.63e-08, 1.23e-04, 
                2.14e-06, 3.80e-05, 2.95e-06, 3.55e-05, 2.82e-07, 
                1.62e-05, 3.16e-07, 3.63e-06, 1.32e-07, 2.29e-06, 
                1.26e-09, 9.77e-08, 1.00e-08, 4.68e-07, 2.45e-07, 
                4.68e-05, 8.32e-08, 1.78e-06, 1.62e-08, 3.98e-08};
   const int rselt[] = {1, 5, 6, 7, 9, 11, 13, 15, 17, 19, 25, 27};
   const int apelt[] = {1, 5, 6, 7, 9, 11, 12, 13, 15, 17, 19, 25, 27};
   const int mkelt[] = {1, 5, 6, 7, 9, 10, 11, 12, 13, 15, 17, 19, 25, 27};
   
   
   *status = 0;

   // Check for invalid switch and itype combinations
   
   if ((swtch == 0 && itype == 1) || (swtch == 0 && itype == 4) ||
          (swtch == 1 && itype == 3))
   {
      char errMsg[] = "***SUMDEM: Invalid switch and itype combination";
      xs_write(errMsg,2);
      return;
   }

   std::ostringstream oss;
   oss << "SUMDEM : itype = " << itype;
   xs_write(const_cast<char *>(oss.str().c_str()), 25);
   
   // set abundances - we use the Anders & Grevesse abundance
   // scale internally so shift from the Solar abundance table in use.
   
   float nabun[TOTEL];
   for (int i=0; i<TOTEL; ++i)
      nabun[i] = 0.0;
      
   if (itype == 1)
   {
      for (int i=0; i<RSEL; ++i)
         nabun[i] = abun[i]*FunctionUtility::getAbundance(string(celt[rselt[i]]))
                        / andgrv[rselt[i]];
   }
   else if (itype == 4)
   {
      for (int i=0; i<APEL; ++i)
         nabun[i] = abun[i]*FunctionUtility::getAbundance(string(celt[apelt[i]]))
                        / andgrv[apelt[i]];
   }
   else
   {
      for (int i=0; i<MKEL; ++i)
         nabun[i] = abun[i]*FunctionUtility::getAbundance(string(celt[mkelt[i]]))
                        / andgrv[mkelt[i]];
   }
   
   // If itype =2 or 3 but switch =2 then we are using one of the mekal
   // models but want to apply the apec code. Need to correct the nabun
   // array for the extra element included in mekal

   if ((itype == 2 || itype == 3) && swtch == 2 )
      for (int i = 5; i<APEL; ++i)
         nabun[i] = nabun[i+1];
      
   // If the calculate option is set...
   if (swtch == 0)
      smdem0_(itype, ear, ne, nabun, dens, z, ninputt, inputt, dem, ifl, 
                photar, *status);

   // Else if the interpolate option is set...
   else if (swtch == 1)
      smdem1_(itype, ear, ne, nabun, dens, z, ninputt, inputt, dem, ifl, 
                photar, *status);

   // Else if the APED interpolate option is set...
   else if (swtch == 2)
      smdem2_(itype, ear, ne, nabun, dens, z, ninputt, inputt, dem, ifl, 
                qtherm, velocity, photar, *status);
   else
   {
      char errMsg[] = "***SUMDEM : Unknown value of switch";
      xs_write(errMsg, 5);
   }
 
   std::vector<double>& tdem = FunctionUtility::tempsDEM();
   std::vector<double>& edem = FunctionUtility::DEM();
   tdem.resize(ninputt);
   edem.resize(ninputt);
   for (int i=0; i<ninputt; ++i)
   {
      tdem[i] = static_cast<double>(inputt[i]);
      edem[i] = static_cast<double>(dem[i]);
   } 
}