c Near-field lensing image stars background to bh in file simage.in. c Copyright 1999 by Robert J. Nemiroff. c May be used non-commercially with proper acknowledgement. program bhsky_image implicit real*8 (a-h,o-z) parameter (ip=10000) dimension x1(ip), y1(ip), x2(ip), y2(ip) dimension A1(ip), A2(ip), color1(ip), color2(ip) dimension aread(ip), dread(ip), sread(ip), bread(ip) common /int/ thint(ip), delint(ip), Aint(ip) c Define files. open (1, file='bhsky1.out', status='unknown') open (2, file='bhsky2.out', status='unknown') open (7, file='sky.out', status='old') open (8, file='star.dat', status='old') open (9, file='star.num', status='old') c Read in bh parameters. write (*,*) ' Before read' read (8,*) read (8,2) Rs read (8,2) Rad read (8,2) Dobs read (8,2) smagcut read (8,2) CthetaRot read (8,2) CphiRot read (8,2) SthetaRot read (8,2) SphiRot close (8) write (*,*) ' Rs (km) = ', Rs write (*,*) ' Rad (km) = ', Rad write (*,*) ' Dobs (km) = ', Dobs write (*,*) ' smagcut = ', smagcut write (*,*) ' CthetaRot = ', CthetaRot write (*,*) ' CphiRot = ', CphiRot write (*,*) ' SthetaRot = ', SthetaRot write (*,*) ' SphiRot = ', SphiRot 2 format (20x, f10.5) c Rotate to put Orion at 0. SphiRot=-SphiRot+90.0d0 c Rot to put Orion just off center. SphiRot = SphiRot - 10.0 c Rot to Hamburg latitude. c SthetaRot=SthetaRot + 53.51d0 c Read in th, delta data table. pi=dacos(-1.0d0) f=pi/180.0d0 zblue=sqrt(1.0d0-Rs/Dobs) CthetaRot=CthetaRot*pi/180.0d0 CphiRot=CphiRot*pi/180.0d0 SthetaRot=SthetaRot*pi/180.0d0 write (*,*) ' zblue = ', zblue 5 format (1h+,i8) c do 10, i=1, 100000 read (7,*, end=20) thint(i), delint(i), Aint(i) 10 continue c 20 nth=i-1 close (7) c Read in star stuff. write (*,*) ' Reading in star data . . .' c do 50 i=1, ip read (9,*,end=75) aread(i), dread(i), sread(i), bread(i) 50 continue c 75 nptot=i-1 close (9) c npt1=0 npt2=0 camnorm=sqrt(1.0d0-cos(pi/4.0d0)) amagcut=10.0d0**(-smagcut/2.512d0+2.0d0) c earth=asin(Rad/Dobs)*180.0d0/pi Rearthapp=Rad*sqrt((1.0d0-Rs/Dobs)/(1.0d0-Rs/Rad)) earth=atan2(Rearthapp,Dobs)*180.0d0/pi write (*,*) ' Earth angle subtended = ', earth c z is observer direction. write (*,*) ' <=== Processing point #' c START STAR LOOP. do 100 i=1, nptot if (100*(i/100).eq.i) write (*,5) i c Iterate unlensed star image. c Rotate sky in phi. alpha=aread(i)+SphiRot delta=dread(i) smag=sread(i) bv=bread(i) c xin0=-sin(delta*f+pi/2.0d0)*cos(alpha*f+pi/2.0d0) zin0=sin(delta*f+pi/2.0d0)*sin(alpha*f+pi/2.0d0) yin0=-cos(delta*f+pi/2.0d0) c Rotate sky in theta. xin=xin0 yin=yin0*cos(SthetaRot) - zin0*sin(SthetaRot) zin=yin0*sin(SthetaRot) + zin0*cos(SthetaRot) c psize=10.0d0**(-smag/2.512d0+2.0d0) r3=sqrt(xin**2 + yin**2 + zin**2) r2=sqrt(xin**2 + yin**2) c Theta for first image. th1=acos(zin/r3)*180.0d0/pi c if (th1.lt.110.0d0) goto 70 c Find tabled del for given theta. call interpD(th1, del1, nth) if (del1.lt.earth) goto 70 call interpA(th1, Amp1, nth) c Assume no pinhole dimming. apixel=sqrt(psize*Amp1) if (apixel.lt.amagcut) goto 70 c Define new spherical coords. thin=del1*pi/180.0d0 phin=atan2(yin,xin) c Rotate Camera in phi. phin=phin + CphiRot c Transform to cartesian coords. xin_orig=sin(thin)*cos(phin) yin_orig=sin(thin)*sin(phin) zin_orig=cos(thin) c Rotate Camera in theta. xin_view=xin_orig yin_view=yin_orig*cos(CthetaRot) - zin_orig*sin(CthetaRot) zin_view=yin_orig*sin(CthetaRot) + zin_orig*cos(CthetaRot) if (zin_view.lt.0.0d0) goto 70 r2in_view=sqrt(xin_view**2 + yin_view**2) c Back to spherical coords. thin_view=dacos(zin_view) phin_view=atan2(yin_view, xin_view) c thin_view outside photo? thcomp=thin_view*180.0d0/pi if (thcomp.gt.70.0d0) goto 70 c npt1=npt1+1 c Flat-field Camera Formalism. x1(npt1)=sqrt(1.0d0-cos(thin_view))/camnorm 1 *(xin_view/r2in_view) y1(npt1)=sqrt(1.0d0-cos(thin_view))/camnorm 1 *(yin_view/r2in_view) c Again test: outside photo? if (abs(x1(npt1)).gt.1.0d0.or.abs(y1(npt1)).gt.1.0d0) then npt1=npt1-1 goto 70 endif c A goes as image radius. A1(npt1)=apixel if (apixel.lt.0.54d0) then color1(npt1)=(apixel**2/0.54d0)*100.0d0 + 200.0d0 if (color1(npt1).gt.254.0d0) color1(npt1)=254.0d0 if (color1(npt1).lt.201.0d0) color1(npt1)=201.0d0 else color1(npt1)=(bv+0.25d0)/2.0d0*100.0d0 *zblue**6 if (color1(npt1).gt.100.0d0) color1(npt1)=100.0d0 if (color1(npt1).lt.1.0d0) color1(npt1)=1.0d0 endif c Theta for second image. 70 th2=360.0d0-th1 call interpD(th2, del2, nth) if (del2.lt.earth) goto 100 call interpA(th2, Amp2, nth) apixel=sqrt(psize*Amp2) if (apixel.lt.amagcut.and.smag.gt.smagcut) goto 100 c Define new spherical coords. thin=del2*pi/180.0d0 phin=atan2(yin,xin) c Rotate Camera in phi. phin=phin + CphiRot c Transform to cartesian coords. xin_orig=sin(thin)*cos(phin) yin_orig=sin(thin)*sin(phin) zin_orig=cos(thin) c Invert image. xin_orig=-xin_orig yin_orig=-yin_orig c Rotate camera in theta. xin_view=xin_orig yin_view=yin_orig*cos(CthetaRot) - zin_orig*sin(CthetaRot) zin_view=yin_orig*sin(CthetaRot) + zin_orig*cos(CthetaRot) if (zin_view.lt.0.0d0) goto 100 r2in_view=sqrt(xin_view**2 + yin_view**2) c Back to spherical coords. thin_view=dacos(zin_view) phin_view=atan2(yin_view, xin_view) c thin_view outside photo? thcomp=thin_view*180.0d0/pi if (thin_view.gt.70.0d0) goto 100 c npt2=npt2+1 x2(npt2)=sqrt(1.0d0-cos(thin_view))/camnorm 1 *(xin_view/r2in_view) y2(npt2)=sqrt(1.0d0-cos(thin_view))/camnorm 1 *(yin_view/r2in_view) A2(npt2)=apixel if (abs(x2(npt2)).gt.1.0d0.or.abs(y2(npt2)).gt.1.0d0) then npt2=npt2-1 goto 100 endif if (apixel.lt.0.54d0) then color2(npt2)=(apixel**2/0.54d0)*100.0d0 + 200.0d0 if (color2(npt2).gt.254.0d0) color2(npt2)=254.0d0 if (color2(npt2).lt.201.0d0) color2(npt2)=201.0d0 else color2(npt2)=(bv+0.25d0)/2.0d0*100.0d0 * zblue**6 if (color2(npt2).gt.100.0d0) color2(npt2)=100.0d0 if (color2(npt2).lt.1.0d0) color2(npt2)=1.0d0 endif 100 continue c 300 continue write (1,*) npt1 c Write out image locations. do 400 i=1, npt1 write (1,350) x1(i), y1(i), A1(i), color1(i) 350 format (2x, f10.5, 5x, f10.5, 5x, f10.5, 5x, f10.5) 400 continue c if (npt2.eq.0) then npt2=2 do 410 i=1, 2 x2(i)=-2.0 y2(i)=-2.0 A2(i)= 0.0 color2(i)=0.0 410 continue endif c write (*,*) ' Writing out star data . . .' write (2,*) npt2 c do 450 i=1, npt2 write (2,350) x2(i), y2(i), A2(i), color2(i) 450 continue c 999 stop end c************************************************************************ c Amplitude interpolation program for ns data. subroutine interpA(theta, Aout, nth) implicit real*8 (a-h,o-z) parameter (ip=10000) common /int/ thint(ip), delint(ip), Aint(ip) c do 100 i=1, nth if (thint(i).lt.theta) goto 100 term=(theta-thint(i-1))/(thint(i)-thint(i-1)) Aout=Aint(i-1) + term*(Aint(i)-Aint(i-1)) return 100 continue c Aout=0.0d0 c 999 return end c************************************************************************ c Delta interpolation program for ns data. subroutine interpD(theta, delout, nth) implicit real*8 (a-h,o-z) parameter (ip=10000) common /int/ thint(ip), delint(ip), Aint(ip) c do 100 i=1, nth if (thint(i).lt.theta) goto 100 term=(theta-thint(i-1))/(thint(i)-thint(i-1)) delout=delint(i-1) + term*(delint(i)-delint(i-1)) return 100 continue c delout=-1.0d0 c 999 return end c************************************************************************