; bendfl.pro
; based on BENDFL.FR
;
;---this procedure computes bending angle information from
;   perturbation refractivity files.  the bending angle is normalized
;   by the half-light bending angle for an isothermal atmosphere.
;   see rgf mit #3 pp.91++ and p.128++.
;

; Revisions:
;	2001 Nov 4 - rgf - IDL version with no buffers

	pro bendfl,rfl,bafl,nbuf,ha,hb,hhat,eps,irgo,irmod,wlpsh, $
     		psi,rhat,dxhat,shphw,beta,del,idiag,imax

	common arrays,znu,theta
	pi=!dpi
;
	if(irgo eq 0)then GOTO,LINE110
	dh=1.d0/hhat
	fac=2.d0*sqrt(hhat/pi)*exp(-ha-dh)
;
;---zero the output file
;
	close,1
	openw,1,bafl
	th=dblarr(imax)
	writeu,1,th
	close,1
;
;---open the refractivity file
;
	close,10
	openr,10,rfl
;
;---compute number of buffers required for i/o
;
	ll=1L+(imax-2)/nbuf
	mm=1L+(imax-1)/nbuf
	for l=1,ll do begin ; beginning of original DO 100 loop
	if(idiag eq 1)then print, "reading",l," of",ll," refractivity buffers"
;
;---read the l'th refractivity buffer, and
;   convert from perturbation refractivity to full normalized
;   refractivity
;
		lmax=LONG(min([nbuf,imax-(l-1)*nbuf]))
		k0=(l-1L)*nbuf
		if(idiag eq 1)then print, "lmax=",lmax," k0=",k0
		znuval=0.d0
		for n=1L,lmax do begin ; original do 25 loop
			fk=k0+n
			readu,10,znuval
			znu[n-1]=fac*exp(fk*dh)*(1.d0+eps*znuval)
			if(idiag eq 1)then print, "n=",n," znu(n-1)=",znu[n-1]
		endfor ; end of do 25 loop
;
;---update theta array, a buffer at a time
;
		if(idiag eq 1)then print, "l=",l," mm=",mm
		for m=1L,mm do begin ; beginning of original DO 60 loop
;
;---read the current theta buffer
;
			openr,1,bafl
			dword=4.d0*(m-1)*nbuf
			flskp,1,dword
			mmax=min([nbuf,imax-(m-1)*nbuf])
			if(idiag eq 1) then print, "m=",m," mmax=",mmax
			theta=dblarr(mmax)
			readu,1,theta
			close,1
			i=(m-1)*nbuf
			openu,1,bafl
			flskp,1,dword
			for m1=1L,mmax do begin ; start of original DO 50 loop
				i=i+1
				sum=0.d0
				k1max=min([nbuf,i-1-k0])
				if(idiag eq 1)then print, "m1=",m1," i=",i," k1max=",k1max
				if(k1max lt 1)then goto,LINE45
				k=k0+1L
				fimk=double(i-k)
				fac1=sqrt(max([0.d0,fimk-1.d0]))
				fac2=sqrt(fimk)
				fac3=sqrt(max([0.d0,fimk-2.d0]))
				for k1=1L,k1max do begin ; start of do 40 loop
					add=znu(k1-1)*(-2.d0*fac1+fac2+fac3)
					sum=sum+add
					if(idiag eq 1) then begin
						print, "k1=",k1," k=",k," fimk=",fimk
						print, "fac1=",fac1," fac2=",fac2
						print, " fac3=",fac3," add=",add
						print, "sum=",sum," znu(k1)=",znu(k1-1)
					endif
;
;---update counters and factors for next pass
;
					fimk=fimk-1.d0
					fac2=fac1
					fac1=fac3
					fac3=sqrt(max([0.d0,fimk-2.d0]))
				endfor ; end of DO 40 loop
				theta(m1-1)=theta(m1-1)+sum
;
;---write the updated theta value, even if no changes have been made.
;
LINE45:
				if(idiag eq 1)then print, "m1=",m1," theta(m1-1)=",theta(m1-1)
				writeu,1,theta(m1-1)
			endfor ; end of DO 50 loop
			close,1
		endfor ; end of DO 60 loop
	endfor ; end of original do 100 loop
	close,10
	if idiag eq 1 then begin
		print,'at end of computing bending angle array...'
		STOP
	endif
	return
;
;---this section used for two-dimensional refractivity profiles
;
LINE110:
	close,1
	openw,1,bafl
	fac=-dxhat/sqrt(2.d0*pi*rhat)
	if(idiag eq 1)then begin
		print, "dxhat=",dxhat," hhat=",hhat
		print, "eps=",eps," wlpsh=",wlpsh
		print, "psi=",psi," rhat=",rhat
		print, "shphw=",shphw," beta=",beta
		print, "del=",del," fac=",fac
	endif
;120
	for i=1L,imax do begin ; start of do 180 loop
		zni=max([10.d0,sqrt(2.d0*rhat*double(i)/hhat)/dxhat])
		ni=zni
		ni1=ni+1
		kmax=ni+ni1
		sum=0.d0
		if(idiag eq 1)then  print, "i=",i," ni=",ni
; 130 was here
		for k=1L,kmax do begin ; start of do 160 loop
			j=k-ni1
			if(irmod eq 0)then sum=sum+ $
				twodw(ha,i,j,dxhat,hhat,eps,wlpsh,psi, $
					rhat,shphw,beta,del,idiag)
			if(irmod eq 1)then sum=sum+ $
				twod1(ha,i,j,dxhat, $
					hhat,eps,wlpsh,psi,rhat,shphw,idiag)
		endfor
		thet=sum*fac
		writeu,1,thet

		format170="(' completed',i5,' of ',i5,' points')"
		if(i mod 100 eq 0)then print,format=format170,i,imax
		if(idiag eq 1) then print, "thet = ",thet
	endfor ; end of do 180 loop

	close,1
	return
LINE200:
	print, "end of file reached in do 25 loop of routine bendfl"
	print, "ll=",ll," l=",l," lmax=",lmax
	print, "k=",k," k0=",k0," n=",n
	stop
LINE300:
	print, "end of file reached in d0 35 loop of routine bendfl"
	print, "ll=",ll," l=",l," mm=",mm
	print, "m=",m," mmax=",mmax," m1=",m1
	stop
	end