Table specta to FITS for Xspec?
Randall Smith
rsmith at kracko.harvard.edu
Tue Jul 21 11:53:26 EDT 1998
>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 ?
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
More information about the fitsbits
mailing list