Viewing contents of file '../idllib/astron/contrib/varosi/code/allpro/calc_effics.pro'
;+
; NAME:
;	Calc_Effics
; PURPOSE:
;	Compute efficiency of a grating using List of flux measurements.
; CALLING:
;	Calc_Effics, eff_meas, effics
; INPUTS:
;	eff_meas = array of structures, efficiency measurements
;		collected by pro get_Eff_Data.
; KEYWORDS:
;	REF_MIRROR = array of structures with mirror reflectivity,
;			from function Read_Reflect (default = unity).
;	REF_GRATING = array of structures with grating reflectivity,
;			from function Read_Reflect (default = unity).
; OUTPUTS:
;	effics = array of structures containing absolute & groove efficiency
;		of grating and all measurements/errors used in computation.
; EXTERNAL CALLS:
;	function unique
;	function def_imhd_struct
;	pro copy_struct_inx
;	function spline
;	function stdev
; COMMON BLOCKS:
;	common Calc_Effics, efft	;structure template.
; PROCEDURE:
;	Match wavelengths, average, interpolate reflectances, compute effic.
; HISTORY:
;	Written: Frank Varosi NASA/GSFC 1994.
;-

pro Calc_Effics, eff_meas, effics, REF_MIRROR=ref_mirr, $
				REF_GRATING=ref_grat, INITIALIZE=init
   common Calc_Effics, efft
   common Calc_Effic1, tension

	if N_struct( efft ) NE 1 then begin

		error = { EFF_ERROR_v2,	WaveLength:0.0,		$
					flux_mirror:0.0,	$
					flux_grating:0.0,	$
					ref_mirror:0.0,		$
					ref_grating:0.0,	$
					absolute:0.0,		$
					groove:0.0		}

		imh = def_imhd_struct()

		efft ={	stime:"",		$
			grating:"",		$
			Lamp:"",		$
			object:"",		$
			environment:"",		$
			measurement:"",		$
			comment:"",		$
			filename:"",		$
			directory:"",		$
			flat_field:"",		$
			dev_axes: imh.dev_axes,	$
			dev_pos: imh.dev_pos,	$
			pbytes:0,		$
			psize:0,		$
			res:0,			$
			x0:0, y0:0,		$
			nx:0, ny:0,		$
			ydhms:intarr(5),	$
			counts:0L,		$
			dtime:0.0,		$
			WaveLength:0.0,		$
			flux_mirror:0.0,	$
			N_mirror:0,		$
			flux_grating:0.0,	$
			N_grating:0,		$
			ref_mirror:0.0,		$
			ref_m_file:"",		$
			ref_grating:0.0,	$
			ref_g_file:"",		$
			absolute:0.0,		$
			groove:0.0,		$
			errors:error,		$
			alpha_ang:0.0,		$
			gr_dens:0.0,		$
			gorder:1,		$
			saved:0			}
	  endif

	if keyword_set( init ) then return

	if N_struct( eff_meas ) LE 0 then begin
		message,"no measurements in List",/INFO
		return
	   endif

	wg = where( eff_meas.mirror EQ 0, Ngrat )
	wm = where( eff_meas.mirror EQ 1, Nmirr )

	if (Ngrat LE 0) OR (Nmirr LE 0) then begin
		if (Nmirr LE 0) then message,"no mirror measurements",/INFO
		if (Ngrat LE 0) then message,"no grating measurements",/INFO
		return
	   endif

	if (Ngrat EQ 1) then begin
		wgu = wg
		Nwave = 1
	 endif else $
		wgu = wg( unique( eff_meas(wg).wavelength, Nwave,/SORT,/ORIG ) )

	effics = replicate( efft, Nwave )
	copy_struct_inx, eff_meas, effics, INDEX_FROM=wgu, /RECUR_TANDEM

	if (Nwave EQ Ngrat) then begin

		effics.N_grating = 1
		effics.flux_grating = eff_meas(wgu).flux
		effics.errors.flux_grating = eff_meas(wgu).error

	 endif else begin

		for iw = 0, Nwave-1 do begin

			wavelen = effics(iw).wavelength
			w = where( eff_meas(wg).wavelength EQ wavelen, nw )
			effics(iw).N_grating = nw

			if (nw GT 1) then begin
				w = wg(w)
				effics(iw).errors.flux_grating = $
				   sqrt( stdev( eff_meas(w).flux, fmean )^2 + $
					total( eff_meas(w).error^2 )/nw )
				effics(iw).flux_grating = fmean
			 endif else if (nw EQ 1) then begin
				w = wg(w(0))
				effics(iw).flux_grating = eff_meas(w).flux
				effics(iw).errors.flux_grating=eff_meas(w).error
			  endif
		  endfor
	  endelse

	for iw = 0, Nwave-1 do begin

		wavelen = effics(iw).wavelength
		w = where( eff_meas(wm).wavelength EQ wavelen, nw )
		effics(iw).N_mirror = nw

		if (nw GT 1) then begin
			w = wm(w)
			effics(iw).errors.flux_mirror = $
				sqrt( stdev( eff_meas(w).flux, fmean )^2 + $
					total( eff_meas(w).error^2 )/nw )
			effics(iw).flux_mirror = fmean
		 endif else if (nw EQ 1) then begin
			w = wm(w(0))
			effics(iw).flux_mirror = eff_meas(w).flux
			effics(iw).errors.flux_mirror = eff_meas(w).error
		  endif
	  endfor

	if (Nwave GT 1) then effics = effics( sort( effics.wavelength ) )
	wavl = effics.wavelength
	if N_elements( tension ) NE 1 then tension = 20
	BELL = string(7b)
	errors = effics.errors

	if N_struct( ref_mirr ) GT 0 then begin
		message,"using mirror reflectance file: " + ref_mirr.file,/INFO
		effics.ref_m_file = ref_mirr.file
		wavr = ref_mirr.wavelength
		if min( wavr ) GT min( wavl ) then $
			message,"min wavelength not in range" +BELL,/INFO
		if max( wavr ) LT max( wavl ) then $
			message,"max wavelength not in range" +BELL,/INFO
		effics.ref_mirror = $
			spline( wavr, ref_mirr.reflectance, wavl, tension )
		errors.ref_mirror = $
				spline( wavr, ref_mirr.error, wavl, tension )
	 endif else begin
		message,"default mirror reflectance = 1" + BELL,/INFO
		effics.ref_mirror = 1
		errors.ref_mirror = 0
		effics.ref_m_file = ""
	  endelse

	if N_struct( ref_grat ) GT 0 then begin
		message,"using grating reflectance file: " + ref_grat.file,/INFO
		effics.ref_g_file = ref_grat.file
		wavr = ref_grat.wavelength
		if min( wavr ) GT min( wavl ) then $
			message,"min wavelength not in range" +BELL,/INFO
		if max( wavr ) LT max( wavl ) then $
			message,"max wavelength not in range" +BELL,/INFO
		effics.ref_grating = $
			spline( wavr, ref_grat.reflectance, wavl, tension )
		errors.ref_grating = $
				spline( wavr, ref_grat.error, wavl, tension )
	 endif else begin
		message,"default grating reflectance = 1" + BELL,/INFO
		effics.ref_grating = 1
		errors.ref_grating = 0
		effics.ref_g_file = ""
	  endelse

	effics.absolute = $
		( effics.flux_grating/effics.flux_mirror ) * effics.ref_mirror

	errd2 = ( errors.flux_grating / effics.flux_grating )^2 + $
		( errors.flux_mirror / effics.flux_mirror )^2 + $
		( errors.ref_mirror / effics.ref_mirror )^2

	errors.absolute = effics.absolute * sqrt( errd2 )

	effics.groove = effics.absolute/effics.ref_grating
	errors.groove = effics.groove * sqrt( errd2 + $
			( errors.ref_grating / effics.ref_grating )^2 )
	effics.errors = errors
end