SUBROUTINE REFxion(Ear,Ne,Param,Ifl,Photar,photer)
 
      IMPLICIT NONE
      INTEGER Ne , Ifl
      REAL Ear(0:Ne) , Param(*) , Photar(Ne), photer(Ne)
 
c     Driver for angle-dependent reflection from a neutral medium.
c     See Magdziarz & Zdziarski, 1995, MNRAS.
c
c     The output spectrum is the sum of the input (the original contents
c     of Photar) and the 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     Version with variable iron abundance and new opacities of Balucinska
c     & McCammon (1992, and 1994, private communication). As expected in AGNs,
c     H and He are assumed to be fully ionized.
c
c     This version allows for changes of the vector 'ear' between subsequent
c       calls.
c
c
c     number of model parameters: 6
c     1: scale, scaling factor for reflection; if <0, no direct component
c        (scale=1 for isotropic source above disk)
c     2: redshift, z
c     3: afe
c     4: cosine of inclination angle
c     algorithm:
c          a[E*(1+z)]=E^{-Gamma}*exp(-E/E_c)+scale*reflection
c     Normalization is the photon flux at 1 keV (photons keV^-1 cm^-2 s^-1)
c     of the cutoff power law only (without reflection)
c     and in the earth frame.
c
c     ballantyne  parameters
c     5: log10(Xi)
c     6: gamma

      INCLUDE './xspec.inc'

      INTEGER isptot, ispref, ispinc, ix, iphotaux

      INTEGER i , j , ierr, nesave
      LOGICAL firstcall
      CHARACTER contxt*127

      SAVE isptot, ispref, ispinc, ix
      SAVE firstcall 
      SAVE nesave

      DATA isptot, ispref, ispinc, ix, nesave/5*-1/
      DATA firstcall/.TRUE./
      
      ierr = 0
 

      IF ( firstcall ) THEN
         CALL xwrite('Compton reflection from neutral medium.',5)
         CALL xwrite('If you use results of this model in a paper,',5)
         CALL xwrite
     &  ('please refer to chris'
     &              ,5)
         firstcall = .FALSE.
      ENDIF

c     check on array
      if (photar(ne/2).ne.0.0) then

      IF ( Ne .NE. nesave ) THEN
         CALL udmget(Ne, 6, isptot, ierr)
         contxt = 'Failed to get memory for sptot'
         IF ( ierr .NE. 0 ) GOTO 999
         CALL udmget(Ne, 6, ispref, ierr)
         contxt = 'Failed to get memory for spref'
         IF ( ierr .NE. 0 ) GOTO 999
         CALL udmget(Ne, 6, ispinc, ierr)
         contxt = 'Failed to get memory for sppow'
         IF ( ierr .NE. 0 ) GOTO 999
         CALL udmget(Ne, 6, ix, ierr)
         contxt = 'Failed to get memory for x'
         IF ( ierr .NE. 0 ) GOTO 999
         nesave = Ne
      ENDIF

c   x is in the source frame and in units of m_e c^2.

      DO 100 i = 1 , Ne
         MEMR(ix+i-1) = ((Ear(i)+Ear(i-1))/2./511.)*(1.+Param(2))
 100  CONTINUE

C    Generate the spinc array from the input photar array.

      DO i = 1, Ne
         IF ( ear(i) .NE. ear(i-1) ) THEN
            MEMR(ispinc+i-1) = photar(i) * ((ear(i)+ear(i-1))/2.)**2
     &                          / (ear(i)-ear(i-1))
         ENDIF
      ENDDO

c     x is the energy array (units m_e c^2)
c     spinc is the input spectrum array (E F_E)
c     spref is the reflected spectrum array (E F_E)
c     sptot is the total spectrum array (E F_E), = spinc if no reflection
c     all dimensions = Ne



      CALL b2DOREFAx(Param,MEMR(ix),Ne,MEMR(isptot),MEMR(ispinc),
     &        MEMR(ispref))

      DO i = 1 , Ne
         Photar(i) = MEMR(isptot+i-1) * (ear(i)-ear(i-1))
     &                     / ((ear(i)+ear(i-1))/2.)**2
      ENDDO

 999  CONTINUE

      IF ( ierr .NE. 0 ) THEN
         CALL xwrite(contxt, 10)
      ENDIF

      end if
c     end fo array check

      RETURN
      END
**==powrefa.spg  processed by SPAG 4.50J  at 11:22 on 27 Feb 1996
 
 
      SUBROUTINE b2DOREFAx(Param,X,Jmax,Sptot,Spinc,Spref)
c
      IMPLICIT NONE
      REAL gamma , ec , scale 
      REAL Param(*) , Sptot(*) , Spref(*) , Spinc(*) , X(*)
      INTEGER Jmax , j
 
c--------
C  INPUT:
c the scaling factor for the reflected spectrum;
c 1 corresponds to seeing equal
c contributions from the reflected and direct spectra
      scale = Param(1)
c redshift
c      z_red = Param(2)
c Iron abundance
c      ab0 = LOG10(Param(3))
c Cosine of inclination angle
c      xm = Param(4)
c--------
 
      IF ( scale.NE.0. ) THEN
         CALL b2DOGREENx(param,x,Jmax,Spinc,Spref)
         DO 120 j = 1 , Jmax
            Spref(j) = Spref(j)*ABS(scale)
 120     CONTINUE
      ELSE
         DO 140 j = 1 , Jmax
            Spref(j) = 0
 140     CONTINUE
      ENDIF
      IF ( scale.GE.0 ) THEN
         DO 160 j = 1 , Jmax
            Sptot(j) = Spinc(j) + Spref(j)
 160     CONTINUE
      ELSE
         DO 180 j = 1 , Jmax
            Sptot(j) = Spref(j)
 180     CONTINUE
      ENDIF
 
      RETURN
      END
 
c ------------------------------------------------------------------ c
 
      SUBROUTINE b2DOGREENx(param,X,Jmax,Spinc,Spref)
c     calculates the reflected spectrum

      IMPLICIT NONE
      REAL ddy , dy , dym , Gamma, GRXYH, 
     &     GRXYL , sr , srp , xjd , xjl , xkabs , Xm , 
     &     xmin , xmo
      REAL xnor , xrefmax , xtransh , xtransl , y , y0 , ymin
      INTEGER i , j , Jmax , jmin , jrefmax , jtransh , jtransl , nc , 
     &        ncp
      REAL param(*), X(*) , Spref(*), Spinc(*)
      REAL pm1(19) , pmx(26) , ap(3) 
      REAL SPINTERP
      EXTERNAL SPINTERP
      real xphotar(10000),xear(0:10000),xpar(4),xphoter(10000)
      real sprefp(10000), xphotar1(10000),xerr(10000)
      real scale,z_red,xi,afe
      real norm,ec,e,e1,sum,dloge,de
      integer j2,j10,ifl,count
      

      character*255 filenm
      real xnors,xkabs_pex,gamma_ref,gamma_pex
      real factor,delta
      save xkabs

      data xkabs/1.0/


c     parameters
      scale=param(1)
      z_red=param(2)
      afe=param(3)
      xm=param(4)
      xi=10**param(5)
      gamma=param(6)

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
      xtransl=12.0/511.0
      xtransh=14.0/511.0
      ymin = .03333
c  angle dependent parameters
      xmo = MAX(.05,MIN(.95,Xm))
      CALL PM1Y(xmo,pm1)
      CALL PMXY(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 100 j = 1 , Jmax
         IF ( X(j).GT.xmin ) GOTO 200
         jmin = j
 100  CONTINUE
 200  jtransl = 0
      DO 300 j = MAX(1,jmin) , Jmax
         IF ( X(j).GT.xtransl ) GOTO 400
         jtransl = j
 300  CONTINUE
 400  jtransh = 0
      DO 500 j = MAX(1,jtransl) , Jmax
         IF ( X(j).GT.xtransh ) GOTO 600
         jtransh = j
 500  CONTINUE
 600  jrefmax = 0
      DO 700 j = MAX(1,jtransh) , Jmax
         IF ( X(j).GT.xrefmax ) GOTO 800
         jrefmax = j
 700  CONTINUE
 
 800  DO 900 j = 1 , jmin
         Spref(j) = 0
 900  CONTINUE

c     afe
      xpar(1)=afe

c     take average 2-10 keV index
c      do j=1,jmax,1
c         if ((x(j).gt.3.9e-3).and.(x(j-1).le.3.9e-3)) j2=j
c         if ((x(j).ge.1.9e-2).and.(x(j-1).lt.1.9e-2)) j10=j
c      end do
      
c      xpar(2)=log10(spinc(j10)/(x(j10)**2))-log10(spinc(j2)/(x(j2)**2))
c      xpar(2)=xpar(2)/(log10(x(j10))-log10(x(j2)))
c      xpar(2)=-xpar(2)
c      gamma=xpar(2)
      xpar(2)=gamma
      
c     reflected fraction - set to zero for continuum, 1 for reflection
      xpar(3)=xi

c     redshift
      xpar(4)=0.0

c     set up energy grid again in keV
      if (jmax.gt.10000) then
         write(*,*) 'need bigger arrays'
         stop
      end if
      do j=0,jmax,1
         xear(j)=x(j)*511.0
      end do

      ifl=1
c     reflion
      filenm='./reflionx.mod'
      call xsatbl(xear, jmax, xpar, filenm, ifl, xphotar1, xerr)

c     continuum
      sum=0.0
      ec=300.0
      dloge=log10(1e3/1e-3)/1000.0
      e1=1.0e-3
      do i=1,1000,1
         e=10**(-3.0+dloge*(float(i)-0.5))
         de=10**(-3.0+dloge*float(i))-e1
         e1=e1+de
         sum=sum+(e**(-gamma+1.0))*exp(-e/ec)*de
      end do
      norm=xi*1e15*sum*1e-4/(4.0*3.1415)

c     convert it to factors
      do j=1,jmax,1
         e=10**(0.5*(log10(xear(j))+log10(xear(j-1))))
         xphotar(j)=norm*(e**(-gamma))*exp(-e/ec)
         xphotar(j)=xphotar(j)*(xear(j)-xear(j-1))
         if (xear(j).lt.30.0) then
            xphotar(j)=xphotar1(j)/xphotar(j)
            else
            xphotar(j)=0.0
         end if
         if (xphotar(j).lt.0.0) xphotar(j)=0.0
      end do

c-----------------------
c     then can do reflection - do up to 14 keV
      DO 1000 j = jmin + 1 , jtransh
         Spref(j) = Spinc(j)*xphotar(j)
 1000 CONTINUE
c-----------------------

c     take average 10-12 keV photon index of reflection
      gamma_ref=log10(spref(jtransh)/(x(jtransh)**2))-
     &   log10(spref(jtransl+1)/(x(jtransl+1)**2))
      gamma_ref=gamma_ref/(log10(x(jtransh))-log10(x(jtransl+1)))
      gamma_ref=-gamma_ref

c     albedo-need to iterate. start at high cx so little reflection at xtrans
c     do in log space
      gamma_pex=2.0*gamma_ref

      count=0

c     find the pexriv model with the same slope
      do while ( ((abs(gamma_ref-gamma_pex))/(abs(gamma_ref)).gt.1.e-3)
     &   .and.(count.lt.40) ) 

c       keep track
        count=count+1

c       need to find crosssection  at the 10 keV point
c       log crosssection set to 1.0 first time around, otherwise
c       remembers the old value
        xkabs=10**xkabs
        xnor=xkabs*((10.0/511.0)**3)
        xkabs=log10(xkabs)

        ap(1) = xnor/1.21
        xjl = ALOG(xtransl)
        xjd = ALOG(xtransh) - xjl
        DO 1100 j = jtransl + 1 , jtransh
c          fjc = .5*SIN(3.14159*((ALOG(X(j))-xjl)/xjd-.5))
c          Spref(j) = Spinc(j)*xphotar(j)
          y = 1/X(j)
          dym = y - ymin
          dy = MIN(2.,dym)
          ddy = (dy-pmx(26))/(ncp+1)
          y0 = y - dy
          sr = SPINTERP(y0,jmax,X,Spinc)*GRXYL(pm1,pmx,dy,y0,ap)
          dy = pmx(26)
          y0 = y - dy
          sr = .5*(sr+SPINTERP(y0,jmax,X,Spinc)*GRXYL(pm1,pmx,dy,y0,ap))
          DO 1050 i = 1 , ncp
            dy = dy + ddy
            y0 = y - dy
            sr = sr + SPINTERP(y0,jmax,X,Spinc)*GRXYL(pm1,pmx,dy,y0,ap)
 1050     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 = SPINTERP(y0,jmax,X,Spinc)*GRXYH(pm1,pmx,dy,y0,ap)
            dy = 2
            y0 = y - dy
            srp = .5*(srp+SPINTERP(y0,jmax,X,Spinc)
     &                    *GRXYH(pm1,pmx,dy,y0,ap))
            DO 1060 i = 1 , nc
               dy = dy + ddy
               y0 = y - dy
               srp = srp + SPINTERP(y0,jmax,X,Spinc)
     &                     *GRXYH(pm1,pmx,dy,y0,ap)
 1060       CONTINUE
            sr = sr + srp*ddy
          ENDIF
c         Spref(j) = (.5-fjc)*Spref(j) + (.5+fjc)*sr*X(j)
          Sprefp(j) = sr*X(j)
 1100   CONTINUE

c       take average 12-14 keV index
        gamma_pex=log10(sprefp(jtransh)/(x(jtransh)**2))-
     &   log10(sprefp(jtransl+1)/(x(jtransl+1)**2))
        gamma_pex=gamma_pex/(log10(x(jtransh))-log10(x(jtransl+1)))
        gamma_pex=-gamma_pex

c       if gamma_pex is smaller then need more ionisation so need less 
c       crossection. and converse.
        xkabs=xkabs-0.3*(gamma_ref-gamma_pex)

      end do 

c     renormalise below 14 keV by the difference at 14 keV.
      factor=sprefp(jtransh)/spref(jtransh)

c     catch it when it falls over
      if (factor.ne.abs(factor)) then
         write(*,*) 'NAN detected',sprefp(jtransh),spref(jtransh),factor
         write(*,*) gamma,gamma_ref,gamma_pex,x(1)*511,x(jmax)*511,xkabs
         stop
      end if
      if (spref(jtransh).eq.0.0) then
         write(*,*) 'INF detected',sprefp(jtransh),spref(jtransh),factor
         write(*,*) gamma,gamma_ref,gamma_pex,x(1)*511,x(jmax)*511,xkabs
         stop
      end if

      do j=jmin+1,jtransh,1
         spref(j)=factor*spref(j)
      end do

c-----------------------
      DO 1200 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 = SPINTERP(y0,jmax,X,Spinc)*GRXYL(pm1,pmx,dy,y0,ap)
         dy = pmx(26)
         y0 = y - dy
         sr = .5*(sr+SPINTERP(y0,jmax,X,Spinc)*GRXYL(pm1,pmx,dy,y0,ap))
         DO 1150 i = 1 , ncp
            dy = dy + ddy
            y0 = y - dy
            sr = sr + SPINTERP(y0,jmax,X,Spinc)*GRXYL(pm1,pmx,dy,y0,ap)
 1150    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 = SPINTERP(y0,jmax,X,Spinc)*GRXYH(pm1,pmx,dy,y0,ap)
            dy = 2
            y0 = y - dy
            srp = .5*(srp+SPINTERP(y0,jmax,X,Spinc)
     &                    *GRXYH(pm1,pmx,dy,y0,ap))
            DO 1160 i = 1 , nc
               dy = dy + ddy
               y0 = y - dy
               srp = srp + SPINTERP(y0,jmax,X,Spinc)
     &                     *GRXYH(pm1,pmx,dy,y0,ap)
 1160       CONTINUE
            sr = sr + srp*ddy
         ENDIF
         Spref(j) = sr*X(j)
 1200 CONTINUE
c----------------------
      DO 1300 j = jrefmax + 1 , Jmax
         Spref(j) = 0
 1300 CONTINUE
 
c     first one looks dire
      spref(1)=spref(2)

      RETURN
      END


      function spinterp(invx, jmax, x, spinc)

      INTEGER jmax
      REAL spinterp, invx, x(*), spinc(*)

c Performs binary search on spinc given array x and input invx (=1/x).
c presumes monotonic increasing array x
c efficiency improvements BD 9/2003

      INTEGER ix, ixlo, ixhi
      REAL xtarg
      DATA ixlo/1/ixhi/1/


      xtarg = 1./invx

      IF ( xtarg .LT. x(1) ) THEN
         spinterp = 0.
         RETURN
      ENDIF

      IF ( xtarg .GT. x(jmax) ) THEN
         spinterp = 0.
         RETURN
      ENDIF


      if (xtarg.gt.x(ixlo).and.xtarg.lt.x(ixhi)) then
           continue
      else if ( xtarg.gt.x(ixhi).and.xtarg.lt.x(ixhi+1)) then
c we already dealt with ixhi == jmax
           ixlo = ixhi
           ixhi = ixhi + 1
      else if ( xtarg.lt.x(ixlo).and.xtarg.gt.x(ixlo - 1)) then
c ... and the case where ixlo == 1
           ixhi = ixlo
           ixlo = ixhi - 1 
      else

              ixlo = 0
              ixhi = jmax + 1

              DO WHILE ( (ixhi-ixlo) .GT. 1 )

                 ix = (ixhi + ixlo) / 2
                 IF ( xtarg .GT. x(ix) ) THEN
                    ixlo = ix
                 ELSE
                    ixhi = ix
                 ENDIF

              ENDDO
      end if

      IF ( ixlo .NE. ixhi ) THEN
         spinterp = ( spinc(ixlo)*(x(ixhi)-xtarg) + 
     &                spinc(ixhi)*(xtarg-x(ixlo)) ) / (x(ixhi)-x(ixlo))
      ELSE
         spinterp = spinc(ixlo)
      ENDIF

      RETURN
      END