
      SUBROUTINE irflct(ear, ne, param, ifl, photar, photer)

      INTEGER ne, ifl
      REAL  ear(0:ne), param(*), photar(ne), photer(ne)

c Convolution subroutine for angle-dependent reflection based on the
c Magdziarz and Zdziarski code. The photar array is assumed to contain
c the input spectrum to be reflected. On output it contains the sum of
c the incident and reflected. It is probably necessary to run the extend
c command before using this model since the incident spectrum should be
c calculated to significantly higher energies than covered by the output
c spectrum.

c
c     See Magdziarz & Zdziarski, 1995, MNRAS.
c     See Zdziarski et al., 1995, ApJL January 10 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 incident spectrum 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
c     This version allows for changes of the vector 'ear' between subsequent
c       calls.
c
c
c     number of model parameters:7
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: abundance of elements heavier than He relative to 
c        the solar abundances 
c     4: iron abundance relative to the above  
c     5: cosine of inclination angle
c     6: reflector temperature in K
c     7: reflector ionization parameter = L/nR^2
c     algorithm:
c          a[E*(1+z)]=incident spectrum+scale*reflection
c     Normalization is the photon flux at 1 keV (photons keV^-1 cm^-2 s^-1)
c     of the incident spectrum only (without reflection)
c     and in the earth frame.

      INCLUDE 'xspec.h'

      INTEGER isptot, ispref, ispinc, ix
      SAVE isptot, ispref, ispinc, ix      

      INTEGER i, ierr, nesave
      CHARACTER contxt*72
      SAVE nesave

      logical firstcall
      save firstcall

      DATA isptot, ispref, ispinc, ix, nesave/5*-1/
      data firstcall/.true./

      ierr = 0

      if (firstcall) then
        call xwrite('Compton reflection from ionized medium',5)
        call xwrite('rel_refl<0 gives the reflected spectrum alone',5)
        call xwrite('assumes H and He are fully ionised',5)
        firstcall = .false.
      end if

c Grab the memory for the following temporary arrays 
c     x is the energy array (units m_e c^2)
c     spref is the reflected spectrum array (E F_E)
c     spinc is the incident spectrum array (E F_E)
c     sptot is the total spectrum array (E F_E), = spinc if no reflection
c     all dimensions = ne   

      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 spinc'
         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 Set array of energies as fractions of electron-positron annihilation energy

      DO i = 0, ne
         MEMR(ix+i) = ((Ear(i)+Ear(i-1))/2./511.)*(1.+Param(2))
      ENDDO

c Generate the spinc array from the input photar array.

      DO i = 1, Ne
         MEMR(ispinc+i-1) = photar(i) * (ear(i)+ear(i-1))/2.
      ENDDO

c Do the reflection

      CALL DOIRFA(param, 7, MEMR(ix), Ne, MEMR(isptot), 
     &            MEMR(ispref), MEMR(ispinc))

c Construct the output photar array from sptot.

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

 999  CONTINUE
      IF ( ierr .NE. 0 ) THEN
         CALL xwrite(contxt, 10)
      ENDIF
 
      RETURN
      END


      SUBROUTINE DOIRFA(param, npar, x, jmax, sptot, spref, spinc)
c
      IMPLICIT NONE
      real scale
      real param(*), sptot(*), spref(*), x(*), spinc(*)
      integer jmax,j,npar
 
c--------
      scale=param(1)
c--------

      if(scale.ne.0.)then
         call nngreencj(x,jmax,spref,param,npar,spinc)
         do 20 j=1,jmax
            spref(j)=spref(j)*abs(scale)
 20      continue
         else
         do 34 j=1,jmax
            spref(j)=0.0
 34      continue
      endif
      if (scale.gt.0) then
         do 22 j=1,jmax
            sptot(j)=spinc(j)+spref(j)
 22      continue
         else
         do 23 j=1,jmax
            sptot(j)=spref(j)
  23     continue
      endif
 
      return
      end

c ------------------------------------------------------------------ c
 
      subroutine nngreencj(x,jmax,spref,param,npar,spinc)
c     calculates the reflected spectrum
c     abfe - iron abundance relative to hydrogen
      implicit none 
      REAL ddy,dy,dym,fjc,GNRJ,GRXYH,
     &     GRXYL,NXSIGMABFJ,sr,srp,xjd,xjl,
     &     xkabs,Xm
      REAL xmin,xmo,xnor,xrefmax,xtransh,xtransl,y,y0,ymin 
      INTEGER i,j,Jmax,jmin,jrefmax,jtransh,jtransl,nc,ncp 
      integer npar
      real x(*),spref(*),param(*),spinc(*)
      real pm1(19),pmx(26),ap(3)
      REAL spinterp
      EXTERNAL spinterp

c     use inclination angle
      xm=param(5)

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 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 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=nxsigmabfj(j,x,jtransh,param,npar,xnor,jmax,spinc)
         spref(j)=spinc(j)*gnrj(xkabs,xmo)
 20   continue
c-----------------------
      if(jmin.ge.jtransl) then
         xkabs=nxsigmabfj(1,x,max(1,jtransh),param,npar,xnor,jmax,spinc)
      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=nxsigmabfj(j,x,jtransh,param,npar,xnor,jmax,spinc)
         y=1/x(j)
         spref(j)=spinc(j)*gnrj(xkabs,xmo)
         dym=y-ymin
         dy=min(2.0,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 62 i=1,ncp
            dy=dy+ddy
            y0=y-dy
            sr=sr+spinterp(y0,jmax,x,spinc)*grxyl(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=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 64 i=1,nc
               dy=dy+ddy
               y0=y-dy
               srp=srp+spinterp(y0,jmax,x,spinc)*grxyh(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.0,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 32 i=1,ncp
            dy=dy+ddy
            y0=y-dy
            sr=sr+spinterp(y0,jmax,x,spinc)*grxyl(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=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 34 i=1,nc
               dy=dy+ddy
               y0=y-dy
               srp=srp+spinterp(y0,jmax,x,spinc)*grxyh(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
 
      function nxsigmabfj(in,x,n,param0,npar,asc0,jmax,spinc)

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  

      INCLUDE 'xspec.h'

      INTEGER iopac, ixold

      REAL asc,asc0,NXSIGMABFJ,NXWARMABJ,xas
      real x(*),spinc(*),param0(*),param(30)
      INTEGER i,In,n,nold,npar,ierr,jmax
      CHARACTER contxt*72
      CHARACTER*4 fgsolr, fgsolr0
      EXTERNAL fgsolr

      logical fl1,firstflag

      save iopac,param,asc,ixold,nold,fl1,firstflag,fgsolr0

      data iopac /0/
      data ixold/-1./
      data nold/0/
      data firstflag,fl1/2*.true./
      DATA fgsolr0/'    '/

      if (npar.gt.30) then
         write(*,*) 'too many parameters'
         stop
      end if

c In this version, changes of ear during one xspec session are allowed.

      if(fgsolr().ne.fgsolr0)then
         fl1=.true.
         fgsolr0=fgsolr()
      endif

      do i=1,npar,1
         if (param(i).ne.param0(i)) firstflag=.true.
      end do

      if (firstflag) then 
         do i=1,npar,1
            param(i)=param0(i)
         end do
         fl1=.true.
         xas=.01947
         asc=1.503e2*
     &       nxwarmabj(xas,param,npar,firstflag,jmax,x,spinc)*xas**3
      endif
      if((x(in).ne.MEMR(ixold+in-1)).or.(n.ne.nold))then
         CALL udmget(N, 6, iopac, ierr)
         contxt = 'Failed to get memory for opac'
         IF ( ierr .NE. 0 ) GOTO 999
         CALL udmget(N, 6, ixold, ierr)
         contxt = 'Failed to get memory for xold'
         IF ( ierr .NE. 0 ) GOTO 999
         nold=n
         do i=1,n
            MEMR(ixold+i-1)=x(i)
         enddo
         fl1=.true.
      endif

      if(fl1)then
         do 10 i=1,n
            MEMR(iopac+i-1)=1.503e2*
     &              nxwarmabj(x(i),param,npar,firstflag,jmax,x,spinc)
 10      continue
         fl1=.false.
      endif

      asc0=asc
      nxsigmabfj=MEMR(iopac+in-1)

 999  CONTINUE
      IF ( ierr .NE. 0 ) THEN
         CALL xwrite(contxt, 10)
         NXSIGMABFJ = 0.
      ENDIF
 
      return
      end

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

c!! need to modify this to pass in energies on which spinc is evaluated.

      FUNCTION NXWARMABJ(x0,param,npar,firstflag,jmax,x,spinc)

      INTEGER npar, jmax
      REAL NXWARMABJ, x0, param(npar), x(*), spinc(*)
      LOGICAL firstflag

      INCLUDE 'xspec.h'

      REAL nxclcwrm
      INTEGER isigma, iion, iint, iarec, inum, iratio, ierr
      CHARACTER*72 contxt

      SAVE isigma, iion, iint, iarec, inum, iratio

      DATA isigma, iion, iint, iarec, inum, iratio /6*-1/

      ierr = 0

C Get the memory for the arrays required by NXCLCWRM

      IF ( isigma .EQ. -1 ) THEN
         CALL udmget(20*30*800, 6, isigma, ierr)
         contxt = 'Failed to get memory for sigma array'
         IF ( ierr .NE. 0 ) GOTO 999
         CALL xwrite('Memory acquired for sigma array', 25)
      ENDIF
      IF ( iion .EQ. -1 ) THEN
         CALL udmget(20*30*10, 6, iion, ierr)
         contxt = 'Failed to get memory for ion array'
         IF ( ierr .NE. 0 ) GOTO 999
         CALL xwrite('Memory acquired for ion array', 25)
      ENDIF
      IF ( iint .EQ. -1 ) THEN
         CALL udmget(20*30, 6, iint, ierr)
         contxt = 'Failed to get memory for int array'
         IF ( ierr .NE. 0 ) GOTO 999
         CALL xwrite('Memory acquired for int array', 25)
      ENDIF
      IF ( iarec .EQ. -1 ) THEN
         CALL udmget(20*30, 6, iarec, ierr)
         contxt = 'Failed to get memory for arec array'
         IF ( ierr .NE. 0 ) GOTO 999
         CALL xwrite('Memory acquired for arec array', 25)
      ENDIF
      IF ( inum .EQ. -1 ) THEN
         CALL udmget(20*30, 7, inum, ierr)
         contxt = 'Failed to get memory for num array'
         IF ( ierr .NE. 0 ) GOTO 999
         CALL xwrite('Memory acquired for num array', 25)
      ENDIF
      IF ( iratio .EQ. -1 ) THEN
         CALL udmget(20*30, 7, iratio, ierr)
         contxt = 'Failed to get memory for ratio array'
         IF ( ierr .NE. 0 ) GOTO 999
         CALL xwrite('Memory acquired for ratio array', 25)
      ENDIF

c Call the routine to do the calculation

      nxwarmabj = nxclcwrm(x0,param,npar,firstflag,jmax,x,spinc,
     &                     MEMR(isigma),MEMR(iion),MEMR(iint),
     &                     MEMR(iarec),MEMR(inum),MEMR(iratio))

 999  CONTINUE
      IF ( ierr .NE. 0 ) THEN
         CALL xwrite(contxt, 10)
         NXWARMABJ = 0.
      ENDIF

      RETURN
      END
 

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

	function nxclcwrm(x0,param,npar,firstflag,jmax,x,spinc,sigma,ion,
     &                    int,arec,num, ratio)

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

        integer npar, jmax
        real     nxclcwrm,x0,xi,temp,xe
        real     t4,z2,e1,e2,re,rp,param(npar)
	real     aw(20),abo(20),ab(20),ion(20,30,10),y
        real     int(20,30),arec(20,30),x(*),spinc(*),andgrv(10)
	real     denom,xp,sigmaeff,e(800),sigma(20,30,800)
	real*8   num(20,30),sum,ratio(20,30),mult,mul(30),xil
	integer  z(20),i,j,imax,n,ios,k,iab0
        integer lenact,ilun10,ilun11
	logical  firstflag
        character*72 absfil
        character*4 fgsolr, fgsolr0
        character el_name*2

        real spinterp, fgabnd
        external fgsolr, el_name, spinterp, fgabnd
        
        save z,aw,abo,ab,e,iab0
        data z(1)/0./
        data fgsolr0/'    '/

C If necessary load the arrays from the iond_mm.dat and mansig.dat
C files

      CALL ldapmr(z, aw, abo, ion, sigma, e, iab0, imax)

      xi=param(7)
      temp=param(6)
      
C initialize other arrays

      if((fgsolr().ne.fgsolr0).or.firstflag)then
         fgsolr0=fgsolr()
         do i=1,imax
            abo(i)=fgabnd(el_name(z(i)))
         enddo
c  insert external abundances (ab_, ab0)
         DO i = 1 , imax
            IF ( z(i).GT.2 ) THEN
               ab(i) = abo(i)*param(3)
            ELSE
               ab(i) = 0.
            ENDIF
         ENDDO
         IF ( iab0.GT.0 ) ab(iab0) = abo(iab0)*param(4)
      endif


         if ( firstflag ) then
           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
                    e1=exp(-ion(i,j,5)/t4)
                    e2=exp(-ion(i,j,7)/t4)
                    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
                 xe=e(k)/511.0e3
c!! fix this
                 denom=denom+(spinterp(1/xe,jmax,x,spinc)/xe)*
     &                (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
                       xe=e(k)/511.0e3
                       int(i,j)=int(i,j)+
c!! fix this
     &                    sigma(i,j,k)*(spinterp(1/xe,jmax,x,spinc)/xe)
     &                      /e(k)*(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
                    sum=sum+dexp(mul(j))
                 enddo
                 sum=sum+dexp(-mult)
                 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
                    num(i,j)=dexp(num(i,j))
                 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

        nxclcwrm=sigmaeff

	return
	end
