Viewing contents of file '../idllib/ghrs/pro/calhrs_wav.pro'
pro calhrs_wav,tflag,tnames,dc_offtab,ih,wave,log
;+
;			calhrs_wav
;
; Routine to compute wavelength array from dispersion coefficients.
;
; CALLING SEQUENCE:
;	calhrs_wav,tflag,tnames,ih,wave,log
;
; INPUTS:
;	tflag - integer flag (if greater than 0 thremal motion correction
;		is done)
;	tnames - reference table names (string array)
;		tnames(0) - carrousel calibration table
;		tnames(1) - dispersion coef. table
;
;	dc_offtab - name of dispersion coef. table from which an
;		average sample offset is to be taken.  It should only
;		have a single row and should be for the same grating
;		mode.  The only column used is AVE_OFFSET.  If blank
;		the no offset is added
;
; INPUT/OUTPUT:
;	ih - header arrays (128xN) for the flux.
;		upon output the order number
;		placed into the array
; OUTPUTS:
;	wave - wavelength array npoints by N
;
; OPTIONAL INPUT/OUTPUT:
;	log - processing log
;
; METHOD:
;	The dispersion coefficients are read from table tnames(1).
;	Rows are selected by grating and carrousel position.  in
;	the event that the exact carrousel position of the data
;	is not found, linear interpolation between carrousel positions
;	is used.
;
;	The wavelengths are computed from sample positions using the
;	relation:
;		
;		s = a0 + a1*m*w + a2*m*m*w*w + a2*m + a4*w +
;			a5*m*m*w + a6*m*w*w + a7*m*m*m*w*w*w +
;			a8*m*m + a9*w*w
;
;	where s is the photocathode sample position computed from
;		the starting sample and delta sample stored in the
;		header array
;	      m is the spectral order computed using the carrousel
;		calibration parameters in tname(0) for the echelle
;		modes. It is set to 1 for the first order gratings.
;	      a0,...,a9 are the dispersion coef.  a7 through a9 are
;		experimental and are not used by the STSCI pipeline.
;	      w is the wavelength.
;	
;	The wavelengths are computed by solving this equation using
;	the Newtons iterative method.
;
;	The order numbers are computed using routine CALHRS_ORD.
;
; HISTORY:
;	version 1   D. Lindler   March 89
;	version 1.1 D. Lindler  Jan 1991  Uses a0_offsets to offset
;		wavelengths for FP-SPLIT observations instead of
;		reinterpolation in DCTAB.  Added DC_OFFTAB.
;		added thermal motion correction
;	version 2.0 D. Lindler May 1992  Allows dispersion coefficient
;		table that contains global fit to the dispersion coef.
;-
;---------------------------------------------------------------------------

	VERSION = 2.0

;
; extract some needed info from the header
;
	nreads = n_elements(ih)/128
	szero = float(ih(70:71,*),0,nreads)
	if min(szero) eq 0.0 then begin
		print,'CALHRS_WAV--called without first performing'+ $
				' mapping function'
		retall
	endif
	deltas = float(ih(72:73,*),0,nreads)
	steps = fix(1.0/deltas + 0.5)		;number of substeps
	ydefs=ih(42,*)				;ydeflections
	sclamp=ih(39)*2+ih(38)			;spectral cal lamp
	grating=ih(48)				;grating mode
	if (grating lt 1) or (grating gt 7) then begin
		print,'CALHRS_WAV--Input data is not a spectral mode'
		retall
	endif
	cpos=ih(43,*)
	cpos=(cpos lt 0)*65536L + cpos		;unsigned integer
	gmode=['NDF','G-1','G-2','G-3','G-4','G-5','E-A','E-B']
	gmode=gmode(grating)
	binid=ih(53)
	if binid eq 2 then aperture='LSA' else aperture='SSA'
	if sclamp eq 1 then aperture='SC1'
	if sclamp eq 2 then aperture='SC2'
;
; compute spectral orders
;
	calhrs_ord,tnames(0),gmode,aperture,cpos,ydefs,orders
;
; update log
;
	hist = strarr(3)
	hist(0)='CALHRS_WAV version '+string(version,'(f5.2)')+ $
		': Wavelength computation for '+gmode
	hist(1)='    Carrousel calibration table '+strtrim(tnames(0),2)
	hist(2)='    Dispersion coefficient table '+strtrim(tnames(1),2)
;
; get dispersion coefficient offset if dc_offtab supplied
;
	if strtrim(dc_offtab) ne '' then begin
		table_ext,dc_offtab,'AVE_OFFSET',ave_offset
		ave_offset = ave_offset(0)
		hist =[hist,strarr(2)]
		nhist = n_elements(hist)
		hist(nhist-2) = '    Dispersion Coefficient offset table '+dc_offtab
		hist(nhist-1) = '        Dispersion coefficient offset ='+ $
			    strtrim(ave_offset,2)+' sample units'
	    end else begin
		ave_offset = 0.0
	end
;
; create output wavelength array
;
	ns = 500*steps(0)	;number of points in the output array
	wave=dblarr(ns,nreads)
;
; find table rows in dispersion coef. table for the given grating
;
	dctab = strtrim(tnames(1),0)
	tab_read,dctab,tcb,table
	gmodes=tab_val(tcb,table,'GRATING')
	nrows=n_elements(gmodes)		;number of rows
	match=intarr(nrows)			;flags for matches
	for i=0,nrows-1 do if strtrim(gmodes(i)) eq gmode then match(i)=1
	valid=where(match)			;valid rows
	nvalid=!err				;number of them
	if nvalid lt 1 then begin
		PRINT,'CALHRS--No valid rows in '+dctab+' for grating '+gmode
		retall
	endif
;
; determine if table contains global coefficients for the grating or
; it contains dispersion coefficients.
;
	tab_col,tcb,'F00'
	if !err eq -1 then begin	;It contains disp. coef.
;
; read the dispersion coefficients for the valid rows ----------------------
;
	    cpos_dc=tab_val(tcb,table,'CARPOS',valid)	;carrousel positions
	    dc=dblarr(nvalid,10)
	    for i=0,9 do begin
		col_name='A'+strtrim(i,2)
		if i gt 6 then tab_col,tcb,col_name else !err=1 ;check for
;						presence of optional columns
		if !err gt 0 then dc(0,i)=tab_val(tcb,table,col_name,valid)
	    endfor
	  end else begin
;
; read global coefficients for the grating ----------------------------------
;
	    valid = valid(n_elements(valid)-1)	;row to read
	    colnames = ['CAP_A','CAP_C','MCENTER','F00','F01','F02', $
			'F10','F11','F12','F20','F21','F22','F30', $
			'F31','F32','F40','F41','F42','F50','F51','F52']
	    global_coef = dblarr(21)
	    for i=0,20 do global_coef(i) = tab_val(Tcb,table,colnames(i),valid)
;
; Correct for thermal/time motion
;
	    if tflag then begin
;
;  Thermal motion for THERMISTOR1 and THERMIStOR2
;
		for itemp = 1,2 do begin
		    st = strtrim(itemp,2)
		    tab_col,tcb,'THERMISTOR'+st
		    if (!err ge 0) then begin
		        thermistor = tab_val(tcb,table,'THERMISTOR'+st,valid)
		        tc = tab_val(tcb,table,'TC'+st,valid)
			if strtrim(thermistor,2) ne '' and (tc ne 0.0) then $
			begin
			    temp_dc = tab_val(tcb,table,thermistor,valid)
		            temp_obs = sxpar(log,thermistor)
		    	    if !err lt 0 then begin
			    	print,'CALHRS_WAV- Keyword '+thermistor+ $
				      ' not found in header: thermal' + $
				      ' correction not done' 
			      end else begin
			        cor = tc*(temp_obs-temp_dc)
			        global_coef(3) = global_coef(3) + cor
			        hist = [hist,'Thermal correction for '+ $
					thermistor+' = '+strtrim(cor,2)+ $
					' sample units']
			    end
			end
		    end
		end ; for itemp
;
; time motion
;
		TAB_COL,tcb,'JD0'
		if !err gt 0 then begin
		    dsdt = tab_val(tcb,table,'DSDT',valid)
		    if dsdt ne 0.0 then begin
		    	jd0 = tab_val(tcb,table,'JD0',valid)
		    	jd1 = tab_val(tcb,table,'JD1',valid)
		    	obstime = string(byte(ih(54:65),0,24))
		    	obsjd = jul_date(obstime)
		    	obsjd = obsjd<jd1>jd0		;do not extrapolate
		    	cor = (obsjd - jd0)*dsdt
		    	hist = [hist,'Time Motion for '+obstime+' = '+$
					strtrim(cor,2)+' sample units']
		    	global_coef(3) = global_coef(3) + cor
		    end
		end
;
; dispersion change with temperature
;
		tab_col,tcb,'dthermistor'
		if !err gt 0 then begin
		    thermistor = tab_val(tcb,table,'DTHERMISTOR',valid)
		    dtc = tab_val(tcb,table,'DTC',valid)
		    if strtrim(thermistor,2) ne '' and (dtc ne 0.0) then begin
			temp_dc = tab_val(tcb,table,thermistor,valid)
		        temp_obs = sxpar(log,thermistor)
		    	if !err lt 0 then begin
			    	print,'CALHRS_WAV- Keyword '+thermistor+ $
				      ' not found in header: thermal' + $
				      ' correction not done' 
			      end else begin
			        cor = dtc*(temp_obs-temp_dc)
			        global_coef(6) = global_coef(6) + cor
			        hist = [hist,'Dispersion correction for '+ $
					thermistor+' = '+strtrim(cor,2)+ $
					' sample units/m*lambda']
			end
		    end
		end

	    end ; if tflag
	end
;
; dump hist
;
	if n_elements(log) gt 0 then sxaddhist,hist,log
	if !dump gt 0 then printf,!textunit,hist

;
; compute wavelengths
;
	cpos_last=0				;last carrousel position done
	m_last=0				;last order number
;
; loop on readouts
;
	for iread=0,nreads-1 do begin
	    cpos1=cpos(iread)			;carrousel position for this
						;  readout
	    m=orders(iread)			;order number
;
; create new starting guess if last readout's wavelengths are not
; appropriate
;
	    if (m ne m_last) or (cpos1 ne cpos_last) then begin
		     acoef = [0.0, 3301.815, 4020.502, 4615.660, 5539.488, $
			      33013.449, 62886.313, 63194.086]
		     ccoef = [0.0, 11106.51, 54807.90, 30600.90, 14887.26, $
			      49183.15, 39021.18, 50575.02]
		     wguess = acoef(grating)/m * $
				sin((ccoef(grating)-cpos1)/10430.378)
		     w=replicate(wguess,ns) 
	    end
	    m_last=m
;
; Get disp. coef. if it is a new carrousel position
;
	    a0_offset = dblarr(10)
	    a0_offset(0) = ave_offset

	    if cpos1 ne cpos_last then begin	;need new disp. coef?
		diff = cpos1 - cpos_last	;change in carrousel position
		if (diff lt 4) and (diff gt 0) then begin	;FP-SPLIT
			case gmode of
				'G-1' : magnification = 20.8
				'G-2' : magnification = 20.1
				'G-3' : magnification = 20.1
				'G-4' : magnification = 20.0
				'G-5' : magnification = 22.0
				'E-A' : magnification = 18.7
				'E-B' : magnification = 18.9
			endcase
			a0_offsets = magnification * [0.0, 0.248, 0.483, 0.752]
			a0_offset(0) = ave_offset + a0_offsets(diff)
			hist='   FPSPLIT offset for carrousel position '+ $
				string(cpos1,'(I6)') + ' = ' + $
				string(a0_offsets(diff),'(F6.2)')+ $
				' sample units'
			if n_elements(log) gt 0 then sxaddhist,hist,log
			if !dump gt 0 then printf,!textunit,hist
		    end else begin
			cpos_last=cpos1
;
; compute new dispersion coefficients
;
			calhrs_dc,global_coef,cpos_dc,dc,cpos1,a,log
		end
	    endif
;
; Compute sample position vector
;
	    s=findgen(ns)*deltas(iread)+szero(iread)
;
; compute wavelengths
;
	    calhrs_stow,s,m,a+a0_offset,w
	    if !err lt 0 then begin
		hist=' ***Wavelength computation failed to converge for '+ $
			'readout '+string(iread,'(i5)')
		if n_elements(log) gt 0 then sxaddhist,hist,log
		if !dump gt 0 then printf,!textunit,hist else print,hist
	    endif
	    wave(0,iread)=w
	    ih(50,iread)=m
	endfor
return
end