      REAL FUNCTION ELRAT(R,DELTIM,ELRATI)
*
*        Calculates electron/mean temperature ratio for adiabatic process with
*        Coulomb collisions between ions and electrons. 
*        Relative accuracy of the fit - 7x10**-5 (max. dev.).
*        INPUT: R       - total particle density/mean temperature^3/2 (cgs)
*               DELTIM  - elapsed time (s)
*               ELRATI  - initial electron/mean temperature ratio
*
      REAL R, DELTIM, ELRATI
      REAL TOT,EOT,OOS,OON,COULMB,CONST
      REAL RT,X0,SX0,Y0,Y,Z
      PARAMETER ( TOT = 2./3. )
      PARAMETER ( EOT = 8./3. )
      PARAMETER ( OOS = 1./7. )
      PARAMETER ( OON = 1./9. )
      PARAMETER ( COULMB = 30.0 )
      PARAMETER ( CONST = 2.*COULMB/503. )
      IF( ELRATI .EQ. 1.0 )THEN
         RT = ELRATI
         GO TO 999
      ENDIF
      X0 = ELRATI
      SX0 = SQRT(X0)
      IF( X0 .GT. 0.01 )THEN
         Y0 = ALOG( (1.+SX0)/(1.-SX0) ) - TOT*SX0*(X0+3.) 
       ELSE
         Y0 = 2.*SX0*X0**2*( 0.2 + X0*( OOS + X0*OON ) )
      ENDIF    
      Y = CONST*DELTIM*R + Y0
      IF( Y .LT. 3.95 )THEN
         Z = (2.5*Y)**0.4
         RT = Z*(1.0000271+Z*(-0.28649804+Z*(0.58116559E-02+Z*(
     .          -0.22288787E-02+Z*(0.48568998E-02-Z*0.78976963E-03)))))
       ELSE
         RT = ( TANH( 0.5*(Y+EOT) ) )**2
      ENDIF
 999  ELRAT = RT
      RETURN
      END
*
*
*
      subroutine sedov(emz,tshock,taus,beta,tempr,tau,npt)

      real emz, tshock, taus, beta, tempr(*), tau(*)
      integer npt
*
*        Provides temperature and ionization timescale for Sedov models.
*
*        INPUT:  emz    -  fractional emission measure (instead of 
*                          fractional radius)
*                tshock -  mean temperature at the shock (K)
*                taus   -  ionization parameter nt (s/cm^3)
*                beta   -  electron/mean temperature ratio at the shock
*        OUTPUT: tempr  -  temperatures (K)
*                tau    -  ionization timescale (s/cm^3)
*                npt    -  number of points
*

      integer n
      parameter ( n = 1000 )
      real temprp(11), taup(11)
      real em(n),v(n),rksi(n),zeta(n),p(n),rho(n),g(n),t(n)
      real tmin, tmax, pp, zetaiz, delt, rhop, zetap, gp
      real taumax, heabun, xmu, xmue, convn, r, delg, delf, tl1
      real gp1, tl, sum, tiz, tf
      integer lenn, iz, iz1, k, nptp, i, j, l, lun, block, ierr, hdutyp
      character*255 pathd, filenm
      logical qini, qanyf
      real elrat
      character fgmodf*255
      integer lenact
      external lenact, fgmodf, elrat

      save qini, em, v, rksi, zeta, p, rho, g, t

      data qini /.true./
*        Read Sedov solution.
      if( qini )then
         pathd = fgmodf()
         lenn = lenact(pathd)
         filenm = pathd(:lenn) // 'sedov.fits'
         CALL getlun(lun)
         ierr = 0
         CALL ftopen(lun, filenm, 0, block, ierr)
         CALL ftmrhd(lun, 1, hdutyp, ierr)
         CALL ftgcve(lun, 1, 1, 1, n, 0.0, em, qanyf, ierr)
         CALL ftgcve(lun, 2, 1, 1, n, 0.0, v, qanyf, ierr)
         CALL ftgcve(lun, 3, 1, 1, n, 0.0, rksi, qanyf, ierr)
         CALL ftgcve(lun, 4, 1, 1, n, 0.0, zeta, qanyf, ierr)
         CALL ftgcve(lun, 5, 1, 1, n, 0.0, p, qanyf, ierr)
         CALL ftgcve(lun, 6, 1, 1, n, 0.0, rho, qanyf, ierr)
         CALL ftgcve(lun, 7, 1, 1, n, 0.0, g, qanyf, ierr)
         CALL ftgcve(lun, 8, 1, 1, n, 0.0, t, qanyf, ierr)
         CALL ftclos(lun, ierr)
         CALL frelun(lun)

cc         open(unit=8,file=pathd(:lenn) // 'sedov.dat',
cc     *                          status='old',form='formatted')
cc         read(8,900) 
cc     *      (em(i),v(i),rksi(i),zeta(i),p(i),rho(i),g(i),t(i),i=1,n)
cc 900     format(8E12.5)
cc         close(8)
         qini = .false.
      endif
*        
         call locate(em,n,emz,iz)
         iz1 = iz + 1
         pp = (emz-em(iz))/(em(iz1)-em(iz))
         tmin=(1.-pp)*p(iz)/rho(iz) + pp*p(iz1)/rho(iz1)
         zetaiz=(1.-pp)*zeta(iz) +pp*zeta(iz1)
         tmax=1./zetaiz**3
         npt = 11
*         npt=int( alog(tmax/tmin)/alog(1.02) ) + 2
*         npt=int( alog(tmax/tmin)/alog(1.05) ) + 2
         delt=(tmin-tmax)/float(npt-1)
         do 10 i = 1, npt
            tempr(i) = tmax + float(i-1)*delt
            rhop = (tempr(i)/tmax)**1.5
            call locate(rho,n,rhop,j)
            pp = (rhop-rho(j))/(rho(j+1)-rho(j))
            zetap = zeta(j)*(1.-pp) + zeta(j+1)*pp
            gp = g(j)*(1.-pp) + g(j+1)*pp
            tau(i) = max(zetaiz**2.5*gp,0.0)
 10      continue
         taumax = tau(npt)
*      if( beta .lt. 1. )then
*           Assumes cosmic H and He abundances (only trace of heavy
*                                                   elements).
         heabun = 1.0/10**1.01
         xmu = (1.+4.*heabun)/(2.+3*heabun)
         xmue = (1.+2.*heabun)/(1.+4*heabun)
         convn = 1./(xmu*xmue)
         r = 4.*convn*taus*zetaiz**7*(1./zetaiz**2.5-1.)/tshock**1.5
         delg=0.1*tau(npt)/zetaiz**2.5
         k = 1
         tau(k) = 0.0
         tempr(k) = beta*tmax
         gp = float(1)*delg
         call locate(g,n,gp,j)
         pp = (gp-g(j))/(g(j+1)-g(j))
         zetap = zeta(j)*(1.-pp) + zeta(j+1)*pp
         rhop = rho(j)*(1.-pp) + rho(j+1)*pp
         delf = (1./zetap**2.5-1.)/(1./zetaiz**2.5-1.)
         tl1 = elrat(r,delf,beta)*tmax*rhop**(2./3.)
         gp1 = gp
         do 30 l=2,10
            gp = float(l)*delg
            call locate(g,n,gp,j)
            pp = (gp-g(j))/(g(j+1)-g(j))
            zetap = zeta(j)*(1.-pp) + zeta(j+1)*pp
            rhop = rho(j)*(1.-pp) + rho(j+1)*pp
            delf = (1./zetap**2.5-1.)/(1./zetaiz**2.5-1.)
            tl = elrat(r,delf,beta)*tmax*rhop**(2./3.)
            k = k + 1
            tempr(k) = tl1
            tau(k) = zetaiz**2.5*gp1
            tl1 = tl
            gp1 = gp
 30      continue
         k = k + 1
         delf = 1.0
         tempr(k) = elrat(r,delf,beta)*tmin
         tau(k) = taumax
         npt = k
*      endif
*
      go to 19
*      tauh = 0.5*tau(npt)
*      call locate(tau,npt,tauh,j)
*      if( (tauh-tau(j))/(tau(j+1)-tau(j)) .lt. 0.5 )then
*         nh = j
*       else
*         nh = j + 1
*      endif
*      if( nh .eq. 1 .or. nh .eq. npt )then
         sum = 0.0
         do 40 i = 2,npt
            sum = sum + 0.5*(tempr(i-1)+tempr(i))*(tau(i)-tau(i-1))
 40      continue
         tiz = sum/tau(npt)
*         tau(1) = 0.0
*         tau(2) = tau(npt)
*         tempr(1) = 2.*tiz - tempr(npt)
*         tempr(2) = tempr(npt)
**         tempr(1) = tempr(2)
*         npt = 2
         tf = tempr(npt)
*         npt = 5
         nptp = 11
         taup(1) = 0.0
         temprp(1) = tempr(1)
         do 42 i=2,nptp-1
            taup(i) = float(i-1)/float(nptp-1)*taumax
            call locate(tau,npt,taup(i),j)
            pp = (taup(i)-tau(j))/(tau(j+1)-tau(j))
            temprp(i) = tempr(j)*(1.-pp) + tempr(j+1)*pp
 42      continue
         taup(nptp) = taumax
         temprp(nptp) = tempr(npt)
         npt = nptp
*         tm = 0.5*( tf + tiz )
*         ti = 2.*tiz - tm
*         tempr(3) = tf
*         tempr(2) = 2.*tm - tf
*         tempr(1) = 2.*ti - tempr(2)
*         tau(2) = 0.5*taumax
*         tau(3) = taumax
*         tau(1) = 0.0
*         npt = 3
*       else
*         tf = tempr(npt)
*         sum = 0.0
*         do 45 i = nh+1,npt
*            sum = sum + 0.5*(tempr(i-1)+tempr(i))*(tau(i)-tau(i-1))
* 45      continue
*         tm = sum/( tau(npt) - tau(nh) )
*         sum = 0.0
*         do 50 i = 2,nh
*            sum = sum + 0.5*(tempr(i-1)+tempr(i))*(tau(i)-tau(i-1))
* 50      continue
*         ti = sum/tau(nh)
*         tempr(3) = tf
*         tempr(2) = 2.*tm - tf
*         tempr(1) = 2.*ti - tempr(2)
*
*         tauh=tau(nh)
*         tau(2) = tauh
*         tau(3) = taumax
*         tau(1) = 0.0
*         npt = 3
*      endif
 19   continue
*        Normalize temperature and ionization timescale to physical units.
      do 20 i = 1,npt
         tempr(i) = tshock*tempr(i)
         tau(i) = taus*tau(i)
 20   continue
      return
      end
*
*
*
      SUBROUTINE LOCATE( XX, N, X, J )

      REAL XX(*)
      REAL X
      INTEGER N, J

      INTEGER JL, JU, JM

      JL = 0
      JU = N + 1
  10  IF( JU - JL .GT. 1 )THEN
         JM = ( JU + JL )/2
         IF( (XX(N).GT.XX(1)) .EQV. (X.GT.XX(JM)) )THEN
            JL = JM
          ELSE
            JU = JM
         ENDIF
         GO TO 10
      ENDIF
*        Routine altered here.
      J = MAX(JL,1)
      RETURN
      END
*
*
*
*
        subroutine contm(t,ei,nt,ear,ne,z,workphot)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c	This subroutine calculates X-ray continua for a set of
c                                            temperatures.
c
c	Input:  t          -   electron temperatures (K)
c               nt         -   number of temperatures
c               ei         -   ionic concentrations by number 
c                              (but could also be
c                              emission integral of each ion)
c                              Hydrogen and helium are completely ionized.
c                              1:       H^+
c                              2:       He^++
c                              3-9:     C ions (including neutrals)
c                              10-17:   N
c                              18-26:   O
c                              27-37:   Ne
c                              38-50:   Mg
c                              51-65:   Si
c                              66-82:   S
c                              83-103:  Ca
c                              104-130: Fe
c                              131-159: Ni
c               ear         -  energy bins
c               ne          -  number of energy bins
c               z           -  redshift
c       Output: workphot(ibin,ielem)    -  bin continuum fluxes for each 
c                               element ielem (-1 for H^+, 0 for He^++,
c                               1-10 for heavy elements)
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*       t     = temperature (K)
*       ei    = ionization fractions
        real keV
        parameter ( keV = 1.60219e-9 )
        integer ne
        real ear(0:ne)
        integer nion,nionp
        parameter (nion=157,nionp=nion+2)
        integer nt, it
        real ei(nionp,*)
        real t(*), z
c       nfreq = total number of continuum frequency points
c       ifreq = current index of continuum frequency point
c	xfreq = dummy real variable
c	freq(nfreq) = grid of photon energies (keV)c		
        integer nfreq, ifreq
        parameter (nfreq=69)
        real freq(nfreq)
*
*          nelem         -  number of heavy elements
*          ielem         -  curent index of elements (-1 for H^+,
*                           0 for He^++, 1-10 for heavy elements
        integer nelem, ielem
        parameter (nelem=10)
c
c	accumulated spectra
c
c			cont(ifreq,ielem) = continnum only xray spectrum
c					(ergs/sec/keV)
c	cont(ifreq,ielem) = continuum
c
       real cont(nfreq,-1:nelem)

      INTEGER numcmbins
      integer nmax
      parameter ( nmax = 10000 )
      INTEGER zstart(nmax),indexin(nmax)
      REAL wl(nmax),wr(nmax)
      real workphot(nmax,-1:nelem)

c        energies in keV at which continuum is calculated
      data freq /.0544,.069,.086,.1,.122,.123,.166,.167,.217,.218,
     *  .3,.33,.34,.36,.37,.488,.49,.522,.523,.666,.668,.706,.708,
     *  .87,.872,1.,1.156,1.158,1.359,1.365,1.582,1.75,1.77,
     *  1.949,1.97,2.04,2.05,2.39,2.44,2.66,2.68,3.223,3.225,
     *  3.48,3.5,5.128,5.13,5.43,5.5,6.1,6.9,7.7,8.82,8.83,9.27,9.28,
     *  10.27,10.29,10.78,10.8,13.6,17.1,21.5,27.7,34.2,43.0,54.1,
     *  68.1,85.8/

c	   zero out spectra
        do 100 ifreq=1,nfreq
           do 110 ielem = -1, nelem
              cont(ifreq,ielem)=0.0
 110       continue
 100    continue
*          Continuum emission
        do 115 it = 1, nt
           call emisc(t(it),ei(1,it),cont,freq)
 115    continue
        do 120 ifreq=1,nfreq
           do 130 ielem = -1, nelem
              cont(ifreq,ielem)=cont(ifreq,ielem)/(keV*freq(ifreq))
 130       continue
 120    continue
*          Initialize rebinning
        call znibin(nfreq,freq,ne,ear(1),z,
     &            numcmbins,indexin,zstart,wl,wr)
        do 20 ielem = -1,nelem
          call zrebin(cont(1,ielem),workphot(1,ielem),ne,numcmbins,
     &               indexin,wl,wr,zstart,freq)
 20     continue
        return
        end

c ************************************************************************
      SUBROUTINE dosedov(nt, delarr, nionp, ei, ntemp, eisum, t, fout,
     &                   ntempr, tempr, tau, ts, tau0, beta, nzmax, 
     &                   ntmp, ionel, ionstage)

      INTEGER nt, nionp, ntemp, ntempr

      REAL delarr(nt), ei(nionp, nt), eisum(nionp, ntemp)
      REAL tempr(ntempr), tau(ntempr), t(ntemp)
      REAL fout(nionp)

      INTEGER ionel(nionp), ionstage(nionp)

      REAL ts, tau0, beta

      INTEGER nzmax, ntmp


      integer nem, iem, ntt
      parameter ( nem = 400 )
      parameter ( ntt = 5 )


c       the lower and upper limiting values of log temperature,
c               the increment in log temp.
c
      real tll,tlu,dtl
      parameter (tll=4.0e0,tlu=9.0e0,dtl=4.0e-3)

      integer lt, it, i
      integer ion, iion
      real eitot, sumem, sumt, dela
      integer npt
      real emb, eme, emz, delem, deleml
      real tt




*           Ionic concentrations for solar abundances.
*           Zero out ionization fractions and emission measures.

         do 198 it = 1,nt
            delarr(it) = 0.0
            do 199 ion = 1,nionp
               ei(ion,it) = 0.0
 199        continue
 198     continue

         deleml = 3./float(nem/2-1)
         do 200 iem = 1, nem/2-1
            emb = 10**( -4. + float(iem-1)*deleml )
            eme = 10**( -4. + float(iem)*deleml )
            emz = 0.5*( emb + eme )
            delem = eme - emb
            call sedov(emz,ts,tau0,beta,tempr,tau,npt)
            call ionsneqr(tempr,tau,npt,nzmax,nionp,fout,ionel,
     &                    ionstage)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c       calculate temperature index
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
            lt=(alog10(tempr(npt))-tll)/dtl+1.5
            if (lt.lt.1)  lt=1
            if (lt.gt.nt) lt=nt
            it = lt
            do 1200 iion=1,nionp
 1200       ei(iion,it)=fout(iion)*delem+ei(iion,it)
            delarr(it) = delarr(it) + delem
 200     continue
         deleml = (4.+alog10(0.9))/float(nem/2-1)
         do 210 iem = 1, nem/2-1
            emb = 1. - 10**( -4. + float(iem)*deleml )
            eme = 1. - 10**( -4. + float(iem-1)*deleml )
            emz = 0.5*( emb + eme )
            delem = eme - emb
            call sedov(emz,ts,tau0,beta,tempr,tau,npt)
            call ionsneqr(tempr,tau,npt,nzmax,nionp,fout,ionel,
     &                    ionstage)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c       calculate temperature index
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
            lt=(alog10(tempr(npt))-tll)/dtl+1.5
            if (lt.lt.1)  lt=1
            if (lt.gt.nt) lt=nt
            it = lt
            do 1210 iion=1,nionp
 1210       ei(iion,it)=fout(iion)*delem+ei(iion,it)
            delarr(it) = delarr(it) + delem
 210    continue
*          Rebin temperature zones for spectral calculations.
        ntmp = 0
        do 115 i=1,nt-ntt+1,ntt
           eitot = 0.0
           do 116 it=i,i+ntt-1
              if( delarr(it) .gt. 0.0 ) eitot = 1.0
 116       continue
           if( eitot .gt. 0.0 )then
              ntmp = ntmp + 1
              sumem = 0.0
              sumt = 0.0
              do 117 ion=1,nionp
                 eisum(ion,ntmp) = 0.0
 117          continue
              do 118 it=i,i+ntt-1
                 tt=it-1
                 tt=tll+dtl*tt
                 tt=10.0e0**tt
                 dela = delarr(it)
                 sumem = dela + sumem
                 sumt = dela*tt + sumt
                 do 119 ion=1,nionp
                    eisum(ion,ntmp) = ei(ion,it) + eisum(ion,ntmp)
 119             continue
 118          continue
              t(ntmp) = sumt/sumem
           endif
 115    continue

        RETURN
        END

