subroutine xskdline (ccheck,ear,ne,param,ifl,photar,photer) implicit none integer ni,ngp,nc,ng,nr,np,nb parameter (ni=2,ngp=5,nc=2,ng=100,nr=100,np=5,nb=100) integer ne, ifl real ear(0:ne), param(9), photar(ne), photer(ne) logical ccheck c This model calculates the line profile from a Keplerian Accretion disk around c both Schwarzschild & Kerr (arbitrary spin) black holes. c c It requires the use of an extensive set of transfer functions, which must be c downloaded from the XSPEC website. These transfer functions should be placed c in $LHEASOFT/../spectral/xspec/manager/XSKDLINE/ c c It is not called directly from XSPEC, rather, it is called by the files c xskdlN (additive model), xskdcN (convolution model, via xskdconv). Here, c N is either 0 or 1 and denotes different modes of operation for the model. c For further information how each of the modes works, consult xskdlN.f, c xskdcN.f, xskdconv.f as appropriate. c c Finally, the user should note that chatter>=12 will case the model to dump c the current line parameters to the screen double precision rval double precision r1,r2 double precision a,as,obs double precision risco,rms double precision ec,emin,emax double precision rmin,rmax,index integer rtype,mtype,ftype,fts,nes character*80 buff logical qfirst double precision de(0:nb-1) double precision e(0:nb-1),n(nb) logical nmask(0:nb-1) logical check1,check2,check4 integer loc,i save nmask,n,e,de,r1,r2 save ec,a,as,obs,rmin,rmax,index save rtype,mtype,ftype,fts,nes save qfirst external rms DATA qfirst/.TRUE./ ftype = int(param(9)) If (qfirst) then nes = ne fts = ftype as = dble(param(2)) If (ftype.EQ.0) then emin = dble(ear(0)) emax = dble(ear(ne-1)) ec = dble(param(1)) a = dble(param(2)) rmin = dble(param(3)) rmax = dble(param(4)) obs = dble(param(5)) index = dble(param(6)) rtype = int(param(7)) mtype = int(param(8)) Else emin = dble(ear(0)) emax = dble(ear(ne-1)) ec = dble(param(1)) a = dble(param(2)) rmin = dble(param(3)) rmax = dble(param(4)) obs = dble(param(5)) index = dble(param(6)) rtype = int(param(7)) mtype = int(param(8)) End If write(*,*)'This is the XSKDLine model by Kris Beckwith' End If check1=.FALSE. If (qfirst) check1=.TRUE. If (ftype.NE.fts) check1=.TRUE. If (ftype.EQ.0) then If (rmin.NE.dble(param(3))) check1=.TRUE. If (rmax.NE.dble(param(4))) check1=.TRUE. If (obs.NE.dble(param(5))) check1=.TRUE. If (index.NE.dble(param(6))) check1=.TRUE. If (rtype.NE.int(param(7))) check1=.TRUE. If (mtype.NE.int(param(8))) check1=.TRUE. Else if ((ftype.EQ.1).OR.(ftype.EQ.2)) then If (a.NE.dble(param(2))) check1=.TRUE. If (rmin.NE.dble(param(3))) check1=.TRUE. If (rmax.NE.dble(param(4))) check1=.TRUE. If (obs.NE.dble(param(5))) check1=.TRUE. If (index.NE.dble(param(6))) check1=.TRUE. If (rtype.NE.int(param(7))) check1=.TRUE. If (mtype.NE.int(param(8))) check1=.TRUE. Else if (ftype.EQ.3) then If (a.NE.dble(param(2))) check1=.TRUE. If (rmax.NE.dble(param(4))) check1=.TRUE. If (obs.NE.dble(param(5))) check1=.TRUE. If (index.NE.dble(param(6))) check1=.TRUE. If (rtype.NE.int(param(7))) check1=.TRUE. If (mtype.NE.int(param(8))) check1=.TRUE. If (rmin.NE.dble(param(3))) then risco = rms(dble(param(2))) If ((rmin.NE.dble(param(3))).AND. & (dble(param(3)).GE.risco)) check1=.TRUE. If ((rmin.GT.risco).AND. & (dble(param(3)).LT.risco)) check1=.TRUE. End If End If !write(*,*)'Does transfer function require calculation?',check1 check2=.NOT.(check1) If (check1) then a = dble(param(2)) rmin = dble(param(3)) rmax = dble(param(4)) obs = dble(param(5)) index = dble(param(6)) rtype = int(param(7)) mtype = int(param(8)) ftype = int(param(9)) !write(*,*)'Line parameters',rmin,rmax,obs,index,rtype,mtype check4=.FALSE. rval = rms(0.998d0) If ((ftype.EQ.0).AND.(rtype.EQ.2) & .AND.(rmin.LT.rval)) check4=.TRUE. If ((ftype.EQ.1).AND.(rtype.EQ.2) & .AND.(a.GT.0.998d0)) check4=.TRUE. If (check4) then rtype=0 write(*,*)'WARNING - Zero Torque Inner Boundary condition &is invalid for a = 1.0, resetting radial emissivity type' End If call lineinit (nmask,n,e,de,nb) !write(*,*)'Initialised line' check2=.FALSE. call obsdet (check2,nmask,n,e,de,r1,r2,rmin,rmax,a,as,obs, & index,rtype,mtype,ftype,nb,ni,ngp,ng,nr,nc,np) !write(*,*)'Calculated line for inclination, result',check2 param(2) = real(a) param(3) = real(rmin) param(4) = real(rmax) param(5) = real(obs) param(6) = real(index) param(7) = real(rtype) param(8) = real(mtype) !write(*,*)'Passed back parameters' End If check4=.FALSE. If (check2) check4=.TRUE. If (nes.NE.ne) check4=.TRUE. If (ec.NE.dble(param(1))) check4=.TRUE. If (emin.NE.dble(ear(0))) check4=.TRUE. If (emax.NE.dble(ear(ne-1))) check4=.TRUE. !write(*,*)'Does line require rebinning?',check4 If (check4) then nes = ne emin = dble(ear(0)) emax = dble(ear(ne-1)) ec = dble(param(1)) !write(*,*)'Line parameters',nes,emin,emax,ec qfirst=.FALSE. call rebinp2 (ne,ec,ear,photar,nmask,n,e,de,nb) !write(*,*)'Rebinned line onto XSPEC grid' End If as = a fts=ftype If (ccheck) then write(buff,'(a,f6.2)')'Spin=',a call xwrite(buff,12) write(buff,'(a,f6.2)')'Rmin(G)=',r1 call xwrite(buff,12) write(buff,'(a,f6.2)')'Rmax(G)=',r2 call xwrite(buff,12) End If return end c ************************************************************************* c * These are the subroutines for xskdline * c ************************************************************************* subroutine lineinit (nmask,n,e,de,ne) implicit none integer ne double precision n(ne) double precision e(0:ne-1) double precision de(0:ne-1) logical nmask(0:ne-1) integer i do i=1,ne-1 n(i-1)=0.0d0 e(i-1)=0.0d0 de(i-1)=0.0d0 nmask(i-1)=.FALSE. end do return end subroutine obsdet (checkout,nmask,n,e,de,r1,r2,rmin,rmax,a,as,obs, & index,rtype,mtype,ftype,nb,ni,ngp,ng,nr,nc,np) implicit none integer nb,ni,ngp integer ng,nr,nc,np integer rtype,mtype,ftype double precision r1,r2 double precision rmin,rmax double precision de(0:nb-1) double precision a,as,obs,index double precision e(0:nb-1),n(nb) logical checkout,nmask(0:nb-1) integer no,nb1 parameter (no=2,nb1=100) double precision oval(no) double precision gb(no,2) double precision g(no,0:nb1-1) double precision dg(no,0:nb1-1) double precision nphot(no,nb1) logical omask(no),check1 logical gmask(no,0:nb1-1) call ovaldet (omask,oval,obs,no) !write(*,*)'Inclinations to use',oval(1),oval(2) check1=.FALSE. call olinedet (check1,gmask,nphot,g,dg,omask,oval,r1,r2,rmin,rmax, & a,as,index,rtype,mtype,ftype,no,nb1,ni,ngp,ng,nr,nc,np) !write(*,*)'Determined lines for inclinations, result',check1 checkout=.FALSE. If (check1) then call olmergep1 (omask,gmask,gb,nphot,g,dg,no,nb1) !write(*,*)'Placed lines onto common grid' call olmergep2 (checkout,omask,oval,obs,nmask,n,e,de,gmask, & gb,nphot,g,dg,nb,no,nb1) !write(*,*)'Interpolated between lines, result',checkout End If return end c ************************************************************************* c * These are the subroutines for obsdet * c ************************************************************************* subroutine ovaldet (omask,oval,obs,no) implicit none integer no double precision oval(no),obs logical omask(no) double precision o o=dble(nint(obs)) !write(*,*)'Integer form of inclination',o,obs If (o.EQ.obs) then oval(1)=obs oval(2)=0.0d0 omask(1)=.TRUE. omask(2)=.FALSE. Else if (o.LT.obs) then oval(1)=o oval(2)=o+1.0d0 omask(1)=.TRUE. omask(2)=.TRUE. Else if (o.GT.obs) then oval(2)=o oval(1)=o-1.0d0 omask(1)=.TRUE. omask(2)=.TRUE. End If return end subroutine olinedet (checkout,nmask,n,e,de,omask,oval,r1,r2, & rmin,rmax,a,as,index,rtype,mtype,ftype,no,nb,ni,ngp,ng,nr,nc,np) implicit none integer ng,nr,nc,np integer no,nb,ni,ngp integer rtype,mtype,ftype double precision r1,r2 double precision rmin,rmax double precision a,as,index logical checkout double precision oval(no) double precision n(no,nb) double precision e(no,0:nb-1) double precision de(no,0:nb-1) logical omask(no) logical nmask(no,0:nb-1) integer nb1 parameter (nb1=100) double precision obs double precision g(0:nb1-1) double precision dg(0:nb1-1) double precision nphot(nb1) logical gmask(0:nb1-1) logical check1,check2 integer n1,n2 integer i,j n1=0 n2=0 do i=1,no check1=omask(i) If (check1) then n1=n1+1 obs=oval(i) check2=.FALSE. If (ftype.EQ.0) call lineinit0 (check2,gmask,nphot,g,dg, & r1,r2,rmin,rmax,a,obs,index,rtype,mtype,nb,ni,ngp,ng,nr,nc,np) If (ftype.EQ.1) call lineinit1 (check2,gmask,nphot,g,dg, &r1,r2,rmin,rmax,a,as,obs,index,rtype,mtype,nb,ni,ngp,ng,nr,nc,np) If (ftype.EQ.2) call lineinit2 (check2,gmask,nphot,g,dg, &r1,r2,rmin,rmax,a,as,obs,index,rtype,mtype,nb,ni,ngp,ng,nr,nc,np) !write(*,*)i,' Calculated line for fit type, result',check1 If (check2) then n2=n2+1 omask(i)=.TRUE. do j=1,nb n(i,j)=nphot(j) e(i,j-1)=g(j-1) de(i,j-1)=dg(j-1) nmask(i,j-1)=gmask(j-1) end do !write(*,*)i,' Stored line profiles' End If End If end do checkout=.FALSE. If (n1.EQ.n2) checkout=.TRUE. !write(*,*)'Were valid inclinations found?',checkout,n1,n2 return end subroutine olmergep1 (omask,nmask,gb,n,e,de,no,nb) implicit none integer no,nb double precision gb(no,2) double precision n(no,nb) double precision e(no,0:nb-1) double precision de(no,0:nb-1) logical omask(no) logical nmask(no,0:nb-1) integer nb1 parameter (nb1=100) double precision gmin,gmax,x,y double precision nphot(nb1),dg1 double precision g(nb1),dg(nb1) logical gmask(nb1),check1 integer i,j If (omask(1).AND.omask(2)) then gmin=1.0d2 gmax=-1.0d2 do i=1,no gb(i,1)=e(i,0) gb(i,2)=e(i,nb-1) gmin=min(gmin,e(i,0)) gmax=max(gmax,e(i,nb-1)) end do !write(*,*)'Limits on energy grid',gmin,gmax dg1=(gmax-gmin)/dble(nb-1) !write(*,*)'Bin width',dg1 do i=1,no do j=1,nb dg(j)=de(i,j-1) gmask(j)=nmask(i,j-1) nphot(j)=n(i,j)/de(i,j-1) g(j)=e(i,j-1)+0.5d0*de(i,j-1) end do !write(*,*)i,' Setup interpolation arrays' e(i,0)=gmin do j=1,nb n(i,j)=0.0d0 de(i,j-1)=dg1 nmask(i,j-1)=.FALSE. If (j.GT.1) e(i,j-1)=e(i,j-2)+dg1 !write(*,*)i,j,' Initialised global grid' y=0.0d0 check1=.FALSE. x=e(i,j-1)+0.5d0*dg1 call linterpolate (check1,y,x,gmask,g,nphot,nb1) !write(*,*)i,j,' Interpolation',x,y,check1 If (check1) then n(i,j)=y*dg1 nmask(i,j-1)=.TRUE. End If !write(50,*)i,j-1,e(i,j-1),de(i,j-1),n(i,j) end do !write(50,*) end do End If return end subroutine olmergep2 (checkout,omask,oval,obs,nmask,n,e,de,gmask, & gb,nphot,g,dg,nb,ni,ng) implicit none integer nb,ni,ng double precision gb(ni,2) double precision oval(ni),obs double precision nphot(ni,ng) double precision de(0:nb-1) double precision g(ni,0:ng-1) double precision dg(ni,0:ng-1) double precision e(0:nb-1),n(nb) logical omask(ni),gmask(ni,0:nb-1) logical checkout,nmask(0:nb-1) integer nb1,ni1,ng1 parameter (nb1=100,ni1=2,ng1=100) double precision work1a(nb1) double precision work1d(nb1) double precision work1e(nb1) double precision work1f(nb1) double precision gmin,gmax,de1 double precision work1b(0:nb1-1) double precision work1c(0:nb1-1) double precision rfrac,risco,rms,x,y double precision aval(ni1),rval(ni1) logical mask1a(ni1) logical mask1b(0:nb1-1) logical mask1c(nb1),check1,check2 integer i,j,nvalid rfrac=1.0d0 If (omask(1).AND.omask(2)) rfrac=(obs-oval(1))/(oval(2)-oval(1)) !write(*,*)'Fractional distance between inclinations',rfrac gmin=g(1,0) gmax=g(1,ng-1) de1=(gmax-gmin)/dble(nb1-1) !write(*,*)'Lower, upper limits of line',gmin,gmax !write(*,*)'Bin width',de1 nvalid=0 work1b(0)=gmin do i=1,nb1 work1a(i)=0.0d0 work1c(i-1)=de1 mask1b(i-1)=.FALSE. If (i.GT.1) work1b(i-1)=work1b(i-2)+de1 If (omask(1).AND.omask(2)) then If (gmask(1,i-1).AND.gmask(2,i-1)) then nvalid=nvalid+1 mask1b(i-1)=.TRUE. work1a(i)=min(nphot(1,i),nphot(2,i))+ & dabs(rfrac*(nphot(2,i)-nphot(1,i))) Else If (gmask(1,i-1).AND.((work1b(i-1).LT.gb(2,1)).OR. & (work1b(i-1).GT.gb(2,2)))) then nvalid=nvalid+1 mask1b(i-1)=.TRUE. work1a(i)=rfrac*nphot(1,i) Else if (gmask(2,i-1).AND.((work1b(i-1).LT.gb(1,1)).OR. & (work1b(i-1).GT.gb(1,2)))) then nvalid=nvalid+1 mask1b(i-1)=.TRUE. work1a(i)=rfrac*nphot(2,i) End If End If !write(48,*)i,work1b(i-1),nphot(1,i),work1a(i),nphot(2,i), ! & gmask(1,i-1),mask1b(i-1),gmask(2,i-1) Else if (omask(1)) then If (gmask(1,i-1)) then nvalid=nvalid+1 work1a(i)=nphot(1,i) mask1b(i-1)=.TRUE. End If Else if (omask(2)) then If (gmask(2,i-1)) then nvalid=nvalid+1 work1a(i)=nphot(2,i) mask1b(i-1)=.TRUE. End If End If end do check1=.FALSE. If (nvalid.GT.0) check1=.TRUE. !write(*,*)'Were valid points found in line?',check1 nvalid=0 If (check1) then gmin=work1b(0) gmax=work1b(nb1-1) If (omask(1).AND.omask(2)) then gmin=min(gb(1,1),gb(2,1))+dabs(rfrac*(gb(2,1)-gb(1,1))) gmax=max(gb(1,1),gb(2,2))+dabs(rfrac*(gb(2,2)-gb(1,2))) do i=1,nb1 mask1c(i)=mask1b(i-1) work1f(i)=work1c(i-1) work1d(i)=work1a(i)/work1c(i-1) work1e(i)=work1b(i-1)+0.5d0*work1c(i-1) !write(49,*)work1e(i),work1d(i)*work1c(i-1) end do !write(49,*)'no no' !write(*,*)'Setup interpolation arrays' End If !write(*,*)'Lower, upper limits on line',gmin,gmax de1=(gmax-gmin)/dble(nb1-1) !write(*,*)'Bin width',de1 e(0)=gmin do i=1,nb n(i)=0.0d0 de(i-1)=de1 nmask(i-1)=.FALSE. If (i.GT.1) e(i-1)=e(i-2)+de1 If (omask(1).AND.omask(2)) then y=0.0d0 check2=.FALSE. x=e(i-1)+0.5d0*de1 call linterpolate (check2,y,x,mask1c,work1e,work1d,nb1) !write(*,*)i,' Line interpolation',x,y,check2 If (check2) then n(i)=y*de1 nvalid=nvalid+1 nmask(i-1)=.TRUE. End If Else If (mask1b(i-1)) then n(i)=work1a(i) nvalid=nvalid+1 nmask(i-1)=.TRUE. End If End If !write(49,*)e(i-1),n(i) end do !write(49,*)'no no' End If checkout=.FALSE. If (nvalid.GT.0) checkout=.TRUE. !write(*,*)'Were valid points found in line?',checkout return end subroutine lineinit0 (checkout,nmask,n,e,de,r1,r2,rmin,rmax, & a,obs,index,rtype,mtype,nb,ni,ngp,ng,nr,nc,np) implicit none integer ngp,ng,nr,nc,np integer rtype,mtype,nb,ni double precision r1,r2 double precision rmin,rmax double precision de(0:nb-1) double precision a,obs,index double precision e(0:nb-1),n(nb) logical checkout,nmask(0:nb-1) integer nb1,ni1,ngp1,ng1,nr1,nc1,np1 parameter (nb1=100,ni1=2,ngp1=5,nc1=2,ng1=100,nr1=100,np1=5) double precision gb(ni1,2) double precision g2(ni1,0:nb1-1) double precision dg2(ni1,0:nb1-1) double precision nphot2(ni1,nb1) double precision nphot1(ni1,ngp1,nc1,ng1) double precision g1(ni1,ngp1,nc1,0:ng1-1) double precision dg1(ni1,ngp1,nc1,0:ng1-1) double precision tf(ni1,ngp1,ng1+1,nr1+1,nc1,np1) logical gmask2(ni1,0:nb1-1) logical gpmask(ni1,ngp1),amask(ni1) logical gmask1(ni1,ngp1,nc1,0:ng1-1) logical tfmask(ni1,ngp1,ng1+1,nr1+1,nc1) logical check1,check2 character*80 buff integer ftype r1=rmin r2=rmax call adet (a,rmin) !write(*,*)'Spin of black hole',a ftype=0 call tfset (check1,amask,gpmask,tfmask,tf,rmin,rmax,a,obs, & ftype,ni1,ngp1,ng1,nr1,nc1,np1) !write(*,*)'Completed transfer function read, result',check1 checkout=.FALSE. If (check1) then call ldet (check2,gmask1,nphot1,g1,dg1,amask,gpmask,tfmask,tf, & index,a,rtype,mtype,ni,ngp,ng,nr,nc,np) !write(*,*)'Determined line profiles, result',check2 If (check2) then call lmergep1 (gmask2,gb,nphot2,g2,dg2,amask,gpmask, & gmask1,nphot1,g1,dg1,nb1,ni1,ngp1,nc1,ng1) !write(*,*)'Merged line profiles onto common energy grid' call lmergep2 (checkout,nmask,n,e,de,amask, & gmask2,gb,nphot2,g2,dg2,a,nb,ni1,nb1) !write(*,*)'Determined line profile for spin' End If End If return end subroutine lineinit1 (checkout,nmask,n,e,de,r1,r2,rmin,rmax, & a,as,obs,index,rtype,mtype,nb,ni,ngp,ng,nr,nc,np) implicit none integer ngp,ng,nr,nc,np integer rtype,mtype,nb,ni double precision r1,r2 double precision rmin,rmax double precision de(0:nb-1) double precision a,as,obs,index double precision e(0:nb-1),n(nb) logical checkout,nmask(0:nb-1) integer nb1,ni1,ngp1,ng1,nr1,nc1,np1 parameter (nb1=100,ni1=2,ngp1=5,nc1=2,ng1=100,nr1=100,np1=5) double precision gb(ni1,2) double precision risco,rms double precision g2(ni1,0:nb1-1) double precision dg2(ni1,0:nb1-1) double precision nphot2(ni1,nb1) double precision nphot1(ni1,ngp1,nc1,ng1) double precision g1(ni1,ngp1,nc1,0:ng1-1) double precision dg1(ni1,ngp1,nc1,0:ng1-1) double precision tf(ni1,ngp1,ng1+1,nr1+1,nc1,np1) logical gmask2(ni1,0:nb1-1) logical gpmask(ni1,ngp1),amask(ni1) logical gmask1(ni1,ngp1,nc1,0:ng1-1) logical tfmask(ni1,ngp1,ng1+1,nr1+1,nc1) logical check1,check2,check4 character*80 buff integer ftype check1=.FALSE. call rconvert (check1,r1,r2,rmin,rmax,a) !write(*,*)'Converted radii to gravitational units',r1,r2,check1 If (check1) then risco = rms(a) If (r1.LT.risco) r1 = risco !write(*,*)'New r_{ms}',risco !write(*,*)'Revised inner edge of disc',r1 !write(*,*)'Spin of black hole',a !write(*,*)'Inner edge of disc, Rin(G)',r1 !write(*,*)'Outer edge of disc, Rout(G)',r2 ftype=1 call tfset (check2,amask,gpmask,tfmask,tf,r1,rmax,a,obs, & ftype,ni1,ngp1,ng1,nr1,nc1,np1) !write(*,*)'Completed transfer function read, result',check2 checkout=.FALSE. If (check2) then call ldet (check4,gmask1,nphot1,g1,dg1,amask,gpmask,tfmask, & tf,index,a,rtype,mtype,ni,ngp,ng,nr,nc,np) !write(*,*)'Determined line profiles, result',check4 If (check4) then call lmergep1 (gmask2,gb,nphot2,g2,dg2,amask,gpmask, & gmask1,nphot1,g1,dg1,nb1,ni1,ngp1,nc1,ng1) !write(*,*)'Merged line profiles onto common energy grid' call lmergep2 (checkout,nmask,n,e,de,amask, & gmask2,gb,nphot2,g2,dg2,a,nb,ni1,nb1) !write(*,*)'Determined line profile for spin' End If End If End If return end subroutine lineinit2 (checkout,nmask,n,e,de,r1,r2,rmin,rmax, & a,as,obs,index,rtype,mtype,nb,ni,ngp,ng,nr,nc,np) implicit none integer ngp,ng,nr,nc,np integer rtype,mtype,nb,ni double precision r1,r2 double precision rmin,rmax double precision de(0:nb-1) double precision a,as,obs,index double precision e(0:nb-1),n(nb) logical checkout,nmask(0:nb-1) integer nb1,ni1,ngp1,ng1,nr1,nc1,np1 parameter (nb1=100,ni1=2,ngp1=5,nc1=2,ng1=100,nr1=100,np1=5) double precision gb(ni1,2) double precision risco,rms double precision g2(ni1,0:nb1-1) double precision dg2(ni1,0:nb1-1) double precision nphot2(ni1,nb1) double precision aval(ni1),rval(ni1) double precision nphot1(ni1,ngp1,nc1,ng1) double precision g1(ni1,ngp1,nc1,0:ng1-1) double precision dg1(ni1,ngp1,nc1,0:ng1-1) double precision tf(ni1,ngp1,ng1+1,nr1+1,nc1,np1) logical gmask2(ni1,0:nb1-1) logical gpmask(ni1,ngp1),amask(ni1) logical gmask1(ni1,ngp1,nc1,0:ng1-1) logical tfmask(ni1,ngp1,ng1+1,nr1+1,nc1) logical check1,check2,check4 integer ftype call aclose (amask,rval,aval,a,ni1,ni1) write(*,*)'Closest values in a',aval(1),aval(2) check1=amask(1).AND.amask(2) If (check1) then r1=dabs(aval(1)-a) r2=dabs(aval(2)-a) If (r1.LE.r2) amask(2)=.FALSE. If (r1.GT.r2) amask(1)=.FALSE. If (amask(1)) a=aval(1) If (amask(2)) a=aval(2) If (amask(1)) as=aval(1) If (amask(2)) as=aval(2) write(*,*)'Spin parameters',a,aval(1),aval(2),amask(1),amask(2) End If check1=.FALSE. call rconvert (check1,r1,r2,rmin,rmax,a) write(*,*)'Converted radii to gravitational units',r1,r2,check1 If (check1) then risco = rms(a) If (r1.LT.risco) r1 = risco write(*,*)'New r_{ms}',risco write(*,*)'Revised inner edge of disc',r1 ftype=2 call tfset (check2,amask,gpmask,tfmask,tf,r1,rmax,a,obs, & ftype,ni1,ngp1,ng1,nr1,nc1,np1) write(*,*)'Completed transfer function read, result',check2 checkout=.FALSE. If (check2) then call ldet (check4,gmask1,nphot1,g1,dg1,amask,gpmask,tfmask, & tf,index,a,rtype,mtype,ni,ngp,ng,nr,nc,np) write(*,*)'Determined line profiles, result',check4 If (check4) then call lmergep1 (gmask2,gb,nphot2,g2,dg2,amask,gpmask, & gmask1,nphot1,g1,dg1,nb1,ni1,ngp1,nc1,ng1) write(*,*)'Merged line profiles onto common energy grid' call lmergep2 (checkout,nmask,n,e,de,amask, & gmask2,gb,nphot2,g2,dg2,a,nb,ni1,nb1) write(*,*)'Determined line profile for spin' End If End If End If return end c ************************************************************************* c * These are the subroutines for lineinit * c ************************************************************************* subroutine adet (a,rmin) implicit none double precision a,rmin double precision a1,a2 double precision aval(101) double precision rval(101) double precision scale,rms logical check1,mask1(101) integer i,n n=101 a1=1.0d0 a2=0.0d0 scale=1.0d-2 call gridinit (aval,a1,a2,scale,n) !write(*,*)'Setup spin grid' do i=1,n mask1(i)=.TRUE. rval(i)=rms(aval(i)) end do !write(*,*)'Initialised r_{ms} grid' check1=.FALSE. call linterpolate (check1,a,rmin,mask1,rval,aval,n) !write(*,*)'Interpolated value of a',rmin,a,check1 return end subroutine rconvert (checkout,r1,r2,rmin,rmax,a) implicit none double precision r1,r2,rmin,rmax,a logical checkout integer n parameter (n=1000) double precision xmin,xmax double precision risco,rms,rout double precision x(n),y(n),scale logical mask1(n),check1,check2 integer i xmin=1.0d0 xmax=1.0d2 rout=1.0d2 scale=1.0d-3 risco = rms(a) call gridinit (x,xmin,xmax,scale,n) call gridinit (y,risco,rout,scale,n) !write(*,*)'Initialised radial conversion grids' do i=1,n mask1(i)=.TRUE. !write(60,*)i,x(i),y(i),mask1(i) end do !write(*,*)'Initialised interpolation mask' check1=.FALSE. check2=.FALSE. call linterpolate (check1,r1,rmin,mask1,x,y,n) call linterpolate (check2,r2,rmax,mask1,x,y,n) !write(*,*)'Gravitational radii',r1,r2,check1,check2 checkout=check1.AND.check2 !write(*,*)'Were conversions succesfull?',checkout return end subroutine tfset (checkout,amask,gpmask,tfmask,tf,rmin,rmax,a,obs, & ftype,ni,ngp,ng,nr,nc,np) implicit none integer ftype integer ni,ngp,ng,nr,nc,np double precision rmin,rmax,a,obs double precision tf(ni,ngp,ng+1,nr+1,nc,np) logical tfmask(ni,ngp,ng+1,nr+1,nc) logical gpmask(ni,ngp),amask(ni) logical checkout integer ni1,ngp1,ng1,nr1,nc1,np1 parameter (ni1=2,ngp1=5,nc1=2,ng1=100,nr1=100,np1=5) double precision aval(ni1),rval(ni1),r1,r2,a1 double precision work5(ngp1,ng1+1,nr1+1,nc1,np1) logical check1,check2,mask1(ngp1) logical mask4(ngp1,ng1+1,nr1+1,nc1) integer n1,n2 integer i,j,k,l,m,n call aclose (amask,rval,aval,a,ni,ni1) !write(*,*)'Closest values in a',aval(1),aval(2) check1=amask(1).AND.amask(2) If ((ftype.EQ.2).AND.check1) then r1=dabs(aval(1)-a) r2=dabs(aval(2)-a) If (r1.LE.r2) amask(2)=.FALSE. If (r1.GT.r2) amask(1)=.FALSE. !write(*,*)'Spin parameters',aval(1),aval(2),amask(1),amask(2) End If n1=0 n2=0 do i=1,ni check1=amask(i) If (check1) then r1=rmin r2=rmax n1=n1+1 a1=aval(i) check2=.FALSE. If (ftype.EQ.0) r1=rval(i) call tfread (check2,mask1,mask4,work5,r1,r2,a1,obs, & ngp1,nc1,ng1,nr1,np1) !write(*,*)i,' Read in transfer function, result',check2 amask(i)=check2 If (amask(i)) then call rminset (mask1,mask4,work5,r1,a1, & ngp1,nc1,ng1,nr1,np1) !write(*,*)i,' Regridded inner edge' call rmaxset (mask1,mask4,work5,r2,a1, & ngp1,nc1,ng1,nr1,np1) !write(*,*)i,' Regridded outer edge' n2=n2+1 do j=1,ngp gpmask(i,j)=mask1(j) If (mask1(j)) then do k=1,ng+1 do l=1,nr+1 do m=1,nc tfmask(i,j,k,l,m)=mask4(j,k,l,m) If (mask4(j,k,l,m)) then do n=1,np tf(i,j,k,l,m,n)=work5(j,k,l,m,n) end do End If end do end do end do End If end do !write(*,*)i,' Placed transfer function in storage arrays' End If End If end do checkout=.FALSE. If (n1.EQ.n2) checkout=.TRUE. !write(*,*)'Were transfer functions read in successfully?',checkout return end subroutine ldet (checkout,gmask,nphot,g,dg,amask,gpmask,tfmask,tf, & index,a,rtype,mtype,ni,ngp,ng,nr,nc,np) implicit none integer rtype,mtype integer ni,ngp,ng,nr,nc,np double precision rmin,rmax double precision a,obs,index double precision nphot(ni,ngp,nc,ng) double precision g(ni,ngp,nc,0:ng-1) double precision dg(ni,ngp,nc,0:ng-1) double precision tf(ni,ngp,ng+1,nr+1,nc,np) logical tfmask(ni,ngp,ng+1,nr+1,nc) logical gpmask(ni,ngp),amask(ni) logical gmask(ni,ngp,nc,0:ng-1) logical checkout integer ni1,ngp1,ng1,nr1,nc1,np1 parameter (ni1=2,ngp1=5,nc1=2,ng1=100,nr1=100,np1=5) double precision work3a(ngp1,nc1,ng1) double precision work3b(ngp1,nc1,0:ng1-1) double precision work3c(ngp1,nc1,0:ng1-1) double precision work5(ngp1,ng1+1,nr1+1,nc1,np1) logical mask1(ngp1),check1 logical mask3(ngp1,nc1,0:ng1-1) logical mask4(ngp1,ng1+1,nr1+1,nc1) integer n1,n2 integer i,j,k,l,m,n n1=0 n2=0 do i=1,ni check1=amask(i) If (check1) then n1=n1+1 do j=1,ngp mask1(j)=gpmask(i,j) do k=1,ng+1 do l=1,nr+1 do m=1,nc mask4(j,k,l,m)=tfmask(i,j,k,l,m) do n=1,np work5(j,k,l,m,n)=tf(i,j,k,l,m,n) end do end do end do end do end do !write(*,*)i,' Extracted transfer function' amask(i)=.FALSE. call lcalc (amask(i),work3b,work3a,work3c, & mask3,mask1,mask4,work5, & index,a,rtype,mtype,ngp1,nc1,ng1,nr1,np1) !write(*,*)i,' Calculated line profiles, result',amask(i) If (amask(i)) then n2=n2+1 do j=1,ngp If (gpmask(i,j)) then do k=1,nc do l=1,ng nphot(i,j,k,l)=work3a(j,k,l) g(i,j,k,l-1)=work3b(j,k,l-1) dg(i,j,k,l-1)=work3c(j,k,l-1) gmask(i,j,k,l-1)=mask3(j,k,l-1) !write(44,*)i,j,k,l,dg(i,j,k,l-1) end do end do End If end do !write(*,*)i,' Stored line profile' End If End If end do checkout=.FALSE. If (n1.EQ.n2) checkout=.TRUE. !write(*,*)'Were line profiles calculated successfully?',checkout return end subroutine lmergep1 (gmask2,gb,nphot2,g2,dg2,amask,gpmask, & gmask1,nphot1,g1,dg1,nb,ni,ngp,nc,ng) implicit none integer nb,ni,ngp,nc,ng double precision gb(ni,2) double precision g2(ni,0:nb-1) double precision dg2(ni,0:nb-1) double precision nphot2(ni,nb) double precision nphot1(ni,ngp,nc,ng) double precision g1(ni,ngp,nc,0:ng-1) double precision dg1(ni,ngp,nc,0:ng-1) logical gmask1(ni,ngp,nc,0:ng-1) logical gpmask(ni,ngp),amask(ni) logical gmask2(ni,0:nb-1) integer nb1,ni1,ngp1,nc1,ng1 parameter (nb1=100,ni1=2,ngp1=5,nc1=2,ng1=100) double precision gmin,gmax double precision work1a(nb1) double precision work1b(0:nb-1) double precision work1c(0:nb-1) double precision work3a(ngp1,nc1,ng1) double precision work3b(ngp1,nc1,0:ng1-1) double precision work3c(ngp1,nc1,0:ng1-1) logical mask1a(ngp) logical mask1b(0:nb-1) logical mask3(ngp1,nc1,0:ng1-1) integer i,j,k,l gmin=1.0d2 gmax=-1.0d2 do i=1,ni gb(i,1)=0.0d0 gb(i,2)=0.0d0 If (amask(i)) then gb(i,1)=1.0d2 gb(i,2)=-1.0d2 do j=1,ngp If (gpmask(i,j)) then do k=1,nc gmin=min(gmin,g1(i,j,k,0)) gmax=max(gmax,g1(i,j,k,ng-1)) gb(i,1)=min(gb(i,1),g1(i,j,k,0)) gb(i,2)=max(gb(i,2),g1(i,j,k,ng-1)) end do End If end do End If end do !write(*,*)'Lower, upper limits of line',gmin,gmax do i=1,ni If (amask(i)) then do j=1,ngp mask1a(j)=gpmask(i,j) If (gpmask(i,j)) then do k=1,nc do l=1,ng work3a(j,k,l)=nphot1(i,j,k,l) work3b(j,k,l-1)=g1(i,j,k,l-1) work3c(j,k,l-1)=dg1(i,j,k,l-1) mask3(j,k,l-1)=gmask1(i,j,k,l-1) !write(46,*)i,j,k,l,work3b(j,k,l-1), ! & work3a(j,k,l),work3c(j,k,l-1), ! & mask3(j,k,l-1),mask1a(j) end do end do End If end do !write(*,*)i,' Extracted line profile' call rebinp1 (mask1b,work1a,work1b,work1c,gmin,gmax, & mask1a,mask3,work3a,work3c,work3b, & nb1,ngp1,ng1,nc1) !write(*,*)i,' Rebinned line profiles' call lsmooth (mask1b,work1a,work1b,work1c,mask1a, & mask3,work3b,nb1,ngp1,ng1,nc1) !write(*,*)i,' Smoothed line boundaries' do j=1,nb nphot2(i,j)=work1a(j) g2(i,j-1)=work1b(j-1) dg2(i,j-1)=work1c(j-1) gmask2(i,j-1)=mask1b(j-1) !write(40,*)i,j,work1b(j-1),work1c(j-1), ! & work1a(j),mask1b(j-1) !write(41,*)i,j,g2(i,j-1),dg2(i,j-1), ! & nphot2(i,j),gmask2(i,j-1) ! If (gmask2(i,j-1)) !write(42,*)g2(i,j-1),nphot2(i,j) end do !write(*,*)i,' Stored line profiles' !write(40,*) !write(41,*) !write(42,*)'no no' End If end do return end subroutine lmergep2 (checkout,nmask,n,e,de,amask, & gmask,gb,nphot,g,dg,a,nb,ni,ng) implicit none integer nb,ni,ng double precision gb(ni,2) double precision nphot(ni,ng) double precision de(0:nb-1),a double precision g(ni,0:ng-1) double precision dg(ni,0:ng-1) double precision e(0:nb-1),n(nb) logical amask(ni),gmask(ni,0:nb-1) logical checkout,nmask(0:nb-1) integer nb1,ni1,ng1 parameter (nb1=100,ni1=2,ng1=100) double precision work1a(nb1) double precision work1d(nb1) double precision work1e(nb1) double precision work1f(nb1) double precision gmin,gmax,de1 double precision work1b(0:nb1-1) double precision work1c(0:nb1-1) double precision rfrac,risco,rms,x,y double precision aval(ni1),rval(ni1) logical mask1a(ni1) logical mask1b(0:nb1-1) logical mask1c(nb1),check1,check2 integer i,nvalid risco=rms(a) !write(*,*)'r_{ms}(a)',risco call aclose (mask1a,rval,aval,a,ni1,ni1) !write(*,*)'Closest values of r_{ms}(a)',rval(1),rval(2) rfrac=1.0d0 If (amask(1).AND.amask(2)) rfrac=(risco-rval(1))/(rval(2)-rval(1)) !write(*,*)'Fractional distance between points in r_{ms}',rfrac gmin=g(1,0) gmax=g(1,ng-1) de1=(gmax-gmin)/dble(nb1-1) !write(*,*)'Lower, upper limits of line',gmin,gmax !write(*,*)'Bin width',de1 nvalid=0 work1b(0)=gmin do i=1,nb1 work1a(i)=0.0d0 work1c(i-1)=de1 mask1b(i-1)=.FALSE. If (i.GT.1) work1b(i-1)=work1b(i-2)+de1 If (amask(1).AND.amask(2)) then If (gmask(1,i-1).AND.gmask(2,i-1)) then nvalid=nvalid+1 mask1b(i-1)=.TRUE. work1a(i)=min(nphot(1,i),nphot(2,i))+ & dabs(rfrac*(nphot(2,i)-nphot(1,i))) Else If (gmask(1,i-1).AND.((work1b(i-1).LT.gb(2,1)).OR. & (work1b(i-1).GT.gb(2,2)))) then nvalid=nvalid+1 mask1b(i-1)=.TRUE. work1a(i)=rfrac*nphot(1,i) Else if (gmask(2,i-1).AND.((work1b(i-1).LT.gb(1,1)).OR. & (work1b(i-1).GT.gb(1,2)))) then nvalid=nvalid+1 mask1b(i-1)=.TRUE. work1a(i)=rfrac*nphot(2,i) End If End If !write(48,*)i,work1b(i-1),nphot(1,i),work1a(i),nphot(2,i), ! & gmask(1,i-1),mask1b(i-1),gmask(2,i-1) Else if (amask(1)) then If (gmask(1,i-1)) then nvalid=nvalid+1 work1a(i)=nphot(1,i) mask1b(i-1)=.TRUE. End If Else if (amask(2)) then If (gmask(2,i-1)) then nvalid=nvalid+1 work1a(i)=nphot(2,i) mask1b(i-1)=.TRUE. End If End If end do check1=.FALSE. If (nvalid.GT.0) check1=.TRUE. !write(*,*)'Were valid points found in line?',check1 nvalid=0 If (check1) then gmin=work1b(0) gmax=work1b(nb1-1) If (amask(1).AND.amask(2)) then gmin=min(gb(1,1),gb(2,1))+dabs(rfrac*(gb(2,1)-gb(1,1))) gmax=max(gb(1,1),gb(2,2))+dabs(rfrac*(gb(2,2)-gb(1,2))) do i=1,nb1 mask1c(i)=mask1b(i-1) work1f(i)=work1c(i-1) work1d(i)=work1a(i)/work1c(i-1) work1e(i)=work1b(i-1)+0.5d0*work1c(i-1) end do !write(*,*)'Setup interpolation arrays' End If !write(*,*)'Lower, upper limits on line',gmin,gmax de1=(gmax-gmin)/dble(nb1-1) !write(*,*)'Bin width',de1 e(0)=gmin do i=1,nb n(i)=0.0d0 de(i-1)=de1 nmask(i-1)=.FALSE. If (i.GT.1) e(i-1)=e(i-2)+de1 If (amask(1).AND.amask(2)) then y=0.0d0 check2=.FALSE. x=e(i-1)+0.5d0*de1 call linterpolate (check2,y,x,mask1c,work1e,work1d,nb1) !write(*,*)i,' Line interpolation',x,y,check2 If (check2) then n(i)=y*de1 nvalid=nvalid+1 nmask(i-1)=.TRUE. End If Else If (mask1b(i-1)) then n(i)=work1a(i) nvalid=nvalid+1 nmask(i-1)=.TRUE. End If End If !write(*,*)i,e(i-1),de(i-1),n(i),nmask(i-1) end do End If checkout=.FALSE. If (nvalid.GT.0) checkout=.TRUE. !write(*,*)'Were valid points found in line?',checkout return end subroutine aclose (amask,rval,aval,a,ni,ni1) implicit none integer ni,ni1 double precision rval(ni1),aval(ni1),a logical amask(ni) double precision aquant(12),rms,tol integer i do i=1,ni aval(i)=0.0d0 rval(i)=0.0d0 amask(i)=.FALSE. end do !write(*,*)'Initialised closest spin arrays' do i=1,12 If (i.EQ.1) aquant(i)=0.0d0 If (i.EQ.2) aquant(i)=0.1501819899670580d0 If (i.EQ.3) aquant(i)=0.2940047205828976d0 If (i.EQ.4) aquant(i)=0.4305113630911228d0 If (i.EQ.5) aquant(i)=0.5584815599030040d0 If (i.EQ.6) aquant(i)=0.6763195721910966d0 If (i.EQ.7) aquant(i)=0.7818758452482436d0 If (i.EQ.8) aquant(i)=0.8721520262014855d0 If (i.EQ.9) aquant(i)=0.9428090365807393d0 If (i.EQ.10) aquant(i)=0.9874884803808748d0 If (i.EQ.11) aquant(i)=0.998d0 If (i.EQ.12) aquant(i)=1.0d0 end do !write(*,*)'Initialised transfer function spins' tol=1.0d-5 do i=1,12 If (dabs(a-aquant(i)).LT.tol) then amask(1)=.TRUE. aval(1)=aquant(i) rval(1)=rms(aval(1)) !write(*,*)'Closest spin',aval(1) exit Else if (i.LT.12) then If ((aquant(i).LT.a).AND.(aquant(i+1).GT.a)) then If (dabs(a-aquant(i+1)).LT.tol) then amask(1)=.TRUE. aval(1)=aquant(i+1) rval(1)=rms(aval(1)) !write(*,*)'Closest spin',aval(1) Else amask(1)=.TRUE. amask(2)=.TRUE. aval(1)=aquant(i) aval(2)=aquant(i+1) rval(1)=rms(aval(1)) rval(2)=rms(aval(2)) !write(*,*)'Closest spins',aval(1),aval(2) End If exit End If End If end do return end c ************************************************************************* c * These are the generic subroutines * c ************************************************************************* subroutine tfread (checkout,gpmask,tfmask,tf, & rmin,rmax,a,obs,ngp,nc,ng,nr,np) implicit none integer ngp,nc,ng,nr,np logical checkout,gpmask(ngp) logical tfmask(ngp,ng+1,nr+1,nc) double precision rmin,rmax,a,obs double precision tf(ngp,ng+1,nr+1,nc,np) integer ng1,nr1,nc1 parameter (nc1=2,ng1=100,nr1=100) double precision rms,r1,r2,tol double precision v(ng1+1,nr1+1,nc1) double precision w(ng1+1,nr1+1,nc1) double precision x(ng1+1,nr1+1,nc1) double precision y(ng1+1,nr1+1,nc1) double precision z(ng1+1,nr1+1,nc1) character*256 path logical check1 logical check2 logical check4 logical m3(ng1+1,nr1+1,nc1) integer nvalid,loc integer start,status integer i,j,k,l start=1 nvalid = 0 tol = 1.0d-6 gpmask(1)=.FALSE. If (dabs(a).LT.tol) then start=2 gpmask(1)=.FALSE. End If !write(*,*)'Line paramaters',rmin,rmax,a,obs do i=start,ngp r1 = rms(a) If (i.EQ.1) r2 = 6.0d0 If (i.EQ.2) r1 = 6.0d0 If (i.EQ.2) r2 = 1.0d1 If (i.EQ.3) r1 = 1.0d1 If (i.EQ.3) r2 = 2.0d1 If (i.EQ.4) r1 = 2.0d1 If (i.EQ.4) r2 = 5.0d1 If (i.EQ.5) r1 = 5.0d1 If (i.EQ.5) r2 = 1.0d2 !write(*,*)i,' Lower limits of radial grouping',r1,rmin,r2 !write(*,*)i,' Upper limits of radial grouping',r1,rmax,r2 check1=.FALSE. check2=.FALSE. check4=.FALSE. If ((rmin.LT.r1).AND.(rmin.LT.r2).AND. & (rmax.GT.r1).AND.(rmax.GT.r2)) check1=.TRUE. If ((rmin.GE.r1).AND.(rmin.LT.r2)) check2=.TRUE. If ((rmax.GT.r1).AND.(rmax.LE.r2)) check4=.TRUE. If (check2.OR.check4) check1=.TRUE. !write(*,*)i,' Does group contribute?',check1,check2,check4 gpmask(i)=.FALSE. If (check1) then call pathinit (path,r1,r2,a,obs) c write(*,*)i,' Path to file ',path check2=.FALSE. call tfopen (check2,m3,v,w,x,y,z,path,ng1,nr1,nc1) c write(*,*)i,' Read in transfer function' If (check2) then gpmask(i)=.TRUE. do j=1,ng+1 do k=1,nr+1 do l=1,nc tf(i,j,k,l,1)=v(j,k,l) tf(i,j,k,l,2)=w(j,k,l) tf(i,j,k,l,3)=x(j,k,l) tf(i,j,k,l,4)=y(j,k,l) tf(i,j,k,l,5)=z(j,k,l) tfmask(i,j,k,l)=m3(j,k,l) end do end do end do End If !write(*,*)i,' Read in transfer function, result',gpmask(i) End If If (gpmask(i)) nvalid = nvalid + 1 end do checkout=.FALSE. If (nvalid.GT.0) checkout=.TRUE. !write(*,*)'Was transfer function read succesfull?',checkout return end subroutine rminset (gpmask,tfmask,tf,rmin,a,ngp,nc,ng,nr,np) implicit none integer ngp,nc,ng,nr,np double precision rmin,a logical checkout,gpmask(ngp) logical tfmask(ngp,ng+1,nr+1,nc) double precision tf(ngp,ng+1,nr+1,nc,np) integer nr1 parameter (nr1=101) double precision r1,r2,scale,rms double precision r(nr1),m(nr1) double precision x(nr1),y(nr1) double precision z(nr1),w(3) logical mask1b(3) logical mask1a(nr1) logical check1,check2 logical check4,check5 logical check6,check7 integer i,j,k,l do i=1,ngp If (gpmask(i)) then r1 = rms(a) If (i.EQ.1) r2 = 6.0d0 If (i.EQ.2) r1 = 6.0d0 If (i.EQ.2) r2 = 1.0d1 If (i.EQ.3) r1 = 1.0d1 If (i.EQ.3) r2 = 2.0d1 If (i.EQ.4) r1 = 2.0d1 If (i.EQ.4) r2 = 5.0d1 If (i.EQ.5) r1 = 5.0d1 If (i.EQ.5) r2 = 1.0d2 !write(*,*)i,' Lower,upper limits of radial grouping',r1,r2 check1=.FALSE. If (rmin.GT.r1) check1=.TRUE. !write(*,*)i,' Does grid need reshaping?',r1,rmin,r2,check1 If (check1) then do j=1,nc do k=1,ng r1=0.0d0 check2=.FALSE. do l=1,nr+1 If (tfmask(i,k,l,j)) then check2=.TRUE. r1=tf(i,k,l,j,2) exit End If end do !write(*,*)i,j,k,' Inner boundary',r1,check2 r2=0.0d0 check4=.FALSE. do l=nr+1,1,-1 If (tfmask(i,k,l,j)) then check4=.TRUE. r2=tf(i,k,l,j,2) exit End If end do !write(*,*)i,j,k,' Outer boundary',r2,check4 !write(50,*)i,j,k,r1,r2,check2,check4 check5=.FALSE. check6=.FALSE. check7=.FALSE. If (check2.AND.check4) then If ((r1.GT.rmin).AND.(r2.GT.rmin)) check5=.TRUE. If ((r1.LT.rmin).AND.(r2.GT.rmin)) check6=.TRUE. If ((r1.LT.rmin).AND.(r2.LT.rmin)) check7=.TRUE. End If !write(*,*)i,j,k,' Results of boundary tests', ! & check5,check6,check7 !write(51,*)i,j,k,r1,rmin,r2,check5,check6,check7 If (check6) then do l=1,nr1 r(l)=tf(i,k,l,j,2) m(l)=tf(i,k,l,j,3) x(l)=tf(i,k,l,j,4) y(l)=tf(i,k,l,j,5) mask1a(l)=tfmask(i,k,l,j) end do !write(*,*)i,j,k,' Setup interpolation arrays' scale=1.0d-2 call gridinit (z,rmin,r2,scale,nr1) !write(*,*)i,j,k,'Initialised new radial grid' do l=1,nr1 w(1)=0.0d0 w(2)=0.0d0 w(3)=0.0d0 mask1b(1)=.FALSE. mask1b(2)=.FALSE. mask1b(3)=.FALSE. tf(i,k,l,j,2)=0.0d0 tf(i,k,l,j,3)=0.0d0 tf(i,k,l,j,4)=0.0d0 tf(i,k,l,j,5)=0.0d0 tfmask(i,k,l,j)=.FALSE. call linterpolate (mask1b(1),w(1),z(l), & mask1a,r,m,nr1) call linterpolate (mask1b(2),w(2),z(l), & mask1a,r,x,nr1) call linterpolate (mask1b(3),w(3),z(l), & mask1a,r,y,nr1) If (mask1b(1).AND.mask1b(2) & .AND.mask1b(3)) then tfmask(i,k,l,j)=.TRUE. tf(i,k,l,j,2)=z(l) tf(i,k,l,j,3)=w(1) tf(i,k,l,j,4)=w(2) tf(i,k,l,j,5)=w(3) End If end do Else if (check7) then do l=1,nr+1 tf(i,k,l,j,2)=0.0d0 tf(i,k,l,j,3)=0.0d0 tf(i,k,l,j,4)=0.0d0 tf(i,k,l,j,5)=0.0d0 tfmask(i,k,l,j)=.FALSE. end do End If end do end do End If exit End If end do return end subroutine rmaxset (gpmask,tfmask,tf,rmax,a,ngp,nc,ng,nr,np) implicit none integer ngp,nc,ng,nr,np double precision rmax,a logical checkout,gpmask(ngp) logical tfmask(ngp,ng+1,nr+1,nc) double precision tf(ngp,ng+1,nr+1,nc,np) integer nr1 parameter (nr1=101) double precision r1,r2,scale,rms double precision r(nr1),m(nr1) double precision x(nr1),y(nr1) double precision z(nr1),w(3) logical mask1b(3) logical mask1a(nr1) logical check1,check2 logical check4,check5 logical check6,check7 integer i,j,k,l do i=ngp,1,-1 If (gpmask(i)) then r1 = rms(a) If (i.EQ.1) r2 = 6.0d0 If (i.EQ.2) r1 = 6.0d0 If (i.EQ.2) r2 = 1.0d1 If (i.EQ.3) r1 = 1.0d1 If (i.EQ.3) r2 = 2.0d1 If (i.EQ.4) r1 = 2.0d1 If (i.EQ.4) r2 = 5.0d1 If (i.EQ.5) r1 = 5.0d1 If (i.EQ.5) r2 = 1.0d2 !write(*,*)i,' Lower,upper limits of radial grouping',r1,r2 check1=.FALSE. If (rmax.LT.r2) check1=.TRUE. !write(*,*)i,' Does grid need reshaping?',r1,rmax,r2,check1 If (check1) then do j=1,nc do k=1,ng r1=0.0d0 check2=.FALSE. do l=1,nr+1 If (tfmask(i,k,l,j)) then check2=.TRUE. r1=tf(i,k,l,j,2) exit End If end do !write(*,*)i,j,k,' Inner boundary',r1,check2 r2=0.0d0 check4=.FALSE. do l=nr+1,1,-1 If (tfmask(i,k,l,j)) then check4=.TRUE. r2=tf(i,k,l,j,2) exit End If end do !write(*,*)i,j,k,' Outer boundary',r2,check4 !write(50,*)i,j,k,r1,r2,check2,check4 check5=.FALSE. check6=.FALSE. check7=.FALSE. If (check2.AND.check4) then If ((r1.LT.rmax).AND.(r2.LT.rmax)) check5=.TRUE. If ((r1.LT.rmax).AND.(r2.GT.rmax)) check6=.TRUE. If ((r1.GT.rmax).AND.(r2.GT.rmax)) check7=.TRUE. End If !write(*,*)i,j,k,' Results of oundary tests', ! & check5,check6,check7 !write(51,*)i,j,k,r1,rmax,r2,check5,check6,check7 If (check6) then do l=1,nr1 r(l)=tf(i,k,l,j,2) m(l)=tf(i,k,l,j,3) x(l)=tf(i,k,l,j,4) y(l)=tf(i,k,l,j,5) mask1a(l)=tfmask(i,k,l,j) end do !write(*,*)i,j,k,' Setup interpolation arrays' scale=1.0d-2 call gridinit (z,r1,rmax,scale,nr1) !write(*,*)i,j,k,'Initialised new radial grid' do l=1,nr1 w(1)=0.0d0 w(2)=0.0d0 w(3)=0.0d0 mask1b(1)=.FALSE. mask1b(2)=.FALSE. mask1b(3)=.FALSE. tf(i,k,l,j,2)=0.0d0 tf(i,k,l,j,3)=0.0d0 tf(i,k,l,j,4)=0.0d0 tf(i,k,l,j,5)=0.0d0 tfmask(i,k,l,j)=.FALSE. call linterpolate (mask1b(1),w(1),z(l), & mask1a,r,m,nr1) call linterpolate (mask1b(2),w(2),z(l), & mask1a,r,x,nr1) call linterpolate (mask1b(3),w(3),z(l), & mask1a,r,y,nr1) If (mask1b(1).AND.mask1b(2) & .AND.mask1b(3)) then tfmask(i,k,l,j)=.TRUE. tf(i,k,l,j,2)=z(l) tf(i,k,l,j,3)=w(1) tf(i,k,l,j,4)=w(2) tf(i,k,l,j,5)=w(3) End If end do Else if (check7) then do l=1,nr+1 tf(i,k,l,j,2)=0.0d0 tf(i,k,l,j,3)=0.0d0 tf(i,k,l,j,4)=0.0d0 tf(i,k,l,j,5)=0.0d0 tfmask(i,k,l,j)=.FALSE. end do End If end do end do End If exit End If end do return end subroutine lcalc (checkout,g,nphot,dg,gmask,gpmask,tfmask,tf, & index,a,rtype,mtype,ngp,nc,ng,nr,np) implicit none logical checkout integer rtype,mtype double precision index,a integer ngp,nc,ng,nr,np double precision nphot(ngp,nc,ng) double precision g(ngp,nc,0:ng-1) double precision dg(ngp,nc,0:ng-1) double precision tf(ngp,ng+1,nr+1,nc,np) logical tfmask(ngp,ng+1,nr+1,nc) logical gmask(ngp,nc,0:ng-1) logical gpmask(ngp) double precision gval,xy(2,2,4),area,emm logical check1 integer nvalid,counter integer i,j,k,l,m,n counter=0 do i=1,ngp If (gpmask(i)) then do j=1,nc do k=1,ng nphot(i,j,k) = 0.0d0 gmask(i,j,k-1)=.FALSE. g(i,j,k-1) = tf(i,k,1,j,1) If (k.EQ.ng) g(i,j,k-1)=tf(i,k+1,1,j,1) dg(i,j,k-1) = tf(i,k+1,1,j,1)-tf(i,k,1,j,1) gval = tf(i,k,1,j,1) + 0.5d0 * dg(i,j,k-1) !write(*,*)i,j,k,' Redshift',gval do l=1,nr !write(*,*)i,j,k,l,' Starting iteration' nvalid=0 do m=1,2 do n=1,2 xy(m,n,1)=tf(i,m+k-1,n+l-1,j,4) xy(m,n,2)=tf(i,m+k-1,n+l-1,j,5) xy(m,n,3)=tf(i,m+k-1,n+l-1,j,2) xy(m,n,4)=tf(i,m+k-1,n+l-1,j,3) If (tfmask(i,m+k-1,n+l-1,j)) nvalid=nvalid+1 end do end do !write(*,*)i,j,k,l,' Extracted box' check1=.FALSE. If (nvalid.EQ.4) check1=.TRUE. !write(*,*)i,j,k,l,' Is box valid?',nvalid,check1 If (check1) then call areacalc (area,xy) !write(*,*)i,j,k,l,' Area of box',area call emmcalc (emm,xy,index,a,rtype,mtype) !write(*,*)i,j,k,l,' Emmisivity weighting',emm counter=counter+1 gmask(i,j,k-1)=.TRUE. nphot(i,j,k)=nphot(i,j,k)+ & (gval*gval*gval*gval*area*emm) !write(*,*)i,j,k,l,' Number of photons',nphot(i,j,k) !write(16,*)i,j,k,l,gval,area,emm,nphot(i,j,k) End If end do !write(17,*)i,j,k,g(i,j,k-1),nphot(i,j,k),dg(i,j,k-1) !If (gmask(i,j,k-1)) !write(18,*)g(i,j,k-1),nphot(i,j,k) end do !write(18,*)'no no' end do End If end do checkout=.FALSE. If (counter.GT.0) checkout=.TRUE. !write(*,*)'Were valid points found in line profile?',checkout return end subroutine rebinp1 (nmask,n,e,de,emin,emax,gpmask,gmask, & nphot,dg,g,nb,ngp,ng,nc) implicit none integer nb,ngp,ng,nc double precision de(0:nb-1) double precision e(0:nb-1),n(nb) double precision nphot(ngp,nc,ng) double precision g(ngp,nc,0:ng-1) double precision dg(ngp,nc,0:ng-1) logical gpmask(ngp) logical nmask(0:nb-1) logical gmask(ngp,nc,0:ng-1) integer ng1 parameter (ng1=100) double precision emin,emax,deval,norm double precision eval(ng1),nval(ng1),x,y logical check1 logical mask1(ng1) integer i,j,k c emin=2.0d1 c emax=-2.0d1 c do i=1,ngp c If (gpmask(i)) then c do j=1,nc c emin=min(emin,g(i,j,0)) c emax=max(emax,g(i,j,99)) c end do c End If c end do !write(*,*)'Lower, upper limits of global grid',emin,emax e(0)=emin nmask(0)=.FALSE. deval=(emax-emin)/dble(nb-1) de(0)=deval do i=1,nb n(i)=0.0d0 If (i.LT.nb) then de(i)=deval nmask(i)=.FALSE. e(i)=e(i-1)+deval End If !write(19,*)i,e(i-1),n(i),de(i-1),nmask(i-1) end do !write(*,*)'Initialised global grids' norm = 0.0d0 do i = 1,ngp If (gpmask(i)) then do j=1,nc do k=1,ng mask1(k)=gmask(i,j,k-1) nval(k)=nphot(i,j,k)/dg(i,j,k-1) eval(k)=g(i,j,k-1)+0.5d0*(dg(i,j,k-1)) !write(20,*)i,j,k,eval(k),nval(k), ! & mask1(k),dg(i,j,k-1),gmask(i,j,k-1) end do !write(*,*)i,j,' Initialised interpolation grids' do k=1,nb check1=.FALSE. x=e(k-1)+0.5d0*deval call linterpolate (check1,y,x,mask1,eval,nval,ng) y = y * deval / x If (check1) nmask(k)=.TRUE. If (check1) n(k) = n(k) + y If (check1) norm = norm + y !write(21,*)i,j,k,x,y,n(k),check1 end do end do End If end do !write(*,*)'Rebinned lines onto global grid' do i = 1,nb n(k) = n(k) / norm !write(22,*)i,e(i-1),n(i) end do !write(*,*)'Normalised line' return end subroutine lsmooth (nmask,n,e,de,gpmask,gmask,g,nb,ngp,ng,nc) implicit none integer nb,ngp,ng,nc double precision de(0:nb-1) double precision e(0:nb-1),n(nb) double precision g(ngp,nc,0:ng-1) logical gpmask(ngp) logical nmask(0:nb-1) logical gmask(ngp,nc,0:ng-1) double precision d1,d2,dnde,tol double precision gbound(ngp,nc,2) integer i,j,k tol=1.0d-1 do i=1,ngp If (gpmask(i)) then do j=1,nc gbound(i,j,1)=g(i,j,0) gbound(i,j,2)=g(i,j,ng-1) !write(*,*)i,j,' Limits on line',gbound(i,j,1), ! & gbound(i,j,2) do k=1,nb-4 If (e(k).LE.gbound(i,j,1).AND. & e(k+1).GE.gbound(i,j,1)) then d1=dabs(e(k)-gbound(i,j,1)) d2=dabs(e(k+1)-gbound(i,j,1)) !write(*,*)i,j,k,1,e(k-1),n(k),nmask(k-1) !write(*,*)i,j,k,1,e(k),n(k+1),nmask(k),d1 !write(*,*)i,j,k,1,e(k+1),n(k+2),nmask(k+1),d2 If (d1.LT.d2) then nmask(k)=.FALSE. dnde=dabs((n(k+1)-n(k))/de(k-1)) If (dnde.LT.tol) nmask(k-1)=.FALSE. !write(*,*)i,j,k,1,1,dnde,de(k-1) Else nmask(k)=.FALSE. dnde=dabs((n(k+2)-n(k+1))/de(k+1)) If (dnde.LT.tol) nmask(k+1)=.FALSE. !write(*,*)i,j,k,1,2,dnde,de(k+1),nmask(k+1) End If End If If (e(k).LE.gbound(i,j,2).AND. & e(k+1).GE.gbound(i,j,2)) then d1=dabs(e(k)-gbound(i,j,2)) d2=dabs(e(k+1)-gbound(i,j,2)) !write(*,*)i,j,k,2,e(k-1),n(k),nmask(k-1) !write(*,*)i,j,k,2,e(k),n(k+1),nmask(k),d1 !write(*,*)i,j,k,2,e(k+1),n(k+2),nmask(k+1),d2 If (d1.LT.d2) then nmask(k)=.FALSE. dnde=dabs((n(k+1)-n(k))/de(k-1)) If (dnde.LT.tol) nmask(k-1)=.FALSE. !write(*,*)i,j,k,2,1,dnde,de(k-1) Else nmask(k)=.FALSE. dnde=dabs((n(k+2)-n(k+1))/de(k+1)) If (dnde.LT.tol) nmask(k+1)=.FALSE. !write(*,*)i,j,k,2,2,dnde,de(k+1),nmask(k+1) End If End If end do end do End If end do return end subroutine rebinp2 (ne,ec,ear,photar,nmask,n,e,de,nb) implicit none integer nb,ne real ear(0:ne),photar(ne) double precision ec,de(0:nb-1) double precision e(0:nb-1),n(nb) logical nmask(0:nb-1) double precision x,y,norm double precision eval(nb),nval(nb) logical check1 logical mask1(nb) logical mask1a(ne) integer i do i=1,nb mask1(i)=nmask(i-1) nval(i)=n(i)/de(i-1) eval(i)=ec*(e(i-1)+0.5d0*(de(i-1))) !write(24,*)i,eval(i),nval(i),mask1(i) end do !write(*,*)'Initialised interpolation grids' norm=0.0d0 do i=1,ne y = 0.0d0 check1=.FALSE. photar(i)=0.0e0 x=dble(ear(i-1)+0.5d0*(ear(i)-ear(i-1))) call linterpolate (check1,y,x,mask1,eval,nval,nb) y = y * (ear(i)-ear(i-1)) If (check1) norm = norm + y If (check1) photar(i) = real(y) mask1a(i)=check1 !write(25,*)i,x,y,photar(i),check1 end do !write(*,*)'Rebinned line onto XSPEC grid',norm do i=1,ne If (mask1a(i)) photar(i) = photar(i) / real(norm) c write(26,*)i,ear(i-1),photar(i) c write(27,*)ear(i-1),photar(i) end do !write(*,*)'Normalised line' c write(26,*) c write(27,*) return end subroutine areacalc (area,xy) implicit none double precision area,xy(2,2,4) double precision dx,dy,ds double precision a,b,c,d double precision ang1,ang2,s double precision zero,one,tol logical mask1(5) logical check1,check2 logical check4 integer nvalid zero=0.0d0 one=1.0d0 tol=1.0d-10 area=0.0d0 !write(12,*)xy(1,1,1),xy(1,1,2) !write(12,*)xy(1,2,1),xy(1,2,2) !write(12,*)xy(2,2,1),xy(2,2,2) !write(12,*)xy(2,1,1),xy(2,1,2) !write(12,*)xy(1,1,1),xy(1,1,2) !write(12,*)'no no' nvalid=0 mask1(1)=.FALSE. dx=dabs(xy(1,1,1)-xy(1,2,1)) dy=dabs(xy(1,1,2)-xy(1,2,2)) a=(dx*dx)+(dy*dy) If (a.GT.tol) then nvalid=nvalid+1 mask1(1)=.TRUE. a=dsqrt(a) End If !write(*,*)'Length of side a',a,dx,dy mask1(2)=.FALSE. dx=dabs(xy(1,1,1)-xy(2,1,1)) dy=dabs(xy(1,1,2)-xy(2,1,2)) b=(dx*dx)+(dy*dy) If (b.GT.tol) then nvalid=nvalid+1 mask1(2)=.TRUE. b=dsqrt(b) End If !write(*,*)'Length of side b',b,dx,dy mask1(3)=.FALSE. dx=dabs(xy(2,1,1)-xy(2,2,1)) dy=dabs(xy(2,1,2)-xy(2,2,2)) c=(dx*dx)+(dy*dy) If (c.GT.tol) then nvalid=nvalid+1 mask1(3)=.TRUE. c=dsqrt(c) End If !write(*,*)'Length of side c',c,dx,dy mask1(4)=.FALSE. dx=dabs(xy(1,2,1)-xy(2,2,1)) dy=dabs(xy(1,2,2)-xy(2,2,2)) d=(dx*dx)+(dy*dy) If (d.GT.tol) then nvalid=nvalid+1 mask1(4)=.TRUE. d=dsqrt(d) End If !write(*,*)'Length of side d',d,dx,dy mask1(5)=.FALSE. dx=dabs(xy(1,1,1)-xy(2,2,1)) dy=dabs(xy(1,1,2)-xy(2,2,2)) ds=(dx*dx)+(dy*dy) If (ds.GT.tol) then nvalid=nvalid+1 mask1(5)=.TRUE. ds=dsqrt(ds) End If !write(*,*)'Length of base of triangles',ds check1=.FALSE. If (nvalid.EQ.5) check1=.TRUE. !write(*,*)'Are all lengths valid?',nvalid,check1 !write(14,*)a,b,c,d,ds,nvalid,check1 If (check1) then check2=.FALSE. ang1=(a*a)+(d*d)-(ds*ds) ang1=ang1/((a*d)+(a*d)) If (dabs(ang1).LE.one) check2=.TRUE. If (check2) ang1=dacos(ang1) !write(*,*)i,j,'Value of first angle',ang1 check4=.FALSE. ang2=(b*b)+(c*c)-(ds*ds) ang2=ang2/((b*c)+(b*c)) If (dabs(ang2).LE.one) check4=.TRUE. If (check4) ang2=dacos(ang2) !write(*,*)i,j,'Value of second angle',ang2 If (check2.AND.check4) then s=a+b+c+d s=0.5d0*s !write(*,*)i,j,'Length of semiperimeter',s !write(19,*)ang1,ang2,s area=dcos(0.5d0*(ang1+ang2)) area=area*area area=a*b*c*d*area area=((s-a)*(s-b)*(s-c)*(s-d))-area area=dsqrt(area) !write(*,*)i,j,'Area of trapezoid',area !write(20,*)area End If Else nvalid=0 c do i=1,4 c If (mask1(i)) nvalid=nvalid+1 c end do !write(*,*)'Number of valid edges',nvalid check2=.FALSE. If (nvalid.EQ.3) check2=.TRUE. !write(*,*)'Is shape triangular?',check2 If (check2) then If (.NOT.(mask1(1))) then a=b b=c c=d Else if (.NOT.(mask1(2))) then b=c c=d Else if (.NOT.(mask1(3))) then c=d End If !write(*,*)'Sides of triangle',a,b,c area=(a+b+c)*(b+c-a) area=area*(c+a-b)*(a+b-c) area=dsqrt(area) area=0.25d0*area !write(*,*)'Area of triangles',area End If End If return end subroutine emmcalc (emm,xy,id,a,rt,mt) implicit none integer rt,mt double precision emm,xy(2,2,4),a,id double precision re,me,fztibc integer i,j re=0.0d0 me=0.0d0 do i=1,2 do j=1,2 re=re+xy(i,j,3) me=me+xy(i,j,4) end do end do re = 0.25d0 * re me = 0.25d0 * me !write(*,*)'Averaged values of re, me',re,me emm=1.0d0 If (rt.EQ.1) then emm = emm * (re**(-id)) Else if (rt.EQ.2) then emm = emm * fztibc(re,a) End If If (mt.EQ.1) then emm = emm / me Else if (mt.EQ.2) then emm = emm * ((1.0d0+(2.06d0*me))/2.06d0) Else if (mt.EQ.3) then emm = emm * (me*dlog(1.0d0+(1.0d0/me))) End If !write(*,*)'Emissivity',emm !write(15,*)re,me,id,emm return end subroutine pathinit (path,rmin,rmax,a,obs) double precision rmin,rmax,a,obs character*256 path,fgmodf character*21 stringp1 character*36 stringp2 character*4 s1 character*3 s2 character*9 fmt integer x(4),zero integer pathlength integer lenact external fgmodf,lenact zero=0 stringp1='0000000000000.FITS.gz' !write(*,*)'Input variables',a,obs,rmin,rmax !write(*,*)'Initialised stringp1 ',stringp1 x(1)=nint(a*1.0d3) x(2)=nint(obs) x(3)=nint(rmin) x(4)=nint(rmax) !write(*,*)'Setup values',x fmt='(I4)' s1 = ' ' write(s1,fmt) x(1) stringp1(1:4) = s1 !write(*,*)'First component of file name ',stringp1 fmt='(I3)' s2 = ' ' write(s2,fmt) x(2) stringp1(5:7) = s2 !write(*,*)'Second component of file name ',stringp1 fmt='(I3)' s2 = ' ' write(s2,fmt) x(3) stringp1(8:10) = s2 !write(*,*)'Third component of file name ',stringp1 fmt='(I3)' s2 = ' ' write(s2,fmt) x(4) stringp1(11:13) = s2 !write(*,*)'Fourth component of file name ',stringp1 fmt='(I1)' If (x(1).LT.1000) write(stringp1(1:1),fmt) zero If (x(1).LT.100) write(stringp1(2:2),fmt) zero If (x(1).LT.10) write(stringp1(3:3),fmt) zero If (x(2).LT.100) write(stringp1(5:5),fmt) zero If (x(2).LT.10) write(stringp1(6:6),fmt) zero If (x(3).LT.100) write(stringp1(8:8),fmt) zero If (x(3).LT.10) write(stringp1(9:9),fmt) zero If (x(4).LT.100) write(stringp1(11:11),fmt) zero If (x(4).LT.10) write(stringp1(12:12),fmt) zero !write(*,*)'Correct filename ',stringp1 stringp2(1:14)='XSKDLINE/0000/' stringp2(10:13) = s1 If (x(1).LT.1000) write(stringp2(10:10),fmt) zero If (x(1).LT.100) write(stringp2(11:11),fmt) zero If (x(1).LT.10) write(stringp2(12:12),fmt) zero stringp2(15:36)=stringp1 c write(*,*)'Path to file ',stringp2,stringp1 path = fgmodf() c write(*,*)'Path to Data directory',datadir pathlength = lenact(path) c write(*,*)'Length of path' path(pathlength+1:pathlength+37) = stringp2 c write(*,*)'Full path to file',path return end subroutine tfopen (checkout,mask,g,r,m,x,y,path,ng,nr,nc) implicit none integer ng,nr,nc character*256 path double precision g(ng+1,nr+1,nc) double precision r(ng+1,nr+1,nc) double precision m(ng+1,nr+1,nc) double precision x(ng+1,nr+1,nc) double precision y(ng+1,nr+1,nc) logical checkout,mask(ng+1,nr+1,nc) integer n parameter (n=20402) double precision nullval double precision a(n),b(n),c(n),d(n),e(n) logical m1(n),check1 character errtext*30,errmessage*80 integer unit,status,readwrite integer blocksize,hdutype integer frow,felem,colnum integer i,j,k,l status=0 call ftgiou(unit,status) c write(*,*)'Got unit number',unit,status readwrite=0 call ftopen(unit,path,readwrite,blocksize,status) c write(*,*)'Opened file, status',status call ftmahd(unit,2,hdutype,status) c write(*,*)'Moved to binary table',hdutype,status frow=1 felem=1 colnum=1 nullval=0.0d0 check1=.FALSE. call ftgcvd (unit,colnum,frow,felem,n,nullval,a,check1,status) c write(*,*)'Read in redshifts, result',check1,status colnum=2 check1=.FALSE. call ftgcvd (unit,colnum,frow,felem,n,nullval,b,check1,status) c write(*,*)'Read in radii, result',check1,status colnum=3 check1=.FALSE. call ftgcvd (unit,colnum,frow,felem,n,nullval,c,check1,status) c write(*,*)'Read in angles, result',check1,status colnum=4 check1=.FALSE. call ftgcvd (unit,colnum,frow,felem,n,nullval,d,check1,status) c write(*,*)'Read in alpha coords, result',check1,status colnum=5 check1=.FALSE. call ftgcvd (unit,colnum,frow,felem,n,nullval,e,check1,status) c write(*,*)'Read in beta coords, result',check1,status colnum=6 check1=.FALSE. call ftgcl (unit,colnum,frow,felem,n,m1,check1,status) c write(*,*)'Read in mask, result',check1,status call ftclos(unit, status) call ftfiou(unit, status) If (status.GT.0) then call ftgerr(status,errtext) write(*,*)'FITSIO Error status =',status,': ',errtext call ftgmsg(errmessage) do while (errmessage.NE. ' ') write(*,*)errmessage call ftgmsg(errmessage) end do End If c write(*,*)'Status of FITS read',status checkout=.FALSE. If (status.EQ.0) checkout=.TRUE. c write(*,*)'Was FITS read succesfull?',checkout If (checkout) then i=1 do j=1,ng+1 do k=1,nr+1 do l=1,nc g(j,k,l)=a(i) r(j,k,l)=b(i) m(j,k,l)=c(i) x(j,k,l)=d(i) y(j,k,l)=e(i) mask(j,k,l)=m1(i) i=i+1 end do end do end do c write(*,*)'Copied TF to output arrays' End If return end subroutine gridinit (x,xmin,xmax,scale,n) implicit none integer n double precision x(n),xmin,xmax,scale double precision step integer i x(1)=xmin x(n)=xmax step=(xmax-xmin)*scale do i=2,n-1 x(i)=xmin+(dble(i-1)*step) end do return end double precision function rms (a) implicit none double precision a double precision a2,z1,z2 double precision one,three,third one=1.0d0 three=one+one+one third=one/three a2=a*a !write(*,*)'Setup parameters',one,three,third,a,a2 z1=((one+a)**third)+((one-a)**third) z1=((one-a2)**third)*z1 z1=z1+one !write(*,*)'Value of Z1',z1 z2=a2+a2+a2+(z1*z1) z2=dsqrt(z2) !write(*,*)'Value of Z2',z2 rms=(three-z1)*(three+z1+z2+z2) rms=dsqrt(rms) rms=three+z2-rms !write(*,*)'Location of rms',rms return end subroutine linterpolate (checkout,ytarg,xtarg,xymask,xin,yin,n) implicit none integer n double precision xtarg,ytarg double precision xin(n),yin(n) logical xymask(n),checkout double precision x(n),y(n),frac integer nvalid integer i,j i=1 nvalid=0 do j=1,n x(j)=0.0d0 y(j)=0.0d0 If (xymask(j)) then x(i)=xin(j) y(i)=yin(j) nvalid=nvalid+1 i=i+1 End If end do !write(*,*)'Sorted points into ascending order in x' !write(*,*)'Number of valid points',nvalid ytarg=0.0d0 do i=1,nvalid-1 If ((x(i).LE.xtarg).AND.(x(i+1).GE.xtarg)) then checkout=.TRUE. frac=(xtarg-x(i))/(x(i+1)-x(i)) ytarg=y(i)+frac*(y(i+1)-y(i)) exit End If end do !write(*,*)'Was valid point found?',xtarg,ytarg,checkout return end double precision function fztibc (r,a) implicit none double precision r,a double precision rms,x,pi double precision x0,x1,x2,x3 double precision dm1,dm2,dm3 double precision dm4,dm5,dm6 double precision third pi = dacos(-1.0d0) third = 1.0d0/3.0d0 x = dsqrt(r) x0 = dsqrt(rms(a)) x1 = 2.0d0*dcos(third*(dacos(a) - pi)) x2 = 2.0d0*dcos(third*(dacos(a) + pi)) x3 = -2.0d0*dcos(third*dacos(a)) dm1 = 0.5d0*3.0d0*a*dlog(x/x0) dm2 = dlog((x-x1)/(x0-x1)) dm2 = 3.0d0*(x1-a)*(x1-a)*dm2 dm2 = dm2 / (x1*(x1-x2)*(x1-x3)) dm3 = dlog((x-x2)/(x0-x2)) dm3 = 3.0d0*(x2-a)*(x2-a)*dm3 dm3 = dm3 / (x2*(x2-x1)*(x2-x3)) dm4 = dlog((x-x3)/(x0-x3)) dm4 = 3.0d0*(x3-a)*(x3-a)*dm4 dm4 = dm4 / (x3*(x3-x1)*(x3-x2)) dm5 = x - x0 - dm1 - dm2 - dm3 - dm4 dm6 = (x*x*x) - x - x - x + a + a dm6 = x * x * dm6 dm6 = 0.5d0 * 3.0d0 *dm5 / dm6 fztibc = dm6 / r end function