**==workhorse.spg  processed by SPAG 4.50J  at 14:50 on 30 Jun 1995
c ***************************************************************************
c  This program runs the Raymond-Smith code to produce tabulated data files
c  for use by XSPEC.
c
c  original version   pat
c  modified 18/7/86   kaa
c  modified 24/2/87   fwjh
c  modified 30/6/95   kaa    removed structures and switched to FITS 
c                            format output 
c
c ****************************************************************************
 
      PROGRAM WORKHORSE

 
c
c	WORKHORSE <filename> <start T> <end T> <nsteps>
c
 
      INCLUDE 'rayspec.inc'
 
      INTEGER NELEMENT, MAXTMP
      PARAMETER (NELEMENT=12, MAXTMP=500)
 
      REAL efilear(0:MXBINS), parval(MAXTMP)
      REAL estep(MXBINS) , emin(MXBINS)
      REAL bin1(MXBINS) , bin2(MXBINS) , enorm(MXBINS)
      REAL frac(NELEMENT)
      REAL addivl(6, NELEMENT), intivl(6)
      REAL tstart , tend , tcon , tt , tstep , rrec , dene , t
 
      INTEGER nsteps(MXBINS)
      INTEGER iparse , nret , iflag , ierr , clen , cpage , iunit
      INTEGER ntsteps , ien , ien1 , iset , lenrec , nrec
      INTEGER idens , icx , i , j , k , l, status
 
      CHARACTER filename*64 , string*255, contxt*72
      CHARACTER*8 elements(NELEMENT)
 
      DATA elements/'He' , 'C ' , 'N ' , 'O ' , 'Ne' , 'Mg' ,
     &              'Si' , 'S ' , 'Ar' , 'Ca' , 'Fe' , 'Ni'/
      DATA dene/0./
      DATA idens , icx/1 , 1/
 
c Get the input parameters for this run
 
      CALL XINIRD('Input <filename> <start T/keV> <end T/keV> <nsteps>:'
     &            ,string,iparse)
      IF ( iparse.GE.LEN(string) ) CALL EXIT(0)
      CALL XGTSTR(string,iparse,1,'filename',1,filename,nret,iflag,-1)
      CALL XTEND(filename,'mod')
      CALL XGTLUN(iunit,ierr)
      CALL XGTRL4(string,iparse,1,'start T',1,tstart,nret,iflag,-1)
      CALL XGTRL4(string,iparse,1,'end T',1,tend,nret,iflag,-1)
      CALL XGTINT(string,iparse,1,'ntsteps',1,ntsteps,nret,iflag,-1)
      CALL XNXTRD('enter <start energy> in eV:',string,iparse,ierr)
      CALL XGTRL4(string,iparse,1,'start energy',1,efilear(0),nret,
     &            iflag,-1)
 
c Get any more sets of energy bins and set up the energy bin array
c efilear is in eV.
 
      iset = 0
      DO WHILE ( .TRUE. )
         CALL XNXTRD('enter <energy step> and <no of steps>:',string,
     &               iparse,ierr)
         iset = iset + 1
         CALL XGTRL4(string,iparse,1,'E step',1,estep(iset),nret,iflag,
     &               -1)
         IF ( nret.LE.0 ) THEN
            iset = iset - 1
            GOTO 100
         ENDIF
         CALL XGTINT(string,iparse,1,'nsteps',1,nsteps(iset),nret,iflag,
     &               -1)
         emin(iset) = efilear(ien)
         DO i = 1 , nsteps(iset)
            ien = ien + 1
            IF ( ien.GT.MXBINS ) GOTO 600
            efilear(ien) = efilear(ien-1) + estep(iset)
            enorm(ien) = 1000.*2./(efilear(ien)+efilear(ien-1))/1.60207
 
c  norm by bin mid energy to change from ergs to photons
 
         ENDDO
      ENDDO

 100  CONTINUE

c Set the initial parameter values

      DO i = 1, NELEMENT
         addivl(1,i) = 1.
         addivl(2,i) = -0.01
         addivl(3,i) = 0.
         addivl(4,i) = 0.
         addivl(5,i) = 1000.
         addivl(6,i) = 1000.
      ENDDO
      intivl(1) = 1.
      intivl(2) = 0.01
      intivl(3) = 0.008
      intivl(4) = 0.008
      intivl(5) = 64.0
      intivl(6) = 64.0

c Store the temperatures to be used in PARVAL.

      tstep = (ALOG10(tend)-ALOG10(tstart))/(ntsteps-1)

      DO i = 1, ntsteps
 
         t = ALOG10(tstart) + (i-1)*tstep
         parval(i) = 10**t

      ENDDO

c Turn the energies into keV.

      DO i = 0, ien
         efilear(i) = efilear(i)/1000.
      ENDDO

c Write the header information, parameter extension, and energy
c extension of the table file.

      CALL wftbmd(filename, ' ', 'ray smith', 'emis meas', 1,
     &            NELEMENT, .TRUE., .TRUE., elements, addivl,
     &            'kT (keV)', intivl, ntsteps, 1, MAXTMP, parval,
     &            ien+1, efilear(0), iunit, ierr, status, contxt)

c Set the initial element abundances to a small number
 
      DO i = 1 , NELEMENT
         frac(i) = 1E-7
      ENDDO
 
c Loop round the required temperatures using log steps. The
c R-S code requires T in Kelvin
 
      tcon = 11.6E6
      tstart = tstart*tcon
      tend = tend*tcon
      tstart = ALOG10(tstart)
      tend = ALOG10(tend)
      tstep = (tend-tstart)/(ntsteps-1)

      DO j = 1 , ntsteps
 
         t = tstart + (j-1)*tstep
         tt = 10**t/tcon
         WRITE (*,'(a,f10.5,a)') '+** calculating model with T/keV= ' , 
     &                           tt

c Write the parameter value column

         CALL ftpcle(iunit, 1, j, 1, 1, tt, status)
         contxt = 'Failed to write PARAMVAL column'
         IF ( status .NE. 0 ) GOTO 999

c Calculate zero-metal continuum
 
         ien1 = 1
         DO k = 1 , iset
            CALL RAYSPEC_M(t, dene, frac, nsteps(k), emin(k), 
     &                     estep(k), idens, icx, bin1(ien1))
            ien1 = ien1 + nsteps(k)
         ENDDO
         DO l = 1, ien
            bin1(l) = bin1(l)*enorm(l)
         ENDDO
 
c Write the zero-metal continuum

         CALL ftpcle(iunit, 2, j, 1, ien, bin1, status)
         contxt = 'Failed to write INTPSPEC column'
         IF ( status .NE. 0 ) GOTO 999
 
c Create the heavy metal spectra
 
         DO i = 1 , NELEMENT
            ien1 = 1
            DO k = 1 , iset
               frac(i) = 1.
               CALL RAYSPEC_M(t, dene, frac, nsteps(k), emin(k), 
     &                        estep(k), idens, icx, bin1(ien1))
               frac(i) = 2.
               CALL RAYSPEC_M(t, dene, frac, nsteps(k), emin(k),
     &                        estep(k), idens, icx, bin2(ien1))
               ien1 = ien1 + nsteps(k)
               frac(i) = 1.E-7
            ENDDO
            DO l = 1, ien
               bin2(l) = enorm(l)*(bin2(l)-bin1(l))
            ENDDO

            CALL ftpcle(iunit, i+2, j, 1, ien, bin2, status)
            WRITE(contxt,'(a,i3)') 'Failed to write column ', i+2
            IF ( status .NE. 0 ) GOTO 999

         ENDDO
 
c End loop over temperatures
 
      ENDDO
 
      CALL ftclos(iunit, status)
      CALL FRELUN(iunit)

      CALL leave()
 
 600  CONTINUE
      WRITE (*,*) 'ERROR: too many energy intervals, maximum:' , MXBINS
      
 999  CONTINUE
      IF ( status .NE. 0 ) THEN
         WRITE(*,'(a)') contxt
         WRITE(*,*) status
         IF ( ierr .NE. 0 ) WRITE(*,*) ierr
      ENDIF
 
      END
