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 IMPLICIT NONE INTEGER IFL, NE REAL*4 EAR(0:NE), PARAM(9), PHOTAR(NE), PHOTER(NE) c c---------Initialization-------------- c character*255 datafile 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 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 call getenv('LMODDIR', datafile) 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 do ii=0,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