; olc.pro
; IDL translation of file olc.fr dated 1/9/81
;
;---this procedure computes output geometric optics light curve
;   from bending angle file. see rgf mit #3 pp.90 ++ for derivations.
;
	pro olc,olcfl,bafl,nbuf,nbuf1,ha,hb,hhat,idiag,imax,lightcurve,npts
;
	common arrays,flux,theta
;
	olcfl=''
	read, prompt="Specify output light curve file name : ",olcfl
;

	dh=1.d0/hhat
	read,prompt="Number of output points per scale height along light curve : ",that
	dj=1.d0/that
	at=dj
;
;---compute duration of light curve
;
	du=ha+exp(double(hb))+hb
	ijmax=LONG(du/at+0.5)
	du=at*ijmax

	close,4
	openw,4,olcfl
;
;---file is open and ready for writing.
;   ntot is the total # of buffers required for output
;
	ntot=1L+(ijmax-1L)/nbuf
	for n=1L, ntot do begin ; start of original do 150 loop
		print, "computing",n," of",ntot," output buffers"
		j0=(n-1L)*nbuf
		jjmin=1L
		jjmax=min([nbuf,ijmax-nbuf*(n-1L)])
	 	if(idiag eq 1)then print, "j0=",j0," jjmin=",jjmin," jjmax=",jjmax
		flux[*]=0.d0
;
;---step through entire theta file and compute contribution to this part
;   of the light curve.
;
		mmax=1L+(imax-1L)/nbuf
		if(idiag eq 1)then print, "mmax=",mmax
		for mm=1,mmax do begin ; start of DO 120 loop
			close,10
			openr,10,bafl
			lmax=min([nbuf1,imax-nbuf*(mm-1L)])
			if(idiag eq 1)then print, "lmax=",lmax
;
;---if only one point in theta buffer, it is the last point in the theta
;   file.  don't go through the routine.
;
			if(lmax le 1)then goto, LINE120
			dword=4.d0*nbuf*(mm-1)
			flskp,10,dword
			if(idiag eq 1)then print, "dword=",dword
			theta=dblarr(lmax)
			readu,10,theta
			close,10
			i0=(mm-1)*nbuf
			lmax1=lmax-1
			if(idiag eq 1)then print, "i0=",i0," lmax1=",lmax1
			for ii=1L,lmax1 do begin ; start of do 100 loop
				i=i0+ii
				zji=((i-1L)*dh+theta[ii-1L])/dj
				zji1=(i*dh+theta[ii+1L-1L])/dj
				zjmin=min([zji,zji1])
				zjmax=max([zji,zji1])
				jmin=LONG(1+zjmin-j0)
				jmax=LONG(1+zjmax-j0)
				if idiag eq 1 then begin
					print, "ii=",ii," i=",i
					print, "theta(ii-1)=",theta(ii-1)
					print, "theta(ii+1-1)=",theta(ii+1-1)
					print, "zji=",zji," zji1=",zji1
					print, "zjmin=",zjmin," zjmax=",zjmax
					print, "jmin=",jmin," jmax=",jmax
				endif
				if(jmin gt jjmax or jmax lt jjmin)then goto,LINE100
				jmin1=max([jmin,jjmin])
				jmax1=min([jmax,jjmax])
;
;---normalize to 1000 at full stellar intensity
;
				fphi=dh/(dj*(zjmax-zjmin))
				phifac=fphi*1000.d0
				if idiag eq 1 then begin
					print,"jmin1=",jmin1," jmax1=",jmax1
					print,"fphi=",fphi," phifac=",phifac
				endif
				for j1=jmin1,jmax1 do begin ; do 80 loop
					j=j1+j0
					add=phifac*(min([zjmax,j])-max([zjmin,(j-1)]))
					flux(j1-1)=flux(j1-1)+add
				if idiag eq 1 then begin
					print, "j1=",j1," j0=",j0
					print, "add=",add," flux(j1-1)=",flux(j1-1)
				endif
				endfor ; end of do 180 loop
LINE100:
			endfor ; end of do 100 loop
LINE120:
	endfor ; end of do 120 loop
;
;---this part of the light curve buffer is now complete.  write to disk.
;
		if(idiag eq 1)then print, "jjmin=",jjmin," jjmax=",jjmax," ijmax=",ijmax
		for j=jjmin,jjmax do begin ; start of do 140 loop
			flx=flux[j-1]
			if(idiag eq 1)then print, "j=",j," flx=",flx," flux(j-1)=",flux(j-1)
			writeu,4,flx
		endfor ; end of do 140 loop
		if n eq 1 then lightcurve=flux[0:jjmax-1] else lightcurve=[lightcurve,flux[0:jjmax-1]]
	endfor ; end of do 150 loop
	close,4
	npts=n_elements(lightcurve)
	return
	end