; golc.pro - This is a virtually line-by-line translation of golc.fr
;		which is an ancient fortran program that computes
;		geometric optics light curves from refractivity profiles.
; 	     The original version uses buffers to disk.
; Revisions:
;	2001-11-04 : rgf - original IDL version
;
;
;---This program creates geometric optics occultation light
;   curves from a variety of inputs.  see rgf mit #3 pp.90 ++ for
;   derivations.
;   for 2-dimensional waves, see rgf mit#3 pp.128ff.
;
	COMMON ARRAYS,ar1,ar2
	NBUF=50000L
	NBUF1=NBUF+1
	ar1=fltarr(NBUF)
	ar2=fltarr(NBUF1)
	print,'Program golc run at '+systime()

	rfl=''
	bafl=''
	pfl=''
	olcfl=''
;
;---specify input data option
;
LINE30:
	read,prompt='0-Stop 1-Refractivity 2-Bending angle 3-Phase front : ',igo
	if(igo eq 0) then stop
	read,prompt="Print diagnostics? 0-No 1-Yes : ",idiag
	if igo eq 1 then GOTO,LINE40
	if igo eq 2 then GOTO,LINE60
	if igo eq 3 then GOTO,LINE80
;
;---input is refractivity
;
LINE40:
	read,prompt='0-Two-dimensional waves 1-Use existing file 2-Create refractivity file : ',irgo
	if(irgo eq 1)then prompt="Specify existing refractivity file name : "
	if(irgo eq 2)then prompt="Specify new refractivity file name : "
	if(irgo ne 0)then read,prompt=prompt,rfl
	print, "Specify refractivity profile parameters : "
	read, prompt="File begins how many scale heights above half-light ? ",ha
	read, prompt="File ends how many scale heights below half-light ? ",hb
	read, prompt="Number of atmosphere shells per scale height : ",hhat
	read, prompt="Number of vertical wavelengths per scale height : ",wlpsh
	read, prompt="Phase (degrees) of vertical wave at half-light : ",psideg
	dtr=!dpi/180.d0
	psi=psideg*dtr
	irmod=0 ; by default
;
;---specify 2-d wave parameters if necessary
;
	if(irgo ne 0)then GOTO, LINE50
	print, "Specify refractivity model :"
	print, "  0 - Sum of independent horizontal and vertical waves "
	read,prompt="  1 - Constant amplitude gravity waves : ",irmod
	read,prompt="Ratio of planetary radius to scale height : ",rhat
	read,prompt="Number of scale heights per horizontal integration bin : ",dxhat
	read,prompt="Number of scale heights per horizontal wavelength : ",shphw
	if(irmod ne 0)then GOTO,LINE50
	read,prompt="Phase (degrees) of horizontal wave at mid-point : ",betdeg
	beta=betdeg*dtr
;
;---construct refractivity profile if necessary
;
LINE50:
	refrfl,irgo,rfl,ha,hb,hhat,wlpsh,psi,imax
	read,prompt="Specify new bending angle file name : ",bafl
	read,prompt="Fractional amplitude of vertical wave : ",eps
	if(irmod eq 0) then $
		read,prompt="Fractional amplitude of horizontal wave : ",del
;
;---create bending angle file
;
	bendfl,rfl,bafl,nbuf,ha,hb,hhat,eps,irgo,irmod,wlpsh, $
     		psi,rhat,dxhat,shphw,beta,del,idiag,imax
	GOTO,LINE100
;
;---use existing bending angle file
;
LINE60:
	read,prompt="Specify existing bending angle file name : ",bafl
	read,prompt="Specify number of points per scale height : ",hhat
	read,prompt="File begins how many scale heights above half-light ? ",ha
	read,prompt="File ends how many scale heights below half-light ? ",hb
	imax=1.5+hhat*(ha+hb)
	GOTO,LINE100
;
;---convert from thin screen phase to bending angle file
;
LINE80:
	read,prompt="Specify existing thin screen parent file name : ",pfl
;***eventually, finish this part to create output bending angle
;   file, specifying eps here
	print,'NOT YET IMPLEMENTED!'
	STOP
LINE100:
;
;---at this stage, all input options have converged, and a bending
;   angle file exists which will be transformed into an output light curve
;   in standard data-file format.
;
	olc,olcfl,bafl,nbuf,nbuf1,ha,hb,hhat,idiag,imax,lightcurve,npts
	print,'Light curve is contained in array lightcurve[], with npts=',npts
	;GOTO,LINE30
	end