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 
      save firstflag,x,phot_ref,phot_disk,in_disk,
     >     phot_rs,x_disk,spref,dx,logxmin,logxmax,logdmin,logdmax
      data firstflag/.true./,fl1/.true./ 
      data eps/1e-5/,ec/0./,ifl/0/   
      data in_disk/ndisk*9999./

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 broadening 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 greencj(x,nx,spref,xm,ab_,ab0,sprf,gamma,ec,temp,xi)
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