subroutine spin(EAR,NE,PARAM,IFL,PHOTAR,PHOTER) IMPLICIT NONE INTEGER IFL, NE REAL*4 EAR(0:NE), PARAM(9), PHOTAR(NE), PHOTER(NE) c program kerrdisk c c Set up for installation into XSPEC v.11.2.3. Eliminates need c for delta function in flux integral, DOES NOT include radiation c within ISCO, has Laor limb darkening. Ionization, optical depth c clauses for radiation to be added later. Enables use of a c broken power law to describe the continuum. c c Most recent innovation: linearly interpolates calculation of c transfer function from a grid pre-computed with Speith algorithms c in makekerrtable.f --> kerrtable.dat. Grid in this file is over c a given range of r from rms to 500*rms. Range in a and mu0 c (cos(theta0)) from a 10 x 10 grid from 0.01-0.998, 8-89 degrees. c This grid will be enlarged for the final version of this code. c c created by LB, 2/14/06 c renamed 8/17/06 c c---------Initialization-------------- c real*8 rplus,sumspec integer nradii,ng,readflag,ia,imu0,abins,mu0bins parameter(nradii=50,ng=20,abins=20,mu0bins=20) real*8 a,theta0,mu0,gstar(ng),g(nradii,ng),cosneh(2) 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(10*ne),eeofine(10*ne) real*8 intgmin,intgmax,inttf(ng,2),intmu(ng,2),rad,lgrad,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) integer energy(10000) integer i,j,k,ii,jj,ii1,ii2,igstar2,irad external gauleg cc save readflag save gmin_tab,gmax_tab save trff_tab,cosne_tab cc DATA readflag / 1 / c c---------Read parameters of model---------- 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) if (rmin .lt. rms) then rmin=rms endif c Set constants: pi = 4.d0*atan(1.d0) c Input of observation angle theta0: theta0=theta0*pi/180.d0 mu0=cos(theta0) c c Read in from the kerrtable.dat file c cc if (readflag .ne. 0) then cc readflag=0 c c Open file c c c USER MUST SPECIFY FULL PATH TO KERRTABLE.DAT HERE -- FILE c SHOULD BE IN THE LOCAL MODELS DIRECTORY ($LMODDIR) c open (8,file= & '/n/artemis/lwb/software/xspec_local_models/kerrtable.dat' & ,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) cc 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)) 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 Interpolate transfer function c do i=1,nradii temp=(1.0-aintfac)*(1.0-mu0intfac) * $ gmin_tab(i,ia,imu0) temp=temp+aintfac*(1.0-mu0intfac)* $ gmin_tab(i,ia+1,imu0) temp=temp+(1.0-aintfac)*mu0intfac* $ gmin_tab(i,ia,imu0+1) temp=temp+aintfac*mu0intfac* $ gmin_tab(i,ia+1,imu0+1) gmin(i)=temp temp=(1.0-aintfac)*(1.0-mu0intfac) * $ gmax_tab(i,ia,imu0) temp=temp+aintfac*(1.0-mu0intfac)* $ gmax_tab(i,ia+1,imu0) temp=temp+(1.0-aintfac)*mu0intfac* $ gmax_tab(i,ia,imu0+1) temp=temp+aintfac*mu0intfac* $ gmax_tab(i,ia+1,imu0+1) gmax(i)=temp do j=1,ng do k=1,2 temp=(1.0-aintfac)*(1.0-mu0intfac) * $ trff_tab(i,j,k,ia,imu0) temp=temp+aintfac*(1.0-mu0intfac)* $ trff_tab(i,j,k,ia+1,imu0) temp=temp+(1.0-aintfac)*mu0intfac* $ trff_tab(i,j,k,ia,imu0+1) temp=temp+aintfac*mu0intfac* $ trff_tab(i,j,k,ia+1,imu0+1) trff(i,j,k)=temp temp=(1.0-aintfac)*(1.0-mu0intfac) * $ cosne_tab(i,j,k,ia,imu0) temp=temp+aintfac*(1.0-mu0intfac)* $ cosne_tab(i,j,k,ia+1,imu0) temp=temp+(1.0-aintfac)*mu0intfac* $ cosne_tab(i,j,k,ia,imu0+1) temp=temp+aintfac*mu0intfac* $ cosne_tab(i,j,k,ia+1,imu0+1) cosne(i,j,k)=temp enddo enddo enddo c c Set up a radial grid: inversely spaced, and define integration values c for g* via gauleg 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 gauleg(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 gauleg(g1,g2,gstar,wg,ng) c c---------Integrate the line profile------------------------- c c c Generate energy grid and finer grid within it (10 x finer) to effectively get c greater resolution than before without smoothing. We will linearly interpolate c between grid point values later c do ii=0,ne lspec(ii)=0.d0 eeo(ii)=dble(ear(ii)) enddo do ii=1,ne do j=1,10 lspecfine((ii-1)*10+j)=0.d0 intfac=float(j)/10.0 eeofine((ii-1)*10+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 interpolation of the TF as well c as gmin and gmax c irad=nradii-1 do lgrad=dlog10(re(nradii)),dlog10(re(1)),0.2/float(ne) rad=10.0**(lgrad) 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) if (rad .lt. rbreak) then ispec=(rad/rbreak)**(-alp1) else ispec=(rad/rbreak)**(-alp2) endif eem1=lineE*intgmin/(1.0d0+z) eem2=lineE*intgmax/(1.0d0+z) ii1=1 ii2=1 do ii=1,10*ne if (eeofine(ii) .lt. eem1) ii1=ii 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 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 Write out lspecfine to a file for checking c open(7,file='lspecfine.dat',status='unknown') do i=1,ne*10 write(7,*) eeofine(i),lspecfine(i) enddo close(7) c c Bin up lspecfine to give to lspec c do i=1,ne do j=1,10 lspec(i)=lspec(i)+lspecfine((i-1)*10+j) enddo enddo c c---------Output of spectral luminosity----------------- c c Divide each lspec(ii) by the observed energy of that 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 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) write(7,*) eeo(ii),normspec(ii) enddo 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 GAULEG(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