c *****************************************
c Compton reflection from relativistic disk 
      subroutine relrefl(xr,refl,xm,gamma,sprf,  
     + par_refl,nrefl,jref,ioniz) 
c *****************************************
      implicit none
      integer ifl
      integer nrefl,jref,ioniz 
      real xr(jref),refl(jref) 
      real par_refl(nrefl)
      real gamma,ab_,ab0,temp,xi,xm,ec
c parameters number for the model used
      integer ndisk,nem
      PARAMETER(nem=5000)

      parameter(ndisk=6)
      integer nx
      real xmin,xmax
c internal integration parameters
      parameter(xmin=2.e-5,xmax=4.,nx=600)
      integer n_disk
      real disk_xmin, disk_xmax
      parameter(disk_xmin = 0.7, disk_xmax = 1.5, n_disk = 400)
     
      logical firstflag,fl1
      integer i,j,k,jlo,jhi
      real x(nx),x_disk(0:n_disk),dx
      real phot_ref(nx),phot_disk(n_disk),phot_rs(nx),spref(nx) 
      real logxmin, logxmax, logdmin, logdmax
      real pp,r,xlo,xhi,disx
      real in_disk(ndisk)
      real spint,eps,dln  
      external sprf 
      data firstflag/.true./,fl1/.true./ 
      data eps/1e-5/,ec/0./,ifl/0/   
      data in_disk/ndisk*9999./
      save firstflag,x,phot_ref,phot_disk,in_disk,
     >     phot_rs,x_disk,spref,dx

c     Driver for angle-dependent reflection from an exponentially-cutoff
c     power law and ionized medium.
c     It needs to be compiled together with xspexrav.f (neutral reflection).
c
c     See Magdziarz & Zdziarski, 1995, MNRAS.
c     See Zdziarski et al., 1995, ApJL January for description of 
c     calculation of ionization (based on Done et al. 1992).
c     The abundances are of Morrison & McCammon.
c
c     The output spectrum is the sum of the cutoff power law and the 
c     reflection component. 
c     The reflection component alone can be obtained 
c     for scale (see below) = rel_refl < 0. Then the actual 
c     reflection normalization is |scale|. Note that you need to 
c     change then the limits of rel_refl. The range of rel_refl in that case 
c     should exclude zero (as then the direct component appears).
c     
c     If E_c=0, there is no cutoff in the power law. 
c
c     This version allows for changes of the vector 'ear' between subsequent
c       calls.
c     and then calculates the broanening in a schwartzchild metric
c
c Fe abudance
      ab0=par_refl(1)
c  ab_       - abundance of elements heavier than He relative to
c              the solar abundances
      ab_=par_refl(2)
c  xi        - disk ionization parameter = L/nR^2
      xi=par_refl(3)
c  temp      - disk temperature (for ionization balance, Tmax=10^6 K)
      temp=par_refl(4)

      dln=alog(10.)

      if(firstflag)then
      logxmin = log10(xmin)
      logxmax = log10(xmax) 
      logdmin = log10(disk_xmin)
      logdmax = log10(disk_xmax) 
         dx=log10(xmax/xmin)/(nx-1)
         do i=1,nx
            x(i)=xmin*10**((i-1)*dx)
         enddo
         disx = log10(disk_xmax / disk_xmin) / n_disk
         x_disk(0) = disk_xmin
         do i = 1, n_disk
           x_disk(i) = disk_xmin * 10**(i * disx)
         enddo
      firstflag=.false. 
      dx=dx*dln
      endif

c      if(ioniz.eq.1) 
      call greeion(x,spref,xm,ab_,ab0,sprf,gamma,ec,temp,xi,nx)
c      if(ioniz.eq.0) 
c      call greeneu(x,spref,xm,ab_,ab0,sprf,gamma,ec,nx)
c spref= x Fx
         do i=1,nx
         phot_ref(i) = spref(i)/x(i)**2  
         enddo 
       

c beta 
      if(in_disk(2).ne.par_refl(5))then
         in_disk(2)=par_refl(5)
         fl1=.true.
      endif
c Rin 
      if(in_disk(3).ne.par_refl(6))then
         in_disk(3)=par_refl(6)
         fl1=.true.
      endif
c Rout
      if(in_disk(4).ne.par_refl(7))then
         in_disk(4)=par_refl(7)
         fl1=.true.
      endif
c inclination 
      pp=180./3.1415*acos(xm) 
      if(in_disk(5).ne.pp)then
         in_disk(5)=pp 
         fl1=.true.
      endif
c redshift (we are in the source frame) 
      if(in_disk(6).ne.0.)then
         in_disk(6)=0. 
         fl1=.true.
      endif
      
      if(fl1) then
        in_disk(1) = 1.0
        call relline(x_disk, n_disk, in_disk, ifl, phot_disk)
         do i=1, n_disk
           if(phot_disk(i) .ge. 1.0e-8) then
             phot_disk(i) = phot_disk(i) / (x_disk(i) - x_disk(i-1))
           else
             phot_disk(i) = 0.0
           endif
         enddo
      endif

     
        do i=1,nx
           phot_rs(i)=1e-20 
        enddo
        do i = 1, nx
c f'(x)=int dx' f(x') line(x/x')
c x(i)=x 
          xlo = disk_xmin * x(i)
          xhi = disk_xmax * x(i)
          jlo = int(nx * (log10(xlo) - logxmin) / (logxmax-logxmin))
          jhi = int(nx * (log10(xhi) - logxmin) / (logxmax-logxmin))
          jlo = max(jlo, 1)
          jhi = min(jhi, nx)
          do j = jlo, jhi
c x(j)=x'
            r = x(i) / x(j)
            k = int(n_disk * (log10(r) - logdmin) / (logdmax-logdmin))
            if((k.ge.1) .and. (k.le.n_disk)) then
            phot_rs(i) = phot_rs(i) + phot_ref(j) * phot_disk(k)
            endif
          enddo
        enddo
      do i=1,nx
      phot_rs(i)=phot_rs(i)*dx 
      enddo

      fl1=.false.

      do i=1,jref
         pp=spint(xr(i),x,nx,phot_rs)
         refl(i)=pp*xr(i)**2 
      enddo

      return
      end

c ------------------------------------------------------------------- c
c interpolation 
      real function spint(xx,xt,jt,spt)
 
      implicit   none
      integer    jt
      real       xx,xt(jt),spt(jt)
 
      integer    il,ih
      real       sp
      save       ih
      data       ih/2/

      if((xx.gt.xt(1)).and.(xx.lt.xt(jt)))then
         if(xx.lt.xt(ih))then
            ih=2
         endif
 10      if(xx.gt.xt(ih))then
            ih=ih+1
            goto 10
         else
            il=ih-1
c            sp=spt(il)+
c     &         (spt(ih)-spt(il))*
c     &         (xx-xt(il))/(xt(ih)-xt(il))
        sp=spt(il)*exp(log(spt(ih)/spt(il))*
     *     log(xx/xt(il))/log(xt(ih)/xt(il))) 
         endif
      else
         sp=1e-20 
      endif
 
      spint=sp
 
      return
      end


      SUBROUTINE relline(earp,ne,param,ifl,photar)

      implicit none
      INTEGER ne,ifl
      REAL earp(0:ne),param(*),photar(ne)
      real ear(0:10000)

c  model to calculate the line shape for a rotating accretion
c  disk. does not include GR effects. note that if param(2) is
c  set to 10 then do the special case of an accretion disk
c  emissivity.
c  modified 11/30/95 by P.Zycki to include object's redshift
c  corrected 20/02/96 by P.Magdziarz 

c  parameters :
c	1        line energy
c       2        power law index for emissivity (10 for disk)
c       3        inner radius (GM/c**2)
c       4        outer radius (GM/c**2)
c       5        inclination  (degrees)
c       6        redshift

      INTEGER nr
      PARAMETER (nr=100)

      INTEGER n, k, ie
      REAL ri, ro, dlogr, dera, alp,
     &     tgcsi2, cosdel, fra, ra, dra,
     &     rd(nr), enel, spm, rilog10, zfac  
      REAL beta, betal, betah, zpo, zpol, zpoh
      REAL sincl, cincl, rafact, flux, total
      REAL enobl, enobh

      DATA dera /57.295779/

c convert input parameters - inclination set to radians.

      enel = param(1)
      alp = param(2)
      ri = param(3)
      rilog10 = alog10(ri)
      ro = param(4)
      zfac = 1+param(6)

      sincl = sin(param(5)/dera)
      cincl = cos(param(5)/dera)

      do n=1,ne
         ear(n) = zfac*earp(n)
      end do

c set spectrum to 0

      DO n = 1, ne
         photar(n) = 0.
      ENDDO

c trap case where inner radius is greater than outer

cjp
c      IF ( ri .GE. ro ) THEN
c         CALL xerror(
c     &     'Inner radius > outer radius  -  model is invalid', 5)
c         RETURN
c      ENDIF

c trap case of inner radius being less than the last stable orbit.

c      IF  (ri .LT. 6.) THEN
c         CALL xerror('Inner radius < 6 R_g  -  model is invalid', 5)
c         RETURN
c      ENDIF

c calculate log step in radius

      dlogr = (alog10(ro)-rilog10)/float(nr-1)

c calculate radii

      rd(1) = ri
      DO n = 2, nr
         rd(n) = 10.**(rilog10+float(n-1)*dlogr)
      ENDDO

c big loop for radii. note that the radius cannot reach 3 else
c the metric goes singular

      DO n = 1, nr - 1

         ra = (rd(n)+rd(n+1))/2.

         dra = rd(n+1) - rd(n)

         rafact = sqrt(1.-3./ra)

c if power-law index is less than ten use to calculate emissivity

         IF (alp.LT.9.9) THEN

            fra = ra**(alp)

c else use the accretion disk emissivity law.

         ELSE

            fra = (1.-sqrt(6./ra))/ra**3

         ENDIF

c loop over azimuthal angles

         DO k = 1, 179, 2

            beta = float(k)/dera
            betal = float(k-1)/dera
            betah = float(k+1)/dera

c calculate mean redshift (1+z = zpo) for the bin

            tgcsi2 = (sincl*sin(beta))
     &               **2/(1.-(sincl*sin(beta))**2)
            cosdel = sincl*cos(beta)
     &               /sqrt((sincl*cos(beta))**2+cincl**2)
            zpo = (1.+cosdel/sqrt(ra*(1.+tgcsi2)-2.))/rafact

c and the low and high redshifts for the bin

            tgcsi2 = (sincl*sin(betal))
     &               **2/(1.-(sincl*sin(betal))**2)
            cosdel = sincl*cos(betal)
     &               /sqrt((sincl*cos(betal))**2+cincl**2)
            zpol = (1.+cosdel/sqrt(ra*(1.+tgcsi2)-2.))/rafact

            tgcsi2 = (sincl*sin(betah))
     &               **2/(1.-(sincl*sin(betah))**2)
            cosdel = sincl*cos(betah)
     &               /sqrt((sincl*cos(betah))**2+cincl**2)
            zpoh = (1.+cosdel/sqrt(ra*(1.+tgcsi2)-2.))/rafact


c  enobl and enobh are the lower and upper observed energy from
c  this azimuthal and radial bin

            enobl = min(enel/zpol, enel/zpoh)
            enobh = max(enel/zpol, enel/zpoh)
            total = enobh - enobl

c  calculate total emission from this bin

            flux = fra*ra*dra*2./dera/zpo**3

c  find fractions of emission from this bin to place in each energy
c  range.

            IF (enobh .LT. ear(0) .OR. enobl .GT. ear(ne) ) GOTO 10

            DO ie = 1, ne
               IF ( ear(ie) .GE. enobh ) THEN
                  IF ( ear(ie-1) .LE. enobl ) THEN
                     photar(ie) = photar(ie) + flux
                     GOTO 10
                  ELSEIF ( ear(ie-1) .GE. enobh ) THEN
                     GOTO 10
                  ELSE
                     IF ( total .GT. 0. ) THEN
                        photar(ie) = photar(ie) 
     &                     + flux*(enobh-ear(ie-1))/total
                     ENDIF
                  ENDIF
               ELSEIF ( ear(ie) .GE. enobl ) THEN
                  IF ( ear(ie-1) .GE. enobl ) THEN
                     IF ( total .GT. 0. ) THEN
                        photar(ie) = photar(ie) 
     &                     + flux*(ear(ie)-ear(ie-1))/total
                     ENDIF
                  ELSE
                     IF ( total .GT. 0. ) THEN
                        photar(ie) = photar(ie) 
     &                     + flux*(ear(ie)-enobl)/total
                     ENDIF
                  ENDIF
               ENDIF
            ENDDO

 10         CONTINUE

         ENDDO
      ENDDO

c normalise values to total

      spm = 0.
      DO n = 1, ne
         spm = spm + photar(n)
      ENDDO

c write values

      IF ( spm .NE. 0 ) THEN
         DO n = 1, ne
            photar(n) = photar(n)/spm
         ENDDO
      ENDIF

      RETURN
      END
c ------------------------------------------------------------------ c
c  Green's functions for Compton reflection with asymtotic
c  absorption included; output is (dy+y0)*G(mu,dy,y0)*W(mu,dy,y0),
c  for .0316<y0<100.; in two parts according to dy=2.
c  grxyh1 - subroutine for dy>2
c  grxyl1 - subroutine for pmx(26)<dy<2
c  pm1   - vector of angle-dependent coefficients,
c          set-up by pm1yn subroutine
c  pmx   - vector of angle-dependent coefficients,
c          set-up by pmxyn subroutine
c  dy    - wavelength shift dy=1/x-1/x0
c  y0    - incident photon wavelength y0=1/x0
c  ap    - absorption parameters
 
      real function    grxyh1(pm1,pmx,dy,y0,ap)
      implicit         none
      real             pm1(19),pmx(26),dy,y0,ap(3)
      real             pyx(10),p1,p4,gr,aic1
 
      p1=log(y0)
      p4=1/y0
      pyx(5)=.230+.033*p1+.763*y0**(-.507)+.007*y0**(-1.553)
      pyx(6)=75.180*p4
      pyx(7)=-.100+.115*p1+.100*y0**1.623
      if(y0.lt.1.)then
         pyx(8)=pmx(19)+pmx(20)*p1
      else
         pyx(8)=pmx(21)+pmx(22)*p1+pmx(23)*(pmx(24)+p4)**pmx(25)
      endif
      p1=1/(y0+dy)
      gr=pyx(5)*(1+pyx(6)*p1)**(pyx(7)*(1+pyx(8)*p1))
      p1=1+dy
      gr=gr*pm1(18)*dy**(-1.5)*p1*
     &   (1+pm1(14)/p1)**(pm1(15)+pm1(16)*p1**pm1(17))
      gr=gr*exp(aic1(y0,(dy-.15),ap(1)))
 
      grxyh1=gr
 
      return
      end
 
      real function    grxyl1(pm1,pmx,dy,y0,ap)
      implicit         none
      real             pm1(19),pmx(26),dy,y0,ap(3)
      real             pyx(10),p1,p2,p3,p4,gr,aic1,xla0,xla
 
      p1=log(min(y0,10.))
      p2=p1*p1
      p3=p2*p1
      pyx(1)=p1*(pmx( 1)+pmx( 2)*p1+pmx( 3)*p2+pmx( 4)*p3)
      pyx(2)=p1*(pmx( 5)+pmx( 6)*p1+pmx( 7)*p2+pmx( 8)*p3)
      pyx(3)=p1*(pmx( 9)+pmx(10)*p1+pmx(11)*p2+pmx(12)*p3)
      pyx(4)=p1*(pmx(13)+pmx(14)*p1+pmx(15)*p2+pmx(16)*p3
     &                                             +pmx(17)*p3*p1)
      pyx(9)=1+pmx(18)*(.326*((y0*.1)**1.692-1)+p1-2.303)
      if(y0.lt.10.)then
         p1=y0+2
         p1=log(p1/(y0+dy))/log(p1/y0)
         p2=p1*p1
         p3=p2*p1
         gr=exp(pyx(1)+pyx(2)*p1+pyx(3)*p2+pyx(4)*p3)
      else
         p4=1/(10+dy)
         p1=log(12*p4)*5.4848
         p2=p1*p1
         p3=p2*p1
         p4=pyx(9)*(y0+dy)*p4
         gr=p4*exp(pyx(1)+pyx(2)*p1+pyx(3)*p2+pyx(4)*p3)
      endif
      p1=pm1(1)+pm1(2)*dy
      p2=(pm1(3)+pm1(4)*max(0.,(pm1(5)-dy))**pm1(6)+
     &   pm1(7)*(pm1(8)-dy)**pm1(9))**pm1(10)/
     &   (1+exp(pm1(11)*(pm1(12)-dy**pm1(13))))
      gr=gr*min(p1,p2)
      xla0=1/(1+ap(1)*y0**3)
      xla=1/(1+ap(1)*(y0+dy)**3)
      gr=gr*xla0*(ap(2)+(1-ap(2))*(1+(exp(ap(3)*dy)-2)*sqrt(1-xla))
     &     *exp(aic1(y0,dy,ap(1))))
 
      grxyl1=gr
 
      return
      end
 
      real function    aic1(y0,dy,apc)
      implicit         none
      real             y0,dy,apc
      real             yp,a0,a1,c0,c1,c2,d0l,d0h,d1l,d1h
      real             aih,ail
 
      yp=y0+dy
      a0=apc
      a1=apc**.3333
      c0=.5/a1
      c1=-3.4641*c0
      c2= 0.5774
      d0l=a0*y0**3
      d0h=a0*yp**3
      d1l=a1*y0
      d1h=a1*yp
      ail=y0*(3-log(1+d0l))+c0*log((1-d1l+d1l**2)/(1+d1l)**2)
     &                     +c1*atan(c2*(2*d1l-1))
      aih=yp*(3-log(1+d0h))+c0*log((1-d1h+d1h**2)/(1+d1h)**2)
     &                     +c1*atan(c2*(2*d1h-1))
 
      aic1=aih-ail
 
      return
      end
 
c ------------------------------------------------------------------ c
c  Angle dependent coefficients for Green's function for Compton
c  reflection;
c  xm  - cosine of the angle between normal to the slab and the
c        observation direction
c  pmx - output coefficients vector
 
      subroutine       pmxyn(xm,pmx)
      implicit         none
      real             xm,pmx(26)
      real             p2,p3,p4
 
              p2=xm*xm
              p3=p2*xm
              p4=p3*xm
      pmx( 1)= .8827-.5885*xm+.7988*p2-.4117*p3
      pmx( 2)= .0517+.1076*xm-.1691*p2+.0678*p3
      pmx( 3)= .0014+.0043*xm
      pmx( 4)=-.0003-.0027*xm+.0021*p2
      pmx( 5)=-1.3259+1.6214*xm-3.7007*p2+2.1722*p3
     &        +exp(20.3*(xm-1.063))
      pmx( 6)=  .0790- .4029*xm+ .6619*p2+ .2210*p3
     &        -exp(17.7*(xm-1.041))
      pmx( 7)=  .0751- .1848*xm+ .4068*p2- .4126*p3
     &        +exp( 9.8*(xm-1.185))
      pmx( 8)=- .0020- .0394*xm+ .1004*p2- .0597*p3
      pmx( 9)= 3.2953-3.6996*xm+7.9837*p2-4.6855*p3
      pmx(10)=- .5278+ .9857*xm-1.7454*p2- .3464*p3
     &        +exp(27.6*(xm-0.959))
      pmx(11)=- .1919+ .5798*xm-1.2879*p2+1.4885*p3
     &        -exp(18.3*(xm-0.986))
      pmx(12)=  .0200+ .0832*xm- .0333*p2- .2370*p3
     &        +exp(16.6*(xm-1.086))
      pmx(13)=-2.2779+2.4031*xm-4.0733*p2+1.9499*p3
     &        -exp(19.6*(xm-0.968))
      pmx(14)=  .4790-1.0166*xm+3.1727*p2-4.0108*p3+3.0545*p4
     &        -exp(30.4*(xm-0.957))
      pmx(15)=  .1122- .3580*xm+ .4985*p2- .3750*p3- .5349*p4
     &        +exp(31.2*(xm-0.972))
      pmx(16)=- .0160- .0471*xm- .0536*p2+ .2006*p3+ .2929*p4
     &        -exp(30.6*(xm-1.009))
      pmx(17)=  .0005+ .0021*xm- .0447*p2+ .1749*p3- .2303*p4
     &        +exp(16.9*(xm-1.130))
      pmx(18)=.006+.089*xm-.102*p2+.056*p3
      pmx(19)=.618-.830*xm
      pmx(20)=.128-.132*xm
      pmx(21)=.632-.875*xm
      pmx(22)=-.672-.286*xm+.717*p2-.481*p3
      pmx(23)=.0126-.0160*xm+.0077*p2
      pmx(24)=.0111+.0030*xm-.0014*p2
      pmx(25)=-2.437-.328*xm-.260*p2+.279*p3
      pmx(26)=1-sqrt(1-(xm-.05)**2)
 
      return
      end
 
c ------------------------------------------------------------------ c
c  Angle dependent coefficients for Green's function for Compton
c  reflection of incident photon at energy x0=1;
c  xm  - cosine of the angle between normal to the slab and the
c        observation direction
c  pm1 - output coefficients vector
 
      subroutine       pm1yn(xm,pm1)
      implicit         none
      real             xm,pm1(19)
      real             p2,p3,p4
 
              p2=.911-.549*(1-xm)**1.471
              p3=.254-.041*xm**(-3.798)
      pm1( 1)=p2-2*p3
      pm1( 2)=p3
      pm1( 3)=.161+.439*xm**.189-.791*xm**7.789
      pm1( 4)=max(0.,(-.871+1.740*xm**.176-1.088*xm**.907))
      pm1( 5)=.934+.054*xm**(-.666)
              p4=min(.524,xm)
      pm1( 6)=4.647-11.574*p4+11.046*p4*p4
      pm1( 7)=.012+.199*xm**.939+.861*xm**8.999
      pm1( 8)=2.150
              p2=xm*xm
              p3=p2*xm
      pm1( 9)=-1.281+2.374*xm-4.332*p2+2.630*p3
      pm1(10)=2.882+.035*xm**(-1.065)
      pm1(11)=30.28+69.29*xm
      pm1(12)=1.037*(1-sqrt(1-.874*p2))**.123
      pm1(13)=.123
              p4=p3*xm
      pm1(14)=56.50-27.54*xm+671.8*p2-1245.*p3+708.4*p4
      pm1(15)=-.897-.826*xm+2.273*p2-.984*p3
      pm1(16)=1.240-1.297*xm
      pm1(17)=-.490+1.181*xm-1.038*p2
      pm1(18)=1.858*(xm**2.033+0.745*xm**1.065)
      pm1(19)=1-sqrt(1-(xm-.05)**2)
 
      return
      end
 
c ------------------------------------------------------------------- c
 
      subroutine greeion(x,spref,xm,ab_,ab0,spp,gamma,ec,temp,xi,jmax)
c     calculates the reflected spectrum
c     ab0 - iron log10 abundance relative to hydrogen
      implicit none 
      integer nem
      PARAMETER(nem=5000)
      REAL ab0,ab_,ddy,dy,dym,Ec,fjc,Gamma,GNRJ1,GRXYH1,
     &     GRXYL1,SIGMAION,SPP,sr,srp,Temp,Xi,xjd,xjl,
     &     xkabs,Xm
      REAL xmin,xmo,xnor,xrefmax,xtransh,xtransl,y,y0,ymin 
      INTEGER i,j,Jmax,jmin,jrefmax,jtransh,jtransl,nc,ncp 
      real x(nem),spref(nem)
      real pm1(19),pmx(26),ap(3)
      external spp
 
c  precision factor for Green's function integration
      ncp=100
c  ranges for different methods
      xmin=2.e-4
      xtransl=0.01957
      xtransh=0.02348
      ymin=.03333
c  angle dependent parameters
      xmo=max(.05,min(.95,xm))
      call pm1yn(xmo,pm1)
      call pmxyn(xmo,pmx)
      ap(2)=0.802-1.019*xmo+2.528*xmo**2-3.198*xmo**3+
     &      1.457*xmo**4+(.030/xmo)**4
      ap(3)=0.381
      xrefmax=1/(ymin+pmx(26))
 
      jmin=0
      do 50 j=1,jmax
         if(x(j).gt.xmin) goto 51
         jmin=j
 50   continue
 51   jtransl=0
      do 52 j=max(1,jmin),jmax
         if(x(j).gt.xtransl) goto 53
         jtransl=j
 52   continue
 53   jtransh=0
      do 54 j=max(1,jtransl),jmax
         if(x(j).gt.xtransh) goto 55
         jtransh=j
 54   continue
 55   jrefmax=0
      do 56 j=max(1,jtransh),jmax
         if(x(j).gt.xrefmax) goto 57
         jrefmax=j
 56   continue
 57   continue

      do 10 j=1,jmin
         spref(j)=0
 10   continue
c-----------------------
      do 20 j=jmin+1,jtransl
         xkabs=sigmaion(j,x,jtransh,gamma,temp,xi,ab_,ab0,xnor)
         spref(j)=spp(1/x(j),gamma,ec)*gnrj1(xkabs,xmo)
 20   continue
c-----------------------
      if(jmin.ge.jtransl) then
         xkabs=sigmaion(1,x,max(1,jtransh),gamma,temp,xi,ab_,ab0,xnor)
      endif
      ap(1)=xnor/1.21
      xjl=alog(xtransl)
      xjd=alog(xtransh)-xjl
      do 60 j=jtransl+1,jtransh
         fjc=.5*sin(3.14159*((alog(x(j))-xjl)/xjd-.5))
         xkabs=sigmaion(j,x,jtransh,gamma,temp,xi,ab_,ab0,xnor)
         y=1/x(j)
         spref(j)=spp(y,gamma,ec)*gnrj1(xkabs,xmo)
         dym=y-ymin
         dy=min(2,dym)
         ddy=(dy-pmx(26))/(ncp+1)
         y0=y-dy
         sr=spp(y0,gamma,ec)*grxyl1(pm1,pmx,dy,y0,ap)
         dy=pmx(26)
         y0=y-dy
         sr=.5*(sr+spp(y0,gamma,ec)*grxyl1(pm1,pmx,dy,y0,ap))
         do 62 i=1,ncp
            dy=dy+ddy
            y0=y-dy
            sr=sr+spp(y0,gamma,ec)*grxyl1(pm1,pmx,dy,y0,ap)
 62      continue
         sr=sr*ddy
         if(dym.gt.2.)then
            ddy=dym-2
            nc=int(ncp*ddy/(dym-pmx(26)))
            ddy=ddy/(nc+1)
            dy=dym
            y0=y-dy
            srp=spp(y0,gamma,ec)*grxyh1(pm1,pmx,dy,y0,ap)
            dy=2
            y0=y-dy
            srp=.5*(srp+spp(y0,gamma,ec)*grxyh1(pm1,pmx,dy,y0,ap))
            do 64 i=1,nc
               dy=dy+ddy
               y0=y-dy
               srp=srp+spp(y0,gamma,ec)*grxyh1(pm1,pmx,dy,y0,ap)
 64         continue
            sr=sr+srp*ddy
         endif
         spref(j)=(.5-fjc)*spref(j)+(.5+fjc)*sr*x(j)
 60   continue
c-----------------------
      do 30 j=jtransh+1,jrefmax
         y=1/x(j)
         dym=y-ymin
         dy=min(2,dym)
         ddy=(dy-pmx(26))/(ncp+1)
         y0=y-dy
         sr=spp(y0,gamma,ec)*grxyl1(pm1,pmx,dy,y0,ap)
         dy=pmx(26)
         y0=y-dy
         sr=.5*(sr+spp(y0,gamma,ec)*grxyl1(pm1,pmx,dy,y0,ap))
         do 32 i=1,ncp
            dy=dy+ddy
            y0=y-dy
            sr=sr+spp(y0,gamma,ec)*grxyl1(pm1,pmx,dy,y0,ap)
 32      continue
         sr=sr*ddy
         if(dym.gt.2.)then
            ddy=dym-2
            nc=int(ncp*ddy/(dym-pmx(26)))
            ddy=ddy/(nc+1)
            dy=dym
            y0=y-dy
            srp=spp(y0,gamma,ec)*grxyh1(pm1,pmx,dy,y0,ap)
            dy=2
            y0=y-dy
            srp=.5*(srp+spp(y0,gamma,ec)*grxyh1(pm1,pmx,dy,y0,ap))
            do 34 i=1,nc
               dy=dy+ddy
               y0=y-dy
               srp=srp+spp(y0,gamma,ec)*grxyh1(pm1,pmx,dy,y0,ap)
 34         continue
            sr=sr+srp*ddy
         endif
         spref(j)=sr*x(j)
 30   continue
c----------------------
      do 40 j=jrefmax+1,jmax
        spref(j)=0
 40   continue
 
      return
      end
 
c ------------------------------------------------------------------ c
c  Reflection factor for nonrelativistic Compton reflection,
c  reflection orders >1 approximated with isotropic phase function,
c  and approximation of H-function by Basko ApJ 223,268,(1978),
c  valid up to 14keV;
c  sig - sigma_bf per hydrogen atom in sigma_Thomson units
c  xm  - cosine of the angle between normal to the slab and the
c        observation direction
 
      real function    gnrj1(sig,xm)
      implicit         none
      real             sig,xm
      real             al0,al1,al2,xl,x2,x4,sq3,xq3
      parameter        (sq3=1.732051)
 
      al0=1/(1+sig/1.21)
      al1=sqrt(1-al0)
      al2=sq3*al1
      xl=log(1+1/xm)
      x2=xm*xm
      x4=x2*x2

      if(al0.lt.0.999) then 
        gnrj1=.5*al0*xm*(.375*((3-2*x2+3*x4)*xl+(3*x2-1)*(.5-xm))+
     & ((1+(1/al1-1)*(1-log(1+al2)/al2))*(1+xm*sq3)/(1+xm*al2)-1)*xl)
      else
        xq3=1+xm*sq3
        gnrj1=.5*al0*xm*(.375*((3-2*x2+3*x4)*xl+(3*x2-1)*(.5-xm))+
     & ((1+sq3/2)*xq3*(1-xq3*sqrt(1-al0))-1)*xl)
      endif
      return
      end
 
c ------------------------------------------------------------------- c
 
      function sigmaion(in,x,n,gamma0,temp0,xi0,ab_0,ab0,asc0)
c     bound-free cross section for ionized matter with cosmic abundances
c     and variable iron abundance (in units of solar iron abundance)
c     per hydrogen atom in units of the Thomson cross section
c     e is in units of 511 keV
      implicit none  
      integer nem 
      PARAMETER(nem=5000) 
      REAL ab,ab0,ab_,ab_0,asc,asc0,gamma,gamma0,SIGMAION,temp,
     &     Temp0,WARMABSJ,xas,xi,Xi0
      INTEGER i,In,n,nold 
      real x(nem),opac(nem),xold(nem)
      logical fl1,firstflag
      save opac,gamma,temp,xi,ab_,ab,asc,xold,nold,fl1,firstflag
      data xold/nem*0./
      data nold/0/
      data gamma,temp,xi,ab_,ab/0.,0.,0.,99999.,99999./
      data firstflag,fl1/2*.true./

c In this version, changes of ear during one xspec session are allowed.
      if((gamma.ne.gamma0).or.(temp.ne.temp0).or.
     &   (xi.ne.xi0).or.(ab_.ne.ab_0).or.(ab.ne.ab0))then
         gamma=gamma0
         temp=temp0
         xi=xi0
         ab=ab0
         ab_=ab_0
         firstflag=.true.
         fl1=.true.
         xas=.01947
         asc=1.503e2*
     &       warmabsj(xas,gamma,temp,xi,ab_,ab,firstflag)*xas**3
      endif

      if((x(in).ne.xold(in)).or.(n.ne.nold))then
         nold=n
         do i=1,n
            xold(i)=x(i)
         enddo
         fl1=.true.
      endif

      if(fl1)then
         do 10 i=1,n
            opac(i)=1.503e2*
     &              warmabsj(x(i),gamma,temp,xi,ab_,ab,firstflag)
 10      continue
         fl1=.false.
      endif

      asc0=asc
      sigmaion=opac(in)
 
      return
      end

c ------------------------------------------------------------------- c

	function warmabsj(x0,gamma,temp,xi,ab_,ab0,firstflag)
c 	calculate cross-sections using Daltabuit and Cox, use cosmic 
c       abundances, and radiative recombination rates from 
c       Aldrovandi and Pequignot 1973 and explicit formula 
c       for recombination to hydrogenlike atoms
c       result in units of 10^-22
        implicit none
        real     warmabsj,x0,gamma,temp,xi,ab_,ab0 
        real     t4,z2,e1,e2,re,rp,ab_d,it4 
	real     aw(20),abo(20),ab(20),ion(20,30,10),y
        real     int(20,30),arec(20,30)
	real     denom,xp,sigmaeff,e(800),sigma(20,30,800)
	real*8   num(20,30),sum,ratio(20,30),mult,mul(30),xil
	real*8   dnum,dmul 
	integer  z(20),i,j,imax,n,ios,k,iab0
	logical  firstflag
c
c     ******new bit for xspec
cjp
c for XSPEC 9
c      include '../../../../../spectral/inc/sitedef.inc'
c for XSPEC 10 when local directory in $XANADU/local 
c      include '../spectral/inc/sitedef.inc'
      character*72 absfil
      integer lenact,ilun10,ilun11,lenn  

      character fgdatd*128
      external lenact, fgdatd

c

c ------------------------------------------------------------------- c
c  set up constants and read in data if first time round

        save z,aw,abo,ab,ion,e,sigma,num,iab0,imax 
        data z(1)/0./

	if(z(1).eq.0) then
           iab0=0
        call getlun(ilun10)


       absfil = fgdatd()
       lenn = lenact(absfil)
       absfil = absfil(:lenn)//'iond_mm.dat'  

c        absfil = XSPECDR(:lenact(XSPECDR))//'iond_mm.dat'
cjp
        call openwr(ilun10, absfil, 'old', ' ', ' ', 0, 0, ios)
c          open(unit=10,file='iond_mm.dat',status='old',iostat=ios)
c          open(unit=11,file='mansig.dat',status='old')
c          ilun10=10
c          ilun11=11

	   i=1
	   read(ilun10,*,iostat=ios) z(i),aw(i),abo(i)
           if(z(i).eq.26) iab0=i
c  NB j=1 is ground state, j=z(i)=fully charged
	   do while (ios.eq.0)
	      do j=1,z(i),1
		 read(ilun10,*) ion(i,j,1),ion(i,j,2),ion(i,j,3),
     &               ion(i,j,4),ion(i,j,5),ion(i,j,6),ion(i,j,7),
     &               ion(i,j,8),ion(i,j,9),ion(i,j,10)
c  change units of coefficients
                 ion(i,j,2)=ion(i,j,2)*1.e+10
                 ion(i,j,4)=ion(i,j,4)*1.e+04
                 ion(i,j,5)=ion(i,j,5)*1.e-04
                 ion(i,j,7)=ion(i,j,7)*1.e-04
	      enddo
	      i=i+1
	      read(ilun10,*,iostat=ios) z(i),aw(i),abo(i)
              if(z(i).eq.26) iab0=i
	   enddo
	   imax=i-1

c       ******new bit for xspec
        close(unit=ilun10)
        call frelun(ilun10)

c  check not overflowing array bounds
           if(imax.gt.20.or.z(imax).gt.30) then 
              write(*,*) 'element array overflowing array bounds'
              write(*,*) 'number of elements= ',imax
              write(*,*) 'maximum atomic number= ',z(imax)
              stop
           end if
c  read in photoionization data from Reilman and Manson ApJ supp 40 1979
        call getlun(ilun11)

      absfil = fgdatd()
      lenn = lenact(absfil)
      absfil = absfil(:lenn)//'mansig.dat'

c        absfil = XSPECDR(:lenact(XSPECDR))//'mansig.dat'
        call openwr(ilun11, absfil, 'old', ' ', ' ', 0, 0, ios)

           do i=1,imax,1
c  now go over all ion states - values come from Reilman
c  ApJ supp 1979 except hydrogenic which is analytic
              do j=1,z(i),1
c  ion identification header
                 read(ilun11,*) 
c  read in over all energies 
c  the table is valid for 5.eV-19907.7eV
                 do k=1,721,1
                    read(ilun11,*) e(k),sigma(i,j,k)
c  change units of sigma
                    sigma(i,j,k)=sigma(i,j,k)/6.6e-27
                 enddo
              enddo
           enddo
c       ******new bit for xspec
        close(unit=ilun11)
        call frelun(ilun11)

c           close(ilun10)
c           close(ilun11)

	endif

c  end of data set up
c ------------------------------------------------------------------- c

        if(firstflag) then
c  insert external iron abundance (ab0)
           ab_d=10.**ab_
           do i=1,imax
              if(z(i).gt.2) ab(i)=abo(i)*ab_d
           enddo
           if(iab0.gt.0) ab(iab0)=ab(iab0)*10.**ab0
           if(xi.gt.1.e-20)then
c  normalized temperature
              t4=temp*.0001
c  calculate recombination coefficients - radiative and dielectric
c  from Aldrovandi+Pequignot 1973 Astro, astrophys 25 137
              do i=1,imax,1
                 do j=1,z(i)-1,1
                    it4=ion(i,j,5)/t4
                    e1=0. 
                    if(it4.lt.40.) e1=exp(-it4) 
                    it4=ion(i,j,7)/t4
                    e2=0. 
                    if(it4.lt.40.) e2=exp(-it4) 
                    arec(i,j)=ion(i,j,2)*(t4**(-ion(i,j,3)))+
     &	               ion(i,j,4)*(t4**(-1.5))*e1*(1.+ion(i,j,6)*e2)
                 enddo
c  do radiative recomb to hydrogenic ion NB assumed gaunt factor=1
c  from Gould+Thakur 1970 ann phys 61 351.
                 z2=z(i)**2
                 y=15.8*z2/t4
                 arec(i,z(i))=1.033e-3*z2*(1.735+log(y)+1/(6.*y))
     &              /sqrt(t4)
              enddo
c  first normalise photon spectrum 
              denom=0.
              do k=1,720,1
                 denom=denom+(e(k)**(1-gamma))*(e(k+1)-e(k))
              enddo
              denom=1/denom
c  calculate the photoionization integrals and population balance
              do i=1,imax,1
                 do j=1,z(i),1
                    int(i,j)=0.0
c  do integral in log10 (eV) from lowest energy
c  to the maximum energy- use same energy grid
c  as the cross-sections are defined on
                    do k=1,720,1
                       int(i,j)=int(i,j)+
     &                    sigma(i,j,k)*(e(k)**(-gamma))*(e(k+1)-e(k))
                    enddo
                    int(i,j)=int(i,j)*denom
c  const=7.958d-2/2.43d+4 
		    ratio(i,j)=dlog(3.2749d-6*int(i,j)/arec(i,j))
	         enddo
	      enddo
c  calculate population numbers
c  set sum=1 as the sum=1+A21+A32A21+....
c  complicated because floating range limit 
              do i=1,imax,1
                 xil=dlog(dble(xi))
                 mult=0.d0
                 do j=z(i),1,-1
                    mul(j)=0.d0
                    do n=j,1,-1
                       mul(j)=mul(j)+ratio(i,n)
                    enddo
                    mul(j)=mul(j)+j*xil
                    if(mul(j).gt.mult) mult=mul(j)
                 enddo
                 do j=z(i),1,-1
                    mul(j)=mul(j)-mult
                 enddo
                 sum=0.d0
                 do j=z(i),1,-1
                    dmul=0.d0 
                    if(dabs(mul(j)).lt.7d1) dmul=dexp(mul(j)) 
                    sum=sum+dmul
                 enddo
                    dmul=0.d0 
                    if(dabs(mult).lt.7d1) dmul=dexp(-mult) 
                 sum=sum+dmul
                 num(i,1)=-mult-dlog(sum)
                 do j=2,z(i)+1,1
                    num(i,j)=num(i,j-1)+ratio(i,j-1)+xil
                 enddo
                 do j=1,z(i)+1,1
                    dnum=0d0 
                    if(dabs(num(i,j)).lt.7d1) dnum=dexp(num(i,j))
                    num(i,j)=dnum
                 enddo
	      enddo				
           else
              do i=1,imax,1
                 num(i,1)=1.
                 do j=2,z(i)+1,1
                    num(i,j)=0.
                 enddo
              enddo
           endif
           firstflag=.false.
        endif
c ------------------------------------------------------------------- c

        xp=x0*5.11e5
	sigmaeff=0.
        if(xp.lt.e(721)) then
c  pick closest energy bin
           k=1
           dowhile(e(k).lt.xp)
              k=k+1
           enddo
           re=(xp-e(k-1))/(e(k)-e(k-1))
           do i=1,imax,1
              do j=1,z(i),1
                 sigmaeff=sigmaeff+sngl(num(i,j))*ab(i)*
     &              (sigma(i,j,k-1)+re*(sigma(i,j,k)-sigma(i,j,k-1)))
              enddo
           enddo
        else
           re=(xp/e(721))**(-3)
           do i=1,imax,1
              rp=ab(i)*re
              do j=1,z(i),1
                 sigmaeff=sigmaeff+sngl(num(i,j))*sigma(i,j,721)*rp
              enddo
           enddo
        endif
        sigmaeff=sigmaeff*6.6e-5

        warmabsj=sigmaeff

	return
	end

c ------------------------------------------------------------------- c

