clcwdem.cxx

#include <xsTypes.h>
#include <XSFunctions/Utilities/FunctionUtility.h>
#include <XSUtil/Utils/XSstream.h>
#include <XSstreams.h>

#include <cmath>
#include <iostream>
#include <cstdlib>
#include <sstream>
#include <valarray>
#include <string>

using namespace std;

namespace
{
  // get number of steps in cooling flow
  size_t get_num_steps()
  {
    size_t numsteps = 20;
    const std::string ntempsval =
      FunctionUtility::getModelString("WDEM_NTEMPS");
    if ( ntempsval != FunctionUtility::NOT_A_KEY() )
      {
	// get number of steps if set
	size_t tempstep;
	istringstream buffer(ntempsval);
	
	buffer >> tempstep;
	if( buffer.fail() )
	  {
	    tcerr << "Invalid value in WDEM_NTEMPS\n";
	  }
	else
	  {
	    numsteps = tempstep;
	  }
      }
    return numsteps;
  }

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

  // wrapper to get spectrum of thermal gas
  void calc_sumdem(const RealArray& energy,
		   const std::valarray<float>& abun,
		   float z,
		   const std::valarray<float>& temps,
		   const std::valarray<float>& dem,
		   RealArray& flux)
  {
    // copy energies to floats
    size_t esize = energy.size();
    valarray<float> ear(esize);
    for(size_t i = 0; i != esize; ++i)
      ear[i] = energy[i];

    // output array
    valarray<float> photar(esize-1);
    photar = 0.;

    int itype = 4; // apec
    int swtch = 2; // apec interpolate
    float dens = 1.;
    int numt = temps.size();
    int earsize = ear.size()-1;
    int ifl = 1;
    int qintherm = 0;
    float velocity = 0.;
    int status = 0;
    sumdem_( &itype, &swtch, &(ear[0]), &earsize,
	     const_cast<float*>(&(abun[0])), &dens, &z,
	     &numt, const_cast<float*>(&(temps[0])),
	     const_cast<float*>(&(dem[0])), &ifl, &qintherm, &velocity,
	     &(photar[0]), &status);

    // write output into Real array
    flux.resize(esize-1);
    for(size_t i = 0; i != (esize-1); ++i)
      flux[i] = photar[i]; 
  }

}

// lmod wdem /home/jss/code/xspec/wdem/
extern "C" void clcwdem(const RealArray& energy,
			const RealArray& parameter,
			int spectrum,
			RealArray& flux,
			RealArray& fluxError,
			const string& init)

{
  const double thigh = parameter[0];
  const double tlow = parameter[1];

  if( tlow >= thigh )
    {
      flux.resize( energy.size()-1 );
      flux = 0;
      return;
    }

  const double index = parameter[2];
  const double redshift = parameter[16];

  //printf("%e %e %e %e\n", thigh, tlow, index, redshift);

  // need abundances as floats :-(
  valarray<float> abun(13);
  for(size_t i = 0; i != 13; ++i)
    abun[i] = parameter[i+3];

  const size_t numsteps = get_num_steps();
  valarray<float> arry_dem(numsteps);
  valarray<float> arry_T(numsteps);

  const double logtlow = log(tlow);
  const double logthigh = log(thigh);
  const double logdelta = (logthigh-logtlow) / numsteps;

  // work out emission measures for each temperature
  double total_em = 0.;
  for(size_t i = 0; i != numsteps; ++i)
    {
      const double T = exp( logtlow + logdelta*(i+0.5) );
      if ( T >= 0.0808 )
	{
	  arry_T[i] = T;
	  const double deltaT = exp(logtlow + logdelta*(i+1)) -
	    exp(logtlow + logdelta*i);
	  const double em = pow( (T/thigh), index) * deltaT;

	  if ( tpout.maxChatter() >= 15 )
	    {
	      tcout << "T=" << T << '\n'
		    << " deltaT=" << deltaT << '\n'
		    << " T/Thigh=" << T/thigh << '\n'
		    << " index=" << index << '\n'
		    << " em=" << em << '\n';
	    }

	  //printf(" %e\n", em);
	  arry_dem[i] = em;
	  total_em += em;
	}
      else
	{
	  arry_T[i] = 0.0808;
	  arry_dem[i] = 0.;
	}
    }

  if(total_em > 1e-30)
    {
      // now normalize to give 1. as total
      arry_dem *= (1./total_em);
    }
  else
    {
      // use isothermal when index gets too high
      arry_dem = 0.;
      arry_dem[numsteps-1] = 1.;
    }

    if( tpout.maxChatter() >= 12 )
      {
	tcout << "wdem (" << numsteps << " temperature steps)\n";
	for(size_t i=0; i!=numsteps; ++i)
	  {
	    tcout << "wdem: T=" << arry_T[i] << ", demval=" << arry_dem[i]
		  << '\n';
	  }
	tcout << "\n";
      }

  // have to copy arrays to floats to send to sumdem
  calc_sumdem(energy, abun, redshift, arry_T, arry_dem,
	      flux);
}

model.dat

wdem          17  0.         1.e20           C_clcwdem add  0
highT   keV     4.    0.0808  0.0808 79.9      79.9       0.001
lowT    keV     0.5   0.0808  0.0808 79.9      79.9       0.01
indx     " "    2.    0.      0.   1000.     1000.        0.01
He       " "    1.    0.      0.   1000.     1000.       -0.01
C        " "    1.    0.      0.   1000.     1000.       -0.01
N        " "    1.    0.      0.   1000.     1000.       -0.01
O        " "    1.    0.      0.   1000.     1000.       -0.01
Ne       " "    1.    0.      0.   1000.     1000.       -0.01
Mg       " "    1.    0.      0.   1000.     1000.       -0.01
Al       " "    1.    0.      0.   1000.     1000.       -0.01
Si       " "    1.    0.      0.   1000.     1000.       -0.01
S        " "    1.    0.      0.   1000.     1000.       -0.01
Ar       " "    1.    0.      0.   1000.     1000.       -0.01
Ca       " "    1.    0.      0.   1000.     1000.       -0.01
Fe       " "    1.    0.      0.   1000.     1000.       -0.01
Ni       " "    1.    0.      0.   1000.     1000.       -0.01
redshift " "    0.    0.      0.   10.       10.         -0.01