c Copyright 1999 by Robert J. Nemiroff (MTU). c Archived as http://ascl.net/nemiroff990310.html c Disclaimer: no known flaws but not guaranteed accurate. c Lensky 6.3: 98Mar3. c All sky now includes magnification bias (optional). c Latest all-sky microlensing version: computes m(tau), m(Nevents). c Computes microlensing rate on the whole sky. program lensky implicit real*8 (a-z) real*8 nms(800,12,10,23), ngi(800,12,10,23), mapp(800,12,10,23) real*8 xin(1000000,7) integer*4 i, linesdsk, linessph integer*4 iterD, iterL2, iterB2, iterMabs integer*4 iD, iL2, iB2, iMabs integer*4 iDos, iDol, iMsabs, iMlabs c common /constants/pi, kmperpc c common /magnitudes/Mabssun, Mabsdim, Mabsbrt c Open data files from galmod.for. open(1, file='/phyfac23/resrch1/nemiroff/galaxy/galmod.dsk', 2 status='old') open(2, file='/phyfac23/resrch1/nemiroff/galaxy/galmod.sph', 2 status='old') open(11, file='lensky.out', status='unknown') c Constants, pc & cgs. pi=dacos(-1.0d0) rad=pi/180.0d0 kmperpc=3.086d13 Aapp=1.341640786d0 c V band. Labssun=3.827d33 Mabssun=4.82d0 c All units in pc. Dolmin=0.0d0 Dosmin=0.0d0 c Dosmax=50000.0d0 Dosmax=15000.0d0 c dD=100.0d0 iterD=int((Dosmax - Dosmin)/dD) c Abs Magnitude limits. mapplim=17.0d0 write (*,*) ' Limiting apparant magnitude = ', mapplim Mabsbrt=-6.0d0 Mabsdim=16.5d0 dMabs=1.0d0 iterMabs=int(Mabsdim - Mabsbrt)/dMabs c Angular limits. B2min=-90.0d0 B2max=90.0d0 dB2=20.0d0 iterB2=int((B2max - B2min)/dB2) + 1 c L2min=0.0d0 L2max=360.0d0 dL2=30.0d0 iterL2=int((L2max - L2min)/dL2) write (*,*) ' Reading in galmod.dsk' c Read in disk galaxy model data. do i=1, 1000000 read (1,*,end=20) (xin(i,j), j=1, 7) c 11 format (f6.0, 1x, f5.0, 1x, f5.0, 1x, f5.1, 1x, f6.2, c 2 1x, e10.3, 1x, e10.3) enddo c 20 linesdsk=i-2 write (*,*) ' Expanding disk array.' nmsin=0.0d0 ngiin=0.0d0 c Expand input array. do i=1, linesdsk Din=xin(i,1) L2in=xin(i,2) B2in=xin(i,3) Mabsin=xin(i,4) mappin=xin(i,5) nmsin=xin(i,6) ngiin=xin(i,7) iD=int((Din - Dosmin)/dD) iL2=int((L2in - L2min)/dL2) + 1 iB2=int((B2in - B2min)/dB2) + 1 iMabs=int((Mabsin - Mabsbrt)/dMabs) + 1 mapp(iD,iL2,iB2,iMabs)=mappin c Comment next lines to isolate spheroid. nms(iD,iL2,iB2,iMabs)=nmsin ngi(iD,iL2,iB2,iMabs)=ngiin enddo c write (*,*) ' Reading in galmod.sph' c Read in spheroid galaxy model data. do i=1, 1000000 read (2,*,end=40) (xin(i,j), j=1, 7) enddo c 40 linessph=i-2 write (*,*) ' Expanding spheroid array.' c Expand input array. c do i=1, linessph Din=xin(i,1) L2in=xin(i,2) B2in=xin(i,3) Mabsin=xin(i,4) mappin=xin(i,5) nmsin=xin(i,6) ngiin=xin(i,7) iD=int((Din - Dosmin)/dD) iL2=int((L2in - L2min)/dL2) + 1 iB2=int((B2in - B2min)/dB2) + 1 iMabs=int((Mabsin - Mabsbrt)/dMabs) + 1 mapp(iD,iL2,iB2,iMabs)=mappin c Comment next lines to isolate disk. nms(iD,iL2,iB2,iMabs)=nms(iD,iL2,iB2,iMabs) + nmsin ngi(iD,iL2,iB2,iMabs)=ngi(iD,iL2,iB2,iMabs) + ngiin enddo c write (*,*) ' Starting computation.' c Reset cumulatives. Nevents=0.0d0 Nstot=0.0d0 Nsdetect=0.0d0 dDol=dD c dDos=10.0d0*dD dDos=dD c Loop through magnitude limits. do 900 mapplim=10.0d0, 20.0d0, 1.0d0 Nstot=0.0d0 Nsdetect=0.0d0 Nevents=0.0d0 c Iterate source positions. do 800 iB2=1, iterB2 B2=B2min + dble(iB2-1)*dB2 B2rad=B2*rad dOmega=sin(B2rad + pi/2.0d0)*(dB2*rad)*(dL2*rad) do 700 iL2=1, iterL2 L2=L2min + dble(iL2-1)*dL2 sigs=0.0d0 c do 600 iDos=1, iterD Dos=Dosmin + dble(iDos)*dDos c Loop source magnitudes. do 500 iMsabs=1, iterMabs Msabs=Mabsbrt + dble(iMsabs-1)*dMabs msapp=mapp(iDos,iL2,iB2,iMsabs) c Comment next line to include Magnif bias. c if (msapp.gt.mapplim) goto 600 lsapp=10.0d0**(-msapp/2.5d0) ns=nms(iDos,iL2,iB2,iMsabs) + ngi(iDol,iL2,iB2,iMsabs) c Number of sources in cell. Nsource=ns*dMabs*Dos**2*dDos*dOmega Nstot=Nstot + Nsource Nsdetect = Nsdetect + Nsource sigs=ns*dMabs*Dos**2*dDos + sigs c Loop lens distance. do 400 iDol=1, iDos-1 Dol=Dolmin + dble(iDol)*dDol Dls=Dos - Dol c zl=Dol*sin(abs(B2*rad)) c Loop lens luminosities. do 300 iMlabs=1, iterMabs Mlabs=Mabsbrt + dble(iMlabs-1)*dMabs mlapp=mapp(iDol,iL2,iB2,iMlabs) llapp=10.0d0**(-mlapp/2.5d0) ltotapp=Aapp*(llapp + lsapp) mtotapp=-2.5d0*log10(ltotapp) c UnComment next line to include Magnif bias. if (mtotapp.gt.mapplim) goto 300 Aabs=Aapp*(1.0d0 + llapp/lsapp) - llapp/lsapp Phiapp=0.5d0 Phiabs=Aabs/sqrt(Aabs**2 - 1.0d0) - 1.0d0 c Giant lenses too large --> only ms lenses. nl=nms(iDol,iL2,iB2,iMlabs) LdLsun=10.0d0**((Mabssun - Mlabs)/2.5d0) c From Mihalas & Binney p. 113. RS=(3.0d0/kmperpc)*(LdLsun**0.3125d0) bl=sqrt(4.0d0*RS*Dol*Dls*Phiabs/Dos) c Number of lenses per source. Nlens=nl*dMabs*pi*bl**2*dDol c Nevents=(Num lenses)*(Num sources). Nevents=Nevents + Nlens*Nsource 300 enddo 400 enddo 500 enddo 600 enddo write (*,*) L2, B2, sigs, dOmega c write (11,*) L2, B2, sigs 700 enddo 800 enddo c Surprise me. tau=Nevents/Nsdetect write (*,*) mapplim, tau, Nevents write (11,*) mapplim, tau, Nevents c write (*,*) mapplim, Nstot c write (11,*) mapplim, Nstot 900 enddo c 999 stop end