program tranlow * * This program applies the transfer function of a rotating black hole * to a given local disc spectrum. * The transfer function is read from fort.1, the local fisk spectrum is read * from fort.2. * The output (spectrum at infinity, % polarization, polarization plane) * is printed in fort.4 * The input spectrum is given as a function of radius and frequency, assuming * an electron scattering inclination dependence of the local disk emission. * The local polarization amplitude is that of electron scattering scaled by the * ratio of electron scattering/total opacity (varies with radius and frequency). * The observed spectrum, and degree and angle of polarization are given as a * function of frequency for the specified polarization. * * The following parameters are used for the transfer function: * ne= # of energy intervals, nrd= # radii intervals, nan= # intervals in * cos(theta) * rbin= the radii values, ebin= energy, anbin= cos(theta), trs= the transfer * matrix elements, all are read from fort.1 * parameter (ne=30, nrd=35, nan=30, nq=ne*nrd*nan*3) * * The following parameters are used for the local disk spectrum: * r=the disk model radius, spec=the flux, freq= the frequency * abs=electron scattering/(es+absorption), it is the relative part of the * scattering in the total opacity (a function of frequency and radius). * all are read from tape2 * real rbin(nrd),ebin(ne),anbin(nan),trs(nrd,nan,ne,3) & ,abs(50,100),spec(50,100),freq(100),r(51),flux(500,3) & ,spl(500,3),spg(500,3),spt(500,3),chan(21) & ,freqh(500),flnewt(500),spnewt(500),acu(50,10) character*8 ftime data trs/nq*0/ data spt,flux /3000*0/ * * chan=the polarization from chandrasekhar * data chan/ .11713, .08979, .07448,.06311,.0541,.04667,.04041 1 ,.03502, .03033, .02619, .02252, .01923, .01627, .01358, .01112 2 ,.00888, .00682, .00492, .00315, .00152, .0 / * * reading the local spectrum * * note: the local spectrum file (fort.2) should start with the following: * obs=cosine(theta), obs=1. for a face on disc. * em9=the bh mass in units of 10e9 solar masses * r=radius/(g*mass/c**2), i.e. radius in gravitational units (=1.476e14, em9=1) * the local spectrum is given as flux (ergs/sec/hz/cm**2) at 100 frequncies * at 50 rings, the 51 ring boundaries are in r(i). note that the frequencies * must be spaced in equal logarithmic steps (this was assumed when the * transfer function was created). * reads transfer matrix from fort.1, * call time(ftime) print*,'program started on ', ftime irot=1 nfr=100 print*,' type in cos inclination from 1 to 0.1' print*,' 1 = face on, 0.1 = nearly edge on' read*, obs read(2,*) oo1,em9 print*, 'BH mass in units of E9 M_sun=',em9 read(2,*) (r(i),i=1,51) print*,'outer disk radius in units of R_g=',r(51) read(2,*) (freq(i),i=1,100) read(2,*) ((abs(i,j),spec(i,j),j=1,100),i=1,50) 104 format(10g13.5) * * reading the transfer function * call time(ftime) write(*,*)'disk spectrum read by ', ftime read(1,41) (rbin(i),i=1,nrd) read(1,41) (anbin(i),i=1,nan) read(1,41) (ebin(i),i=1,ne) 41 format( 8f9.4) 42 format( 8e9.3) do 1 ir=1,nrd do 1 ip=1,3 read(1,*) nr,np if(nr.ne.ir.or.np.ne.ip) stop 'error in reading data' do 1 it=1,nan read(1,*) n1,n2 if(n1.eq.0) goto 1 read(1,42) (trs(ir,it,ie,ip),ie=n1,n2) * * aw=angular width of an inclination bin. the flux accumulated in this * ring is diveded by aw to get the apparent emitted flux (i.e. assuming * the source is emitting isotropically) * if(it.eq.1) then aw=(anbin(2)-anbin(1))/2. elseif(it.eq.nan)then aw=(anbin(nan)-anbin(nan-1))/2. else aw=(anbin(it+1)-anbin(it-1))/2. endif do 2 ie=n1,n2 2 trs(ir,it,ie,ip)=trs(ir,it,ie,ip)/aw 1 continue call time(ftime) write(*,*)'transfer function read ', ftime * * dt1, dt2 = the factors in the interpolation of trs for a given obs * warning!! the lowest possible inclination is myu=0.96, i.e. obs=1 * is actually using myu=0.96, which is the average myu for the lowest * inclination bin * np=3 do 10 i=1,nan-1 dt1=(obs-anbin(i))/(anbin(i+1)-anbin(i)) dt2=(anbin(i+1)-obs)/(anbin(i+1)-anbin(i)) itt=i if(dt1*dt2.ge.0) goto 9 10 continue * * pol-the chandrasekhar polarization for the given obs * 9 do 8 ip=1,20 if(obs.gt.(ip-1)*.05.and.obs.le.ip*.05)then dpol=(obs-(ip-1)*.05)/.05 pol=dpol*chan(ip+1)+(1-dpol)*chan(ip) goto 11 endif 8 continue print*,' obs is out of range obs= ',obs * * begining of transfer application * summing over all radii of the model * first: find the interpolation in radius, i.e. disk model radius (r) vs. transfer function radius (rbin). * 11 do 12 ir=1,50 rm=(r(ir)+r(ir+1))/2. if(rm.le.6.and.irot.eq.0) goto 12 itr=0 do 13 irt=1,nrd-1 dr1=-(rm-rbin(irt))/(rbin(irt)-rbin(irt+1)) dr2=-(rbin(irt+1)-rm)/(rbin(irt)-rbin(irt+1)) if(dr1*dr2.ge.0) then itr=irt goto 14 endif 13 continue if(itr.eq.0) then if(rm.gt.rbin(1)) goto 99 if(rm.lt.rbin(nrd))then dr1=1. dr2=0. itr=nrd-1 else stop ' a bug in interpolating the radius' endif endif * * rbin(nrd) is the inner radius of the transfer function, if rm is smaller * use rbin(nrd). rbin(1) is the outer radius of the transfer function, if * rm is larger then skip the transfer function by going to 99 * * * sp-loads the flux emitted at rm. the loop runs over all energy shifts, for * each shift it runs over 3 polarization components (i,q,u) * 14 do 16 ie=1,ne e=ebin(ie) do 16 ip=1,np * * interpolating trs between 2 angles and between 2 radii * d1=dt1*dr1 d2=dt1*dr2 d3=dt2*dr1 d4=dt2*dr2 xn=trs(itr,itt,ie,ip)*d4+trs(itr+1,itt,ie,ip)*d3 & +trs(itr,itt+1,ie,ip)*d2 + trs(itr+1,itt+1,ie,ip)*d1 if(xn.eq.0) goto 16 * * spl-loads the flux of the three components after the trs transformation * do 17 if=1,nfr spl(if,ip)=spec(ir,if)*xn if(ip.ge.2) spl(if,ip)=spl(if,ip)*abs(ir,if) 17 continue adfr=alog(freq(2)/freq(1)) a=alog(e)/adfr if(a.ge.0)then ia=ifix(a) a1=a-ia else a=-a ia=-ifix(a)-1 a1=-ia-a endif * * ia-the no. of bins shift in frequency corresponding to the energy shift * a1-the fraction of bin shift * spg-the flux shifted in frequency, spt-accumulates shifted flux for all * energy shifts * do 18 if=2,nfr if1=if+ia if(if1.gt.nfr.or.if1.lt.1) goto 18 spg(if1,ip)=spl(if,ip)*(1.-a1)+spl(if-1,ip)*a1 spt(if1,ip)=spg(if1,ip)+spt(if1,ip) 18 continue 16 continue 99 area=(r(ir+1)**2-r(ir)**2) do 20 if=1,nfr * * if r is out of the transfer matrix then we just correct the flux to the * limb darkening*myu, and put in the chndra' pol', since the angle is 90 deg * q is in the 180 direction, i.e. negative. * spnewt- the newtonian (disk rest frame) spectrum c cc spnewt(if)=spec(ir,if)*(.843*obs+1.7355*obs**2) spnewt(if)=spec(ir,if) if(rm.gt.rbin(1)) then spt(if,1) =spec(ir,if)*(.843*obs+1.7355*obs**2) spt(if,2)= -spt(if,1)*pol*abs(ir,if) spt(if,3)=0. endif flnewt(if)=flnewt(if)+spnewt(if)*area*2 do 20 ip=1,np * * multiplying by the area*2 for both sides of the disk, i.e. assuming * isotropic emission. flux-accumulates the flux from all rings. * flnewt-the disk spectrum with no relativistic effects. c flux(if,ip)=flux(if,ip)+spt(if,ip)*area*2. if10=if/10 if(10*if10.eq.if) acu(ir,if10)=flux(if,1) spt(if,ip)=0. 20 continue 12 continue do 21 if=1,nfr * * s1-the observed flux, s2-degree of polarization, s3-angle of polarization c xlm=3.e18/freqh(if) arlog=alog10(6.8446e28*em9*em9) s1=alog10(flux(if,1))+arlog s1newt=alog10(flnewt(if))+arlog alf=alog10(freq(if)) s2=sqrt(flux(if,2)**2+flux(if,3)**2)/flux(if,1) if(flux(if,2).eq.0) flux(if,2)=1.e-8 s3=57.29*.5*atan(flux(if,3)/flux(if,2)) if(s3.lt.-20.) s3=s3+90 c-- alf = log frequency (Hz), s1 = log flux (erg/sec/Hz), s2 = polarized fraction c-- s3 = polarization plane write(4,'(4e12.5)') alf,s1+alf,s2,s3 cc write(4,'(3e12.5)') alf,s1,s1newt ! print fluxes only cc write(4,'(2e12.5)') alf,s1newt ! print fluxes only 21 continue call time(ftime) write(*,*)'end reached on ', ftime write(*,*)'nu, nu*f_nu, %P, P plane, printed in unit 4 (fort.4)' end