Subject: Table specta to FITS for Xspec?

From: Askar.Ibragimov at ksu.ru (Askar Ibragimov)
Date: Tue, 21 Jul 1998 16:23:16 GMT

Is anyone of this conference is user of XSPEC ?

Now I need to fit spectral data with my own model. I have a group of text files,
produced by my program, with theoretical spectra, printed as "Energy, eV ---
Flux, ergs cm-2 s-1 keV-1" (of course, I can make Flux in photons). I'm beginner
in XSPEC and using FITS, can anyone describe me what (as I understood,
procedures on Fortran found in XSPEC distribution) should I do to make file with
compilation of my models, which can be read from XSPEC ?

Askar

---

From: Randall Smith <rsmith at kracko.harvard.edu>
Date: 21 Jul 1998 11:53:26 -0400

>Now I need to fit spectral data with my own model. I have a group of
>text files, produced by my program, with theoretical spectra, printed
>as "Energy, eV --- Flux, ergs cm-2 s-1 keV-1" (of course, I can make
>Flux in photons). I'm beginner in XSPEC and using FITS, can anyone
>describe me what (as I understood, procedures on Fortran found in
>XSPEC distribution) should I do to make file with compilation of my
>models, which can be read from XSPEC ? 

You're in luck--I've dealt with this problem both in C and Fortran.
One solution I can send you is a C program that will read a single
model spectrum (E(i), F(i)) and output it as an XSPEC table model.
This program won't allow you to _fit_ to the model, as it only works
with a single spectrum.  If you would like a copy, please email me.
It's not long.

I have also worked out how to compile a collection of models into a
table model that can then be used to fit data using XSPEC.  (In my
case, I created models of SNR emission as a function of time, and then
fit to the SNR age.)  To do this, you need to use the wftbmd.f routine
that is in the XSPEC distribution.  (If you'd prefer to use C, I've got
a version of this routine that has been converted to C.)  You also
need to have either the fitsio or the cfitsio library compiled.  Then
you can use the following subroutine (included below) to output the
file.  This subroutine requires the name of the output file
(filename), the name of the model that will be used inside of XSPEC
(modelname), the number of different parameters in the model (npar),
and the energy binning.  I use a constant energy spacing, and so just
have it give the minimum energy (binmin) and the size of the bins
(binsyz), but obviously that's easily changed.  Two common blocks are
used, one which contains all the parameters that describe the model
(xspecparms) and one which contains the spectral data itself
(xspecvals).  

Please let me know if you have any questions about how to make this
work.  Good luck with your models!

Yours,
Randall Smith
rsmith at cfa.harvard.edu

------------------------------------------------------------------------

      subroutine output_v9_xspec(filename,modelname,npar,binmin,binsyz)
      
      implicit none
c
c
c   Author: Randall Smith, Smithsonian Astrophysical Observatory
c   Email : rsmith at cfa.harvard.edu
c   Date Released: July 21, 1998
c 
c

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     MAXNPVALS = Maximum # of parameter values for a given parameter.c
c     MAXPARMS = Maximum # of different parameters; XSPEC sets limit. c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      integer MAXNPVALS,MAXPARMS,MAXBINS
      parameter(MAXNPVALS=20,MAXPARMS=2,MAXBINS=3000)
      
      character*12 parmname(MAXPARMS)
      real*4 initval(MAXPARMS), minval(MAXPARMS), lowval(MAXPARMS)
      real*4 highval(MAXPARMS),maxval(MAXPARMS),deltaval(MAXPARMS)
      integer*2 numpvals(MAXPARMS), interp(MAXPARMS)
      real*4 parvalues(MAXPARMS,MAXNPVALS)

      common/xspecparms/parmname, initval, minval, lowval, highval,
     $     maxval,deltaval, numpvals, interp, parvalues
      
      real channels(MAXNPVALS**MAXPARMS, MAXBINS)  ! Here is the data
      integer nbin
      common/xspecvals/channels,nbin  ! Note: huge...Beware...
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      integer MAXNADDPM
      parameter(MAXNADDPM=1)    ! max addl parameters 
      integer*2 NPERBIN
      parameter(NPERBIN=4) ! defined by XSPEC

      character*128 filename,contxt
      character*12 modelname, modunt, infile
      character*12 addnam(MAXNADDPM), intnam(MAXPARMS)
      integer npar, j
      real binmin, binsyz, energy(MAXBINS), data(MAXBINS)
      
      integer nintpm, naddpm, maxtab, nenerg
      integer intntb(MAXPARMS), intmth(MAXPARMS)
      real addivl(6, MAXNADDPM), intivl(6, MAXPARMS)
      real inttab(MAXNPVALS, 1)
      real paramdata(MAXPARMS)
      logical qrdshf, qaddtv
c     output parameters 
      integer ounit, ierr, status
      integer nelements
      integer row, i
      
      nintpm = npar     ! only one parameter, time, currently
      naddpm = 0        ! no additional parameters
      modunt = 'Years'  ! Model units 
      qrdshf = .FALSE.  ! No redshifting
      qaddtv = .TRUE.   ! Additive model
      addnam(1) = ''    ! No additional parameter names
      intnam(1) = 'Time' ! Time is the unit.
      infile = 'spectrum'
      nenerg = nbin + 1
      do i=1,nenerg
         energy(i) = (binmin + i*binsyz)/1000.  ! Energy bins
      end do

      maxtab = MAXNPVALS

      do i=1,MAXPARMS
         intmth(i) = interp(i) ! 1 = log, 0 = linear interpolation
         intivl(1,i) = initval(i)  ! initial value to use
         intivl(2,i) = deltaval(i) ! inital delta step to take when fitting
         intivl(3,i) = minval(i)   ! hard minimum (can't go under this)
         intivl(4,i) = lowval(i)   ! soft minimum (mininum fitted #)
         intivl(5,i) = maxval(i)   ! hard maximum (can't go over this)
         intivl(6,i) = highval(i)  ! soft maximum (maximum fitted #)
         intntb(i) = numpvals(i)   ! # of diff. parm. vals (spectra) avail.

         do j=1,numpvals(i)               ! Set the actual values of the 
            inttab(j,i) = parvalues(i,j)  ! parameters...ie, t=10^5 yrs,
         end do                           ! n = 10^11 cm^-3, T=10^6 K etc.
      end do

      ounit = 89        ! output unit to use.
c
c  This call writes all the necessary header information.
c
      call wftbmd(filename, infile, modelname, modunt,
     $     nintpm, naddpm, qrdshf, qaddtv, addnam, addivl,
     $     parmname, intivl, intntb, intmth, maxtab,
     $     inttab, nenerg, energy, ounit, ierr, status, contxt)

c
c  Now write the model data
c
      do row=1,numpvals(1)
         nelements = 1
         paramdata(1) = parvalues(1,row)
         call ftpcle(ounit, 1, row, 1, nelements, paramdata, status)
         
         nelements = nbin
         do i=1,nbin 
            data(i) = channels(row,i)
         end do
         call ftpcle(ounit, 2, row, 1, nelements, data, status)
      end do
c
c  Now close the file
c
      call ftclos(ounit, status)
      
      return
      end