      subroutine spin(EAR,NE,PARAM,IFL,PHOTAR,PHOTER)
c
c     XSPECv11 Subroutine to produce an emission line profile from a thin
c     keplerian disk around a Kerr black hole with arbitrary spin
c     parameter.   Method relies on the transfer function approach
c     of Cunningham (1975); the line profile can be computed via a
c     simple integral with a analytic integrand apart from the
c     effects of gravitational light bending/lensing.  This are
c     incorporated into a slowly varying transfer function computed
c     via the code of Speith et al. (1995) and sparsely tabulated in the
c     accompanying table "kerrtable.dat".   High quality results are
c     obtained despite the sparse sampling of the TF through
c     linear interpolation of this slowing varying function.
c
c     This code also generates the basic relativistic disk kernel
c     used for the convolution model kerrconv.
c
c     Two features have been hardwired into this code but are simple
c     to change.  Firstly, we assume that the line emissivity is described
c     with a broken powerlaw (see comment starting "EMISSIVITY PROFILE").
c     It is trivial to include any functional form.  Secondly, we have
c     used the same limb-darkening rule as used in Laor (1991; see
c     comment starting "LIMB DARKENING").  This is readily changed to any
c     functional form.
c
c     At the current time, kerrdisk and kerrconv do not allow for emission
c     within the innermost stable circular orbit
c
c     For other details of this code see Brenneman & Reynolds, 2006,
c     ApJ, volume 652, pp.1028.  Please reference this paper if you
c     publish results dereived from this model.
c
c     Developed by Laura Brenneman and Chris Reynolds
c     Dept. of Astronomy, University of Maryland, College Park
c
c
c     Modified Mar 09 for xspec12: calls fgmodf instead of using
c        xspec11's LMODDIR setting.  (CG)
c
c     
      IMPLICIT NONE
      INTEGER IFL, NE
      REAL*4 EAR(0:NE), PARAM(9), PHOTAR(NE), PHOTER(NE)

c
c---------Initialization--------------
c
      character*255 datafile
      character*255 fgmodf
      integer lenact
      integer energy(10000)
      integer i,j,k,ii,jj,ii1,ii2,igstar2,irad,ilgrad,ilgrad_max
      integer nradii,ng,readflag,ia,imu0,abins,mu0bins
      parameter(nradii=50,ng=20,abins=20,mu0bins=20)
      real*8 rplus,sumspec,lgrad
      real*8 a,theta0,mu0,gstar(ng),g(nradii,ng)
      real*8 a_tab(abins),mu0_tab(mu0bins),temp
      real*8 aintfac,mu0intfac
      real*8 trff_tab(nradii,ng,2,abins,mu0bins)
      real*8 cosne_tab(nradii,ng,2,abins,mu0bins)
      real*8 gmin_tab(nradii,abins,mu0bins)
      real*8 gmax_tab(nradii,abins,mu0bins)
      real*8 trffh(2),re(nradii),gmin(nradii)
      real*8 gmax(nradii),trff(nradii,ng,2),cosne(nradii,ng,2)
      real*8 ispec,lspec(ne),normspec(ne),eeo(0:ne)
      real*8 lspecfine(4*ne),eeofine(4*ne)
      real*8 intgmin,intgmax,inttf(ng,2),intmu(ng,2),rad,intfac
      real*8 intgs
      real*8 rms,marginal,gee,gstar2,trf,mu,eem1,eem2
      real*8 rmin,rmax,alp1,alp2,rbreak,lineE,z
      real*8 rmin_grid,rmax_grid
      real*8 pi,r1,r2,r(nradii),wr(nradii),g1,g2,wg(ng)
      external spingauleg, fgmodf
      save a_tab,trff_tab,cosne_tab,gmin_tab,gmax_tab,mu0_tab
c
c     Get parameters of model from xspec
c
      lineE=param(1)
      alp1=param(2)
      alp2=param(3)
      rbreak=param(4)
      a=param(5)
      theta0=param(6)
      rmin=param(7)
      rmax=param(8)
      z=param(9)

      rms=marginal(a)

      rmin=rmin*rms
      rmax=rmax*rms

c
c     Set pi and convert angle variables
c
      pi = 4.d0*atan(1.d0)
      theta0=theta0*pi/180.d0
      mu0=cos(theta0)
c
c     Read in transfer function from the kerrtable.dat file.  We use
c     the value of a_tab(1) to determine whether this is the first call
c     to the code in a given xspec session.  On subsequent calls, we
c     do not need to read the file.
c
      if ( dabs(a_tab(1) - 1.0e-2) .gt. 1.0e-3) then
      datafile = fgmodf()
      datafile=datafile(1:lenact(datafile))//'/kerrtable.dat'
      open (8,file=datafile,status='unknown')
      read(8,*) ii,jj
      read(8,*) a_tab
      read(8,*) mu0_tab
      do i=1,abins
         do j=1,mu0bins
            read(8,*) gmin
            read(8,*) gmax
            read(8,*) trff
            read(8,*) cosne
            do ii=1,nradii
               gmin_tab(ii,i,j)=gmin(ii)
               gmax_tab(ii,i,j)=gmax(ii)
               do jj=1,ng
                  do k=1,2
                     trff_tab(ii,jj,k,i,j)=trff(ii,jj,k)
                     cosne_tab(ii,jj,k,i,j)=cosne(ii,jj,k)
                  enddo
               enddo
            enddo
         enddo
      enddo
      close(8)
      endif
c
c     Work out interpolation factors in spin and mu0 directions
c
      ia=1
      do i=1,abins
         if (a_tab(i) .lt. a) ia=i
      enddo
      aintfac=(a-a_tab(ia))/(a_tab(ia+1)-a_tab(ia))
c
      imu0=1
      do i=1,mu0bins
         if (mu0_tab(i) .lt. mu0) imu0=i
      enddo
      mu0intfac=(mu0-mu0_tab(imu0))/(mu0_tab(imu0+1)-mu0_tab(imu0))
c     
c     attempt to fix low-inclination problem
c
      if (mu0_tab(mu0bins) .lt. mu0) then
         imu0=mu0bins-1
         mu0intfac=1
      endif
c
c     Interpolate transfer function in a and mu0 plane
c
      do i=1,nradii
         gmin(i)=(1.0-aintfac)*(1.0-mu0intfac) * 
     $        gmin_tab(i,ia,imu0)
     $        +aintfac*(1.0-mu0intfac)*
     $        gmin_tab(i,ia+1,imu0)
     $        +(1.0-aintfac)*mu0intfac*
     $        gmin_tab(i,ia,imu0+1)
     $        +aintfac*mu0intfac*
     $        gmin_tab(i,ia+1,imu0+1)
c
         gmax(i)=(1.0-aintfac)*(1.0-mu0intfac) * 
     $        gmax_tab(i,ia,imu0)
     $        +aintfac*(1.0-mu0intfac)*
     $        gmax_tab(i,ia+1,imu0)
     $        +(1.0-aintfac)*mu0intfac*
     $        gmax_tab(i,ia,imu0+1)
     $        +aintfac*mu0intfac*
     $        gmax_tab(i,ia+1,imu0+1)
         do j=1,ng
            do k=1,2
               trff(i,j,k)=(1.0-aintfac)*(1.0-mu0intfac) * 
     $              trff_tab(i,j,k,ia,imu0)
     $              +aintfac*(1.0-mu0intfac)*
     $              trff_tab(i,j,k,ia+1,imu0)
     $              +(1.0-aintfac)*mu0intfac*
     $              trff_tab(i,j,k,ia,imu0+1)
     $              +aintfac*mu0intfac*
     $              trff_tab(i,j,k,ia+1,imu0+1)
c
               cosne(i,j,k)=(1.0-aintfac)*(1.0-mu0intfac) * 
     $              cosne_tab(i,j,k,ia,imu0)
     $              +aintfac*(1.0-mu0intfac)*
     $              cosne_tab(i,j,k,ia+1,imu0)
     $              +(1.0-aintfac)*mu0intfac*
     $              cosne_tab(i,j,k,ia,imu0+1)
     $              +aintfac*mu0intfac*
     $              cosne_tab(i,j,k,ia+1,imu0+1)
            enddo
         enddo
      enddo
c
c     Set up a radial grid: inversely spaced, and define integration values
c     for g* via spingauleg as in radial case 
c
      rmin_grid=rms
      rmax_grid=500*rms
      r1=1.d0/sqrt(rmax_grid)
      r2=1.d0/sqrt(rmin_grid)
      call spingauleg(r1,r2,r,wr,nradii)
      do i=1,nradii
         re(i)=1.d0/(r(i)**2.0)
      enddo
      g1=0.d0
      g2=1.d0
      call spingauleg(g1,g2,gstar,wg,ng)
c
c---------Integrate the line profile-------------------------
c

c
c     Generate energy grid and finer grid within it (4 x finer) 
c     to effectively get greater resolution than before without 
c      smoothing.  We will linearly interpolate between grid point 
c     values later 
c
      eeo(0)=dble(ear(0))
      do ii=1,ne
         lspec(ii)=0.d0
         eeo(ii)=dble(ear(ii))
      enddo
      do ii=1,ne
         do j=1,4
            lspecfine((ii-1)*4+j)=0.d0
            intfac=float(j)/4.0
            eeofine((ii-1)*4+j) = 
     &           intfac*eeo(ii)+(1.0-intfac)*eeo(ii-1)
         enddo
      enddo
c
c     Latest generation of line integration.  Integrates over a
c     large number of radii, using linear radial interpolation 
c     of the TF as well as gmin and gmax 
c
      irad=nradii-1
      do ii=1,nradii
         if (rmin .lt. re(ii)) irad=ii
      enddo
c
c     work out the number of zones for the final radial integral
c     ... currently set to ne per decade of radius
c
      ilgrad_max=3*int(float(ne)*dlog10(re(1)/re(nradii)))
c
c     Perform radial integral...
c
      do ilgrad=1,ilgrad_max
c
c        work out corresponding rrelevant 
c
         rad=rmin * 10 ** ( dfloat(ilgrad-1)*dlog(re(1)/re(nradii))/
     &     dfloat(ilgrad_max-1))

         if ((rad .gt. rmin) .and. (rad .lt. rmax)) then

         if (rad .gt. re(irad)) irad=irad-1
         intfac=(rad-re(irad+1))/(re(irad)-re(irad+1))
         do j=1,ng
            do k=1,2
               inttf(j,k)=intfac*trff(irad,j,k) + 
     &              (1.0-intfac)*trff(irad+1,j,k)
               intmu(j,k)=intfac*cosne(irad,j,k) + 
     &              (1.0-intfac)*cosne(irad+1,j,k)
            enddo
         enddo
         intgmin=intfac*gmin(irad)+(1.0-intfac)*gmin(irad+1)
         intgmax=intfac*gmax(irad)+(1.0-intfac)*gmax(irad+1)
c
c     EMISSIVITY PROFILE is hardwired in here.  Currently a broken
c     powerlaw
c
         if (rad .lt. rbreak) then
            ispec=(rad/rbreak)**(-alp1)
         else 
            ispec=(rad/rbreak)**(-alp2)
         endif
c
         eem1=lineE*intgmin/(1.0d0+z)
         eem2=lineE*intgmax/(1.0d0+z)
         ii1=1
         ii2=1

         do ii=1,4*ne
            if (eeofine(ii) .lt. eem1) ii1=ii
         enddo
         do ii=ii1-1,4*ne
            if (eeofine(ii) .lt. eem2) ii2=ii
         enddo

         if ((ii1 .gt. 1) .and. (ii2 .gt. 1)) then
 2005       do ii=ii1+1,ii2
            gee=eeofine(ii)/(lineE/(1.0d0+z))
            gstar2=(gee-intgmin)/(intgmax-intgmin)
            do k=1,2
               if (gstar2 .le. gstar(1)) then
                  trf=inttf(1,k)
                  mu=intmu(1,k)
               endif
               if (gstar2 .ge. gstar(ng)) then
                  trf=inttf(ng,k)
                  mu=intmu(ng,k)
               endif
               if ((gstar2 .lt. gstar(ng)) .and. 
     &              (gstar2 .gt. gstar(1))) then
                  do j=1,ng-1
                     if (gstar(j) .lt. gstar2) igstar2=j
                  enddo
                  intgs=(gstar2-gstar(igstar2))/
     &                 (gstar(igstar2+1)-gstar(igstar2))
                  trf=intgs*inttf(igstar2,k) + 
     &                 (1.0-intgs)*inttf(igstar2+1,k)
                  mu=intgs*intmu(igstar2,k) + 
     &                 (1.0-intgs)*intmu(igstar2+1,k)
               endif
c
c     Next line actually is the guts of the integral. The
c     LIMB DARKENING is hardwired in here.  It is currently set to
c     mu*(1+2.06*mu)
c
               lspecfine(ii)=lspecfine(ii) +
     &              rad*gee*(2.0*pi*gee)**2.0
     &              *trf*ispec*(1.0+2.06*mu)*rad/
     &              (sqrt(gstar2-gstar2**2.0)*(intgmax-intgmin))
            enddo
         enddo
         endif
         endif
      enddo
c
c     Bin up lspecfine to give to lspec 
c
      do i=1,ne
         do j=1,4
            lspec(i)=lspec(i)+lspecfine((i-1)*4+j)
         enddo
      enddo

c
c---------Output of spectral luminosity-----------------
c
c     Divide each lspec(ii) by the observed energy of that 
c     gridpoint --> ph/cm^2/s units 
c
      do ii=1,ne
         lspec(ii)=lspec(ii)/eeo(ii)
      enddo
      sumspec=0.d0
      do ii=1,ne
         if (lspec(ii) .gt. 0.d0) then
c
c     Sumspec weighted by the energy bin size 
c
            sumspec=sumspec+lspec(ii)*(eeo(ii)-eeo(ii-1))
         endif
      enddo
c      open(7,file='lspec.dat',status='unknown')
      do ii=1,ne
         if (lspec(ii) .gt. 0.d0) then
c
c     Normspec weighted by the energy bin size 
c
            normspec(ii)=lspec(ii)*(eeo(ii)-eeo(ii-1))/sumspec
         else 
            normspec(ii)=0.d0
         endif
         photar(ii)=normspec(ii)
c         write(7,*) eeo(ii),normspec(ii)
      enddo
c      close(unit=7)

      end
c
c---------Subroutines------------------------
c
      function marginal(a)
      implicit none
      real*8 marginal,a,Z1,Z2
      Z1=1.0+(1.0-a**2.0)**0.33*((1.0+a)**0.33+(1.0-a)**0.33)
      Z2=((3.0*a**2.0)+(Z1**2.0))**0.5
      marginal=3.0+Z2-((3.0-Z1)*(3.0+Z1+(2*Z2)))**0.5
      return

      end 


C========================================================================

C  Subroutine to calculate the abscissas and weights for the Gauss-Legendre-
C  Quadrature. The routine is based on the NUMERICAL RECIPES and uses an
C  algorithem of G.B. Rybicki.
C  Input: x1 ,x2: range of integration.
C          n: order of the orthogonal polynomials and the quadrature formula.
C  Output: x = x(n): array of the abscissas.
C          w = w(n): array of the weights.


      SUBROUTINE SPINGAULEG(X1,X2,X,W,N)

      INTEGER N,M,I,J
      REAL*8 X1,X2,X(N),W(N)
      REAL*8 PI,XM,XL,Z,P1,P2,P3,Z1,PP,EPS
      PARAMETER (PI = 3.14159265358979323846D0)
      PARAMETER (EPS=3.D-14)

      M=(N+1)/2
      XM=0.5D0*(X2+X1)
      XL=0.5D0*(X2-X1)
      DO 12 I=1,M
         Z=COS(PI*(I-.25D0)/(N+.5D0))
 1       CONTINUE
         P1=1.D0
         P2=0.D0
         DO 11 J=1,N
            P3=P2
            P2=P1
            P1=((2.D0*J-1.D0)*Z*P2-(J-1.D0)*P3)/J
 11      CONTINUE
         PP=N*(Z*P1-P2)/(Z*Z-1.D0)
         Z1=Z
         Z=Z1-P1/PP
         IF(ABS(Z-Z1).GT.EPS) GOTO 1
         X(I)=XM-XL*Z
         X(N+1-I)=XM+XL*Z
         W(I)=2.D0*XL/((1.D0-Z*Z)*PP*PP)
         W(N+1-I)=W(I)
 12   CONTINUE
      RETURN
      END



