#include <stlToCArrays.h>
#include <xsFortran.h>
#include <functionMap.h>
#include <XSUtil/Utils/XSutility.h>
#include <xsTypes.h>
#include <cmath>

extern "C" void sumdem_(int* itype, int* swtch, float* ear, int* ne, float* abun, float* dens,
        float* z, int* nTemp, float* temps, float* dem, int* ifl, int* qtherm,
        float* velocity, float* photar, int* status);

void cemVMekal(const RealArray& energyArray, const RealArray& params,
        int spectrumNumber, RealArray& flux, RealArray& fluxErr, 
        const string& initString)
{
   /*  
     c
     c XSPEC model subroutine to calculate Continous Emission
     c Measure of the form :
     c   Q(T) = Norm*(T/Tmax)^alpha  
     c See, for example, Schmitt et al. ApJ 365, 704 (1990), 
     c but note that this program yields Schmitt's alpha - 1.0.
     c
     c This program calls 'vmeka' thus allowing one to vary the elemental
     c abundances.
     c 
     c Parameters:
     c    param(1) = slope of CEM, alpha
     c    param(2) = maximum temperature, tmax
     c    param(3) = nH (cm^-3)  Fixed at 1 for most applications
     c    param(4) = He abundance
     c    param(5) = C   "
     c    param(6) = N   "
     c    param(7) = O   "
     c    param(8) = Ne  "
     c    param(9) = Na  "
     c    param(10)= Mg  "
     c    param(11)= Al  "
     c    param(12)= Si  "
     c    param(13)= S   "
     c    param(14)= Ar  "
     c    param(15)= Ca  "
     c    param(16)= Fe  "
     c    param(17)= Ni  " 
     c    param(18) = redshift used in the Mewe-Kaastra plasma model
     c    param(19) = switch(0=calculate MEKAL model, 1=interpolate MEKAL model)
     c
     c K. P. Singh    April 22, 1994
     c
     c Disclaimer: Any resemblance to a real program is purely
     c             coincidental
     c
       Translated from Fortran cevmkl.f by C. Gordon   Dec. 2007
   */
   
   using namespace XSutility;
   
   const Real k = 8.6171e-8;
   const Real alpha = params[0];
   const Real tMax = params[1];
   const Real max = log10(tMax/k);
   int nt = static_cast<int>((max - 5.5)*10.0 + 1.0);
   
   auto_array_ptr<float> apTarr(new float[nt]);
   auto_array_ptr<float> apDem(new float[nt]);
   float* tarr = apTarr.get();
   float* dem = apDem.get();
   
   /*
      c Integrate contributions in form:
      c    f = (sum over i) Wi * Fi * logdeltT
      c where
      c   Wi= (Ti/Tmax)^alpha and Fi is F(Ti) from meka:
      c
      c it is important to do this in uniform Log(T) steps.
      c
      c I am shuffling the definitions of Fi and photar
      c temporarily for ease.
      c
      c Note: Doing the stepping from logT=log(Tmax) in nt 0.1 steps.
      c       This suggestion from Chris Done to ensure that changing
      c       Tmax by less than 0.1 gives a change in chi-squared.
   */
   
   for (int i=nt; i>=1; --i)
   {
      Real logTemp = max - 0.1*(i-1);
      tarr[i-1] = k*pow(10.0, logTemp);
      dem[i-1] = 0.1*pow(tarr[i-1]/tMax, alpha);      
   }
   int swtch = static_cast<int>(params[18]); 
   
   float *ear=0, *pars=0, *photar=0, *photer=0;
   XSFunctions::stlToFloatArrays<float>(energyArray, params, flux, fluxErr,
        ear, pars, photar, photer);
   auto_array_ptr<float> apEar(ear);
   auto_array_ptr<float> apPars(pars);
   auto_array_ptr<float> apPhotar(photar);
   auto_array_ptr<float> apPhoter(photer);
   
   int itype=2, ne=energyArray.size()-1, ierr=0;   
   int qtherm = 0;
   float velocity = 0.0;
   sumdem_(&itype, &swtch, ear, &ne, &pars[3], &pars[2], &pars[17], &nt, 
                tarr, dem, &spectrumNumber, &qtherm, &velocity, photar, &ierr);
   XSFunctions::floatFluxToStl<float>(photar, photer, ne, false, flux, fluxErr);
   if (ierr != 0)
   {
      char msg[] = "***Error returned from SUMDEM called from cevmkl";
      xs_write(msg, 10);
   }
}