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









