c Calculates (b, phi, Amp) table for passing photons from the sky. c Copyright 1999 by Robert J. Nemiroff. c May be used non-commercially with proper acknowledgement. program bhsky_table implicit real*8 (a-h,o-z) character*3 converge c open (1, file='star.dat', status='old') open (2, file='sky.out', status='unknown') c pi=dacos(-1.0d0) iter=1000 eps=1.0d-4 c read (1,*) read (1,10) Rs read (1,10) Rad read (1,10) Dobs 10 format (20x, f10.5) close (1) c 0 degrees is toward the star. Rinf=Rad/sqrt(1.0d0-Rs/Rad) bmin=3.0d0*dsqrt(3.0d0)/2.0d0*Rs Bbig=Dobs/sqrt(1.0d0-Rs/Dobs) sinarg=3.0d0*sqrt(3.0d0)*Rs/2.0d0/Dobs*sqrt(1.0d0-Rs/Dobs) if (sinarg.gt.1.0d0) sinarg=1.0d0 delesc=dasin(sinarg) if (Rad.gt.1.5d0*Rs) then Earg=Rad/Dobs*sqrt((1.0d0-Rs/Dobs)/(1.0d0-Rs/Rad)) else Earg=1.5d0*Rs/Dobs*sqrt((1.0d0-Rs/Dobs)/ 1 (1.0d0-Rs/(1.5d0*Rs))) endif delearth=dasin(Earg) if (Dobs.lt.1.5d0*Rs) then delesc=pi-delesc delearth=pi-delearth endif write (*,*) ' delesc = ', delesc*180.0d0/pi write (*,*) ' delearth = ', delearth*180.0d0/pi converge='off' c ddel is negative. dphi=pi/180.0d0 delmax=pi delmin=delesc ddelbig=(delmin-delmax)/179.4754323247d0 ddelsmall=ddelbig/10.0d0 ddeltiny=ddelbig/50.0d0 deltest=delmax - ddelbig c do 500 id=1800, 0, -1 if (phideg.le.175.0d0) ddel=ddelbig if (phideg.gt.175.0d0.and.phideg.le.185.0d0) 1 ddel=ddelsmall if (phideg.gt.185.0d0.and.phideg.le.330.0d0) 1 ddel=ddelbig if (phideg.gt.300.0d0) ddel=ddeltiny if (dphi.ge.5.0d0*pi/180.0d0) ddel=ddel/2.0d0 deltest=deltest + ddel btest=Bbig*sin(deltest) if (b.eq.0.0d0) b=bmin/1000.0d0 c Star outside photon sphere. 20 if (Rad.ge.1.5d0*Rs) then c Infinity to observer. if (deltest.gt.pi-delesc) then b=btest delnew=deltest s=Rad eobs=dacos(Rad/Dobs) emin=eobs emax=pi/2.0d0 goto 40 endif c Infinity - rtop - observer. if (deltest.gt.delearth.and.deltest.le.pi-delesc) then b=btest delnew=deltest rtop=2.0d0*b/sqrt(3.0d0) 1 *cos(acos(-sqrt(27.0d0)*Rs/2.0d0/b)/3.0d0) s=rtop eobs=dacos(s/Dobs) emin=0.0d0 emax=pi/2.0d0 goto 40 endif c Photon would hit earth: converge. if (deltest.le.delearth) then if (phideg.gt.360.0d0) goto 999 if (dphi.lt.eps) goto 999 converg=1.2d0 delnew=(delold*delearth**(converg-1.0d0))**(1.0d0/converg) b=Bbig*sin(delnew) rtop=2.0d0*b/sqrt(3.0d0) 1 *cos(acos(-sqrt(27.0d0)*Rs/2.0d0/b)/3.0d0) s=rtop eobs=dacos(s/Dobs) emin=0.0d0 emax=pi/2.0d0 converge='on' goto 40 endif c endif c Star inside photon sphere. if (Rad.lt.1.5d0*Rs) then c Observer outside photon sphere. if (Dobs.gt.1.5d0*Rs) then c Infinity - observer. if (deltest.gt.pi-delesc) then if (phideg.gt.360.0d0) goto 999 b=btest delnew=deltest s=Rad eobs=dacos(s/Dobs) emin=eobs emax=pi/2.0d0 goto 40 endif c Infinity - rtop - observer if (deltest.ge.delearth.and.deltest.lt.pi-delesc) then if (phideg.gt.360.0d0) goto 999 b=btest delnew=deltest rtop=2.0d0*b/sqrt(3.0d0) 1 *cos(acos(-sqrt(27.0d0)*Rs/2.0d0/b)/3.0d0) s=rtop eobs=dacos(s/Dobs) emin=0.0d0 emax=pi/2.0d0 goto 40 endif c Photon would hit earth: converge. if (deltest.lt.delearth) then if (phideg.gt.360.0d0) goto 999 if (dphi.lt.eps) goto 999 converg=1.01d0 delnew=(delold*delearth**(converg-1.0d0))**(1.0d0/converg) b=Bbig*sin(delnew) rtop=2.0d0*b/sqrt(3.0d0) 1 *cos(acos(-sqrt(27.0d0)*Rs/2.0d0/b)/3.0d0) s=rtop eobs=dacos(s/Dobs) emin=0.0d0 emax=pi/2.0d0 converge='on' goto 40 endif c endif c Observer inside photon sphere. if (Dobs.le.1.5d0*Rs) then c Infinity - observer. if (deltest.gt.delesc.and.phideg.lt.300.0d0) then b=btest delnew=deltest s=Rad eobs=dacos(s/Dobs) emin=eobs emax=pi/2.0d0 goto 40 endif c No rtop possible. Converge. if (deltest.lt.delesc.or.phideg.ge.300.0d0) then if (phideg.gt.360.0d0) goto 999 if (dphi.lt.eps) goto 999 converg=1.01d0 delnew=(delold*delesc**(converg-1.0d0))**(1.0d0/converg) b=Bbig*sin(delnew) s=Rad eobs=dacos(s/Dobs) emin=eobs emax=pi/2.0d0 converge='on' goto 40 endif c endif c endif c write (*,*) ' You shouldnt be here! ' c 40 sum=0.0d0 if (rtop.lt.Rad) rtop=Rad if (b.eq.0.0d0) b=bmin/1000.0d0 de=(emax-emin)/dble(iter) e=emin-de/2.0d0 c do 100 i=1, iter - 1 e=e+de c Integrate over ranges not near zero. if (e.lt.pi/4.0d0) then top=tan(e)*de bot=s**2/b**2/cos(e)**2 - 1.0d0 + Rs*cos(e)/s else top=de bot=s**2/b**2/sin(e)**2 - cos(e)**2/sin(e)**2 1 + cos(e)**3/sin(e)**2*Rs/s endif c bot=sqrt(bot) term=top/bot c How do I add thee: one of three ways. c 1) No rtop (collision trajectery away). if (delnew.gt.pi-delesc) then sum=sum + term goto 100 endif c 2) rtop exists and photon passes it - c so count this piece twice. if (e.lt.eobs.and.delnew.lt.pi/2.0d0) then sum=sum+2.0d0*term goto 100 endif c 3) rtop exists and photon never passes it - c No deflection this time through - no c sum statement tripped! if (e.gt.eobs) sum=sum+term c 100 continue c phiold=phinew phinew=sum phideg=phinew*180.0d0/pi phiout=180.0d0-phideg dphi=phinew-phiold c delout=delnew*180.0d0/pi ddel=delnew-delold delold=delnew c if (delnew*180.0d0/pi.ge.179.999d0) then Amp=1.0d0 goto 425 endif c top=ddel*dacos(cos(delnew)**2 - sin(delnew)**2) c bot=dphi*dacos(cos(pi-phinew)**2 - sin(pi-phinew)**2) top=ddel*dsin(delnew) bot=dphi*dsin(pi-phinew) blueshift=1.0d0 Amp=abs(top/bot)*blueshift if (phiold.eq.0.0d0) Amp=1.0d0 if (Amp.ge.100.0d0) Amp=100.0d0 c 425 write (*,450) phideg, delout, Amp write (2,450) phideg, delout, Amp 450 format (2x, f10.5, 5x, f10.5, 5x, g12.5) c if (converge.eq.'on') goto 20 500 continue c 999 stop end