;
;+
;
; Computes best isothermal temperature and emission measure from
; data with more than 2 SXT filters.
;

PRO sxt_multi_teem, filters, fluxes, te, em, iso,	$
		date=date,				$
		newflux=newflux,  tecoef=tecoef

;
; INPUT PARAMETERS:
;	filters = vector of filter numbers.  Numbers are the
;		same as GT_FILTB.  Filters must be in order, 
;		thinnest first.
;	fluxes = observed fluxes, ADU per millisecond.  May be
;		2-D or 3-D; last dimension must be the same as
;		filters.  Data must be in the same order as filters.
;
; OUTPUT PARAMETERS:
;	te = log10 of temperature, K.
;	em = log 10 of emission measure, cm^-3
;	iso = "isothermality".  Vector dot product of observed
;		residuals with isothermal residuals.  Noise-free
;		isothermal values correspond to 1.
;
; OPTIONAL INPUT KEYWORDS:
;	date = date of observation in any standard Yohkoh format.
;		Used to account for the entrance filter failure on
;		November 12, 1992.  Should be input for data after
;		that date.  Default is earlier.
;
; OPTIONAL OUTPUT KEYWORDS:
;	newflux = calculated fluxes based on the isothermal T, EM.
;		Must be defined before use.
;	tecoef = vector of weights for computing temperature scaler.
;		
;
; METHOD:
;	In n-dimensional space, the isothermal fluxes for isothermal
;	gas describe a curve.  
;	A crude principal components analysis is performed.  The
;	emission measure is the scale factor.  A temperature
;	vector is found from the isothermal fluxes and used to
;	determine the observed temperatures.  The residuals are
;	compared (vector dot product) with the isothermal
;	vector residuals to determine the quality of the fit.
;
; RESTRICTIONS:
;	There is no weighting of the data.  If the data from the
;	different filters have greatly different errors, the
;	result is questionable.
;
; HISTORY:
;	Written  March 17, 1994  Barry LaBonte
;
;-

;----------------------------------------------------------------
; Handle the inputs
nfilt = N_ELEMENTS(filters)
IF( nfilt LT 3 ) THEN BEGIN
	PRINT, 'Fewer than 3 filters specified, SXT_MULTI_TEEM'
	RETURN
ENDIF
IF( KEYWORD_SET(date) EQ 0 ) THEN date = '1-JAN-92'

; Get array sizes, make trimmed dataset
sz = SIZE(fluxes)
ndim = sz(0)
dim1 = sz(1)
IF( ndim EQ 2 ) THEN BEGIN
	flux = fluxes
	dimin = dim1
  ENDIF ELSE BEGIN
	dim2 = sz(2)
	flux = REFORM(fluxes, dim1*dim2, nfilt)
	dimin = dim1 * dim2
ENDELSE

; Work on nonzero data only
mf = flux(*,0)
FOR i=1,nfilt-1 DO BEGIN
	mf = mf < flux(*,i)
ENDFOR
nonz = WHERE(mf GT 0.)
flux = flux(nonz, *)
dimf = N_ELEMENTS(flux(*,0))


;-------------------------------------------------------------------
; Define reference temperatures, fluxes
; Temperature step size matches the database
delt = 0.05
nref = FIX(2.5 / delt) + 1
xref = FINDGEN(nref)
tref = 5.5 + delt * FINDGEN(nref)
fref = FLTARR(nref, nfilt)
FOR i=0,nfilt-1 DO BEGIN
        fref(*,i) = SXT_FLUX(tref, filters(i), DATE=date)
ENDFOR
lref = ALOG10(fref)
; Emission measure scaler
emscalref = TOTAL(lref,2)/nfilt
; Flux residual vector
residref = lref - REBIN(emscalref, nref, nfilt)

; Solve for temperature scaler for this combination of filters
coeft = FLTARR(nfilt)
FOR i=0,nfilt-1 DO BEGIN
	result = POLY_FIT(tref, residref(*,i), 1)
	coeft(i) = result(0,1)
ENDFOR
coefn = coeft/TOTAL(ABS(coeft))
tescalref = FLTARR(nref)
FOR i=0,nfilt-1 DO BEGIN
	tescalref = tescalref + residref(*,i) * coefn(i)
ENDFOR
; Determine limiting range of temperature
smax = MAX(tescalref, imax)
IF(imax NE nref-1) THEN PRINT, $
	 'Isothermal Temperature undefined for this filter set above log(T) = ', tref(imax)
; Get the correct coefficients
coefc = FLTARR(nfilt)
FOR i=0,nfilt-1 DO BEGIN
        result = POLY_FIT(tescalref, residref(*,i), 1)
        coefc(i) = result(0,1)
ENDFOR

; Isothermal residual vector and its squared norm
isoref = residref - (tescalref # coefc)
normref = TOTAL( isoref^2, 2 )

;----------------------------------------------------------------------
; Work in log space, remove emission measure scaler, compute
; temperature scaler
lflux = ALOG10(flux)
emscalobs = TOTAL( lflux, 2 )/nfilt
residobs = lflux - REBIN(emscalobs, dimf, nfilt)
tescalobs = FLTARR(dimf)
FOR i=0,nfilt-1 DO BEGIN
	tescalobs = tescalobs + residobs(*,i) * coefn(i)
ENDFOR
; Find the valid data
onscale = WHERE(tescalobs LE smax)
obsfinal = residobs - (tescalobs # coefc)
obsnorm = SQRT( TOTAL( obsfinal^2, 2 ) )

temps = FLTARR(dimf)
ems = FLTARR(dimf)
; Do the critical interpolations all at once
temps(onscale) = DSPLINE( tescalref, tref, tescalobs(onscale) )
ems(onscale) = emscalobs - DSPLINE( tref, emscalref, temps(onscale) )
fls = FLTARR(dimf,nfilt)
FOR i=0,nfilt-1 DO BEGIN
	fls(onscale,i) = DSPLINE( tref, isoref(*,i), temps(onscale) )
ENDFOR
norms = DSPLINE( tref, normref, temps(onscale) )

; Now measure the "isothermality"
isos = FLTARR(dimf)
isos(onscale) = TOTAL( obsfinal(onscale,*) * fls(onscale,*), 2 )/norms(onscale)

;-----------------------------------------------------------------------
;Define the output arrays and fill them
te = FLTARR(dimin)
em = FLTARR(dimin)
iso = FLTARR(dimin)
te(nonz) = temps
; Note - SXT_FLUX gives ADU/1000 msec for logEM=44; we use ADU/1msec
em(nonz) = ems + 47.
iso(nonz) = isos
IF( ndim GT 2 ) THEN BEGIN
	te = REFORM(te, dim1, dim2)
	em = REFORM(em, dim1, dim2)
	iso = REFORM(iso, dim1, dim2)
ENDIF
IF( KEYWORD_SET(newflux) NE 0 ) THEN BEGIN
	newflux = FLTARR(dimin, nfilt)
	newf = FLTARR(dimf)
	FOR i=0,nfilt-1 DO BEGIN
		newf(onscale) =  DSPLINE( tref, fref(*,i), temps(onscale) )
		newflux(nonz, i) = newf * 10^ems
	ENDFOR
	IF( ndim GT 2 ) THEN newflux = REFORM(newflux, dim1, dim2, nfilt)
ENDIF
tecoef = coefc

RETURN
END