pro calhrs_rip,tnames,ih,flux,err,log
;+
;*NAME
;			calhrs_rip
;
;*PURPOSE
; Corrects GHRS data for echelle ripple
;
;*CALLING SEQUENCE:
;	calhrs_rip,tnames,ih,flux,err,log
;
; INPUTS:
;	tnames - strarr array containing reference table names
;		tnames(0) = echelle grating parameter table
;		tnames(1) = echelle "fudge" factor table
;	ih - header array (128 x nreads)
;	flux - flux array (npoints x nreads)
;
; OPTIONAL INPUT/OUTPUTS:
;	err - statistical error array. If not supplied or contains
;		a scalar, errors are not propagated.
;	log - processing log
;
; HISTORY:
;	Version 1  D. Lindler  Apr 1989
;	version 2.0 D. Lindler Sept. 9, 1991.  changed ripple function to
;		normalize it to 1.0 at car. pos = 27492(EA) and 39144(EB) and
;		epsilon = 0.0 (center of photocathode)
;	27-JAN-1992	JKF/ACC		- installed in IDL Version 2 (NOTE!)
;-
;---------------------------------------------------------------------------
	VERSION=2.0
;
; get header information
;
	s=size(flux) & ns=s(1) & nreads = n_elements(flux)/ns
	grating=ih(48)
	m=ih(50)		;spectral order
	if (grating lt 6) or (grating gt 7) then return	;is it echelle mode?
	if grating eq 6 then begin
		gmode='E-A'
		gcenter = 27492.0	;car. pos. at center of orders
	    end else begin
		gmode='E-B'
		gcenter = 39144.0
	end
	cpos=ih(43,*)
	cpos=cpos + (cpos lt 0)*65536L	;make unsigned integer
	s0=float(ih(70:71,*),0,nreads)
	deltas=float(ih(72:73,*),0,nreads)
;
; read row from first table for given grating mode
;
	table1=strtrim(tnames(0),2)
	table2=strtrim(tnames(1),2)
	tab_read,table1,tcb,table
	gmodes=tab_val(tcb,table,'GRATING')
	row=-1				;row to read
	for i=0,n_elements(gmodes)-1 do $
			if strtrim(gmodes(i)) eq gmode then row=i
	if row lt 0 then begin
		print,'CALHRS_RIP-- No row for grating '+gmode+' in table ' + $
				table1
		retall
	endif
	f=tab_val(tcb,table,'F',row)		;focal length
	beta=tab_val(tcb,table,'BETA',row)	;blaze angle
	delta=tab_val(tcb,table,'DELTA',row)	;half-angle between collimator and
					;	cross-disperser
	r0=tab_val(tcb,table,'R0',row)	;reference carrousel position
;
; update history
;
	hist=strarr(4)
	hist(0)='CALHRS_RIP version '+string(version,'(f5.2)')+ $
		': Correction for echelle ripple'
	hist(1)='    Echelle grating constant table='+table1
	hist(2)='                  Ripple constants='+table2
	hist(3)='  BETA='+string(beta,'(f7.3)')+'    DELTA='+ $
		string(delta,'(F6.3)')+'     F='+string(f,'(f9.1)')
	if !dump gt 0 then printf,!textunit,hist
	if n_elements(log) gt 0 then sxaddhist,hist,log
;
; convert to radians
;
	beta=beta*3.14159/180.
	delta=delta*3.14159/180.
;
; read rows in the second table for given order and grating mode
;
	tab_read,table2,tcb,table
	orders=tab_val(tcb,table,'SPORDER')	
	good=where(orders eq m)
	ngood=!err
	if ngood gt 0 then begin
		gmodes=tab_val(tcb,table,'GRATING',good)
		valid=bytarr(ngood)
		for i=0,ngood-1 do if strtrim(gmodes(i)) eq gmode then $
				valid(i)=1	;same grating mode?
		valid=where(valid) & ngood=!err
		if ngood gt 0 then good=good(valid)
	endif
	if ngood le 0 then begin
		print,'CALHRS_RIP-- No rows in table '+table2+' for '+ $
			gmode+' order '+strtrim(m,2)
		retall
	endif
;
; read fudge factors for this order
;
	carpos=tab_val(tcb,table,'CARPOS',good)
	a=tab_val(tcb,table,'A',good)
	b=tab_val(tcb,table,'B',good)
;
; loop on readouts
;
	clast=0
	for i=0,nreads-1 do begin
;
; get a and b for this carrousel position using linear interpolation
;
	    	if cpos(i) ne clast then begin
			linterp,carpos,a,cpos(i),aa
			linterp,carpos,b,cpos(i),bb
			clast=cpos(i)
			hist='    Carrousel position='+string(cpos(i),'(i6)')+ $
			     '  order='+string(m,'(i3)')+'    a='+ $
			     string(aa,'(F9.5)')+'   b='+string(bb,'(f9.5)')
			if !dump gt 0 then printf,!textunit,hist
			if n_elements(log) gt 0 then sxaddhist,hist,log
		endif
;
; compute sample positions relative to center of photocathode
;
		samp=(s0(i)-280.0) + findgen(ns)*deltas(i)
;
; compute epsilon
;
		eps=atan(samp/(f/0.05))
		e2=eps/2
		theta=(r0-cpos(i))*2*3.14159/65536. - beta
		x=3.14159*m*cos(theta+beta+delta)*sin(theta+e2)/ $
			sin(theta+beta+e2)
		x = aa*x+bb
		sinc2=(sin(x)/x)^2
		gnorm=cos(theta+beta+delta)/cos(theta+beta-delta+eps)*sinc2
;
; normalize by value at carrousel position = gcenter eps=0.0
;
		theta=(r0-gcenter)*2*3.14159/65536. - beta
		x=3.14159*m*cos(theta+beta+delta)*sin(theta)/ $
			sin(theta+beta)
		x = aa*x+bb
		sinc2=(sin(x)/x)^2
		gnorm_center = cos(theta+beta+delta)/ $
				cos(theta+beta-delta)*sinc2
		gnorm = gnorm/gnorm_center
		flux(0,i)=flux(*,i)/gnorm
		if n_elements(err) gt 1 then err(0,i)=err(*,i)/gnorm
		ih(66,i)=ih(66,i) or 128		;flag as done
	endfor
return
end