C C LIGHT_IN_SLIT in ~crb/ing C C Chris Benn RGO 22 Dec 1988 C C Calculates refractive index at each wavelength for arbitrary C atmosphere, and uses this to calculate wavelength-dependent C total of light entering slit (at arbitrary angle on sky) C from a circularly-symmetric gaussian image at arbitrary C zenith angle. C C Refraction formulae from Filippenko 1982. PASP 94, 715. C C 2001 Apr 1 Unix version C Used several times since 1988 but code not checked C rigorously C C To compile, f77 lightinslit.f -o lightinslit C implicit none integer i,j,k,kk,nel1,nel2 real x1,y1,x2,y2,sum,dx,dy,pi,x,y real p,t,f,ref1,rindex(15),wave(15) real z,secz,tanz,diff(100,15) real alt,sx,sy,angle,shift,sigma real vlight,rlight,vr real suml(15) character*1 ans pi=3.1415926535 call atmosphere(t,p,f) call refract(t,p,f,wave,rindex) print *,'Write table of arcsec differential refraction ', 1 '(secz, wavelength)?' print *,' (Return for no)' read(*,'(a1)') ans if(ans.ne.' ') 1call table(rindex) do kk=1,1000000 C C Differential atm. refraction effects on light gathered at different wavelengths C C Calculate fraction of flux from Gaussian image standard C deviation SIGMA, falling within rectangular slit C SX long by SY wide, at angle ANGLE to the parallactic C angle (i.e. area covered by slit pointing to zenith), C when observing at zenith distance Z, at each of 15 C wavelengths 0.3 - 1.0 microns for which refractive index C calculated above as RINDEX(15). C Image is assumed to be circularly symmetric and of the same C shape at each frequency, and centred on the slit at C wavelength = 0.5 microns C print *,'Slit length ("), width ("), gaussian image FWHM (")?' print *,' (NB slit area will be divided into approximately' print *,' 1000 elements for integration)' read(*,*) sx,sy,sigma sigma=sigma/2.35 print *,'Seeing sigma in arcsec = ',sigma nel1=sqrt(sx/sy*1000.) nel2=sqrt(sy/sx*1000.) print *,'Dividing slit into ',nel1,' * ',nel2,' elements' print *,'Tilt of slit relative to parallactic angle (deg)?' read(*,*) angle angle=angle*pi/180. print *,'Zenith distance (deg)?' print *,' (-ve value will be interpreted as sec(ZD))' read(*,*) z if(z.lt.0) then z=acos(-1/z) else z=z*pi/180. endif print *,'Zenith distance (deg) = ',z*180./pi print *,'sec(ZD) = ',1/cos(z) write(*,108) 108 format(14x,' Coordinates of slit in sigma',/, 1 ' Wavelength Shift" ', 1 ' X1 Y1 X2 Y2 Light') do i=1,15 shift=206265.*(rindex(i)-rindex(5))*tan(z) x2=sx/2.-shift*cos(angle) x2=x2/sigma x1=x2-sx/sigma y2=sy/2.-shift*sin(angle) y2=y2/sigma y1=y2-sy/sigma call sumarea(nel1,nel2,x1,y1,x2,y2,sum) suml(i)=sum write(*,106) wave(i),shift,x1,y1,x2,y2,sum 106 format(f7.2,7x,f6.3,3x,4f7.3,f14.5) enddo C C Calculate fraction of light entering slit in v and in r C vlight=(suml(3)/2.+suml(4)+suml(5)+suml(6)/2.)/3. rlight=(suml(7)/2.+suml(8)+suml(9)/2.)/2. vr=vlight/rlight write(*,116) vlight,rlight,vr 116 format(/,' Light in v, light in r and ratio: ',3f8.4,/) enddo ! kk loop end subroutine sumarea(nel1,nel2,x1,y1,x2,y2,sum) implicit none integer nel1,nel2,i,j real x,y,x1,y1,x2,y2,dx,dy,sum,pi pi=3.1415926535 C C Divide area into NEL1 * NEL2 elements and sum fluxes from C each element C dx=(x2-x1)/float(nel1) dy=(y2-y1)/float(nel2) sum=0 do i=1,nel1 x=x1+(float(i)-0.5)*dx do j=1,nel2 y=y1+(float(j)-0.5)*dy sum=sum+1./2./pi*exp(-0.5*(x**2+y**2))*dx*dy enddo enddo end subroutine atmosphere(t,p,f) C C Obtain temperature (deg C), pressure (mm Hg) and water vapour C pressure (mm Hg) C Data for ALT = 2 km from Filippenko, formulae from Allen Ap. Q C implicit none real t,p,f,alt character*1 ans t=7 p=600 f=8 print *,'Default physical conditions:' print *,' Temperature (deg C) = ',t print *,' Atm. pressure (mm Hg) = ',p print *,' Water vapour (mm Hg) = ',f print *,'These are values for altitude = 2 km, lat = 30 deg' print *,'(La Palma Observatory has alt = 2.3 km, lat = 29 deg)' print *,'You can input new altitude or reset them explicitly' print *,'Do you want to change these values? (return for no)' read(*,'(a1)') ans if(ans.ne.' ') then alt=0 print *,'New value for temperature?' print *,' (Negative value will be interpreted as ', 1 'altitude in km)' read(*,*) t if(t.lt.0) then alt=-t t=13-3*alt endif print *,'New value for pressure?' print *,' (Negative value will be interpreted as ', 1 'altitude in km)' read(*,*) p if(p.lt.0) then alt=-p p=760.*exp(-alt/0.0293/(t+273)) endif print *,'New value for water vapour pressure?' read(*,*) f print *,' Temperature (deg C) = ',t print *,' Atm. pressure (mm Hg) = ',p print *,' Water vapour (mm Hg) = ',f endif end subroutine refract(t,p,f,wave,rindex) implicit none integer i real wave(15),rindex(15),ref1,t,p,f C C Obtain refractive indices at each wavelength (microns) C write(*,100) 100 format(' Wavelength (microns) Refractive index') do i=1,15 wave(i)=0.3+(i-1)*0.05 C C Refractive index n-1 of dry air P=760, T=15 C Edlen 1953. J Opt Soc Am 43, 339. C Coleman, Bozman and Meggers 1960. Tables of Wavenumbers NBS mono. 3 C ref1=64.328+29498.1/(146-1/wave(i)/wave(i))+255.4/ 1 (41.-1./wave(i)/wave(i)) ref1=ref1/1.e6 C C Correct for pressure P (mm of Hg) and temperature T (deg C) C Barrell 1951. J Opt Soc Am 41, 295. C ref1=ref1*p*(1+(1.049-0.0157*t)*p*1.e-6)/ 1 720.883/(1+0.003661*t) C C Correct for water vapour (Barrell 1951. J Opt Soc Am 41, 295) C F = water vapour pressure in mm of Hg C ref1=ref1-(0.0624-0.000680/wave(i)/wave(i))*f/(1+0.003661*t)*1.e-6 rindex(i)=ref1+1 write(*,101) wave(i),rindex(i) 101 format(6x,f7.2,8x,f10.8) enddo end subroutine table(rindex) implicit none integer i,j real rindex(15),secz,tanz,diff(100,15),z real ref5000 real wave(15) open(unit=7,name='light_in_slit.dat') do i=1,15 wave(i)=0.3+(i-1)*0.05 enddo write(7,101) wave 101 format('Differential atmospheric refraction in arcsec',/, 1 ' ',/, 1 'Abs" = refraction at 0.5 microns',//, 1 19x,'Wavelength (microns)',/, 1 ' ZD SecZ Abs" ',15f7.2,/) do i=1,100 secz=1.0+(i-1)*0.05 tanz=sqrt(secz*secz-1) do j=1,15 diff(i,j)=206265*(rindex(j)-rindex(5))*tanz enddo ref5000=206265*(rindex(5)-1)*tanz z=180./3.1415926535*acos(1/secz) write(7,103) z,secz,ref5000,(diff(i,j),j=1,15) 103 format(f6.2,f6.2,f8.2,2x,15f7.2) enddo print *,'Written to LIGHT_IN_SLIT.DAT' end C C Obtain coordinates C C do i=1,1000 C print *,'X1, Y1, X2, Y2?' C read(*,*) x1,y1,x2,y2 C nel1=30 C nel2=30 C call sumarea(nel1,nel2,x1,y1,x2,y2,sum) C print *,'Sum = ',sum C enddo