; Copyright (c) 1998-2009 A.P. Hitchcock  All rights reserved
;+
; NAME:
;		  LINE_FIT
;
; LAST CHANGED: --------------------------------  30-Jul-09
;
; This procedure selects inputs, fits spectrum at each point along a line spectral scan
; to sum of models.  Uses COMMON blocks to store and transfer data. Maximum 5 components !
;
; KEYWORDS
;	AXIS - if set, indicates called from AXIS
;	COEFF - name of file listing sources of reference spectra
;	VERBOSE - if set, additional information usppplied
;
; COMMON
; 	AXIS_COM - AXIS
;	stack_process_com
; 	BSIF_COM - nc data
;
;MODIFICATION HISTORY:
; (23-Dec-99 aph) first developed from stack_fit
; (27-feb-00 aph) add groupID to get_text
; (30-nov-00 cgz) Replaced call of nwin with explicit window statement
; (04-Jun-03 aph) adapt to use ax_par_save & ax_par_load
; (22-sep-06 aph) fix output to display constant & chi-square properly
;                 remove groom (left over from kneedler); use axis_log;
; (30-jul-09 aph) changed analcom to stack_process_com (update to stack_process)
;-

PRO line_fit, axis=axis, coeff = coeff, verbose = verbose
@axis_com
@stack_process_com
@bsif_com

on_error,2

; ----------- read-in linescan data
file = pickfile2(/READ, FILTER='*.axb', /LPATH, DEFPATH=defpath)
if strlen(file) GT 0 THEN BEGIN  ; bail-out if no filename
	tmp_ls = axb_load(file=file)		;  read-in linescan
	t=ax_name(file)
	fileroot = t(0) + t(1)
	if keyword_set(axis) then $
		fileroot = get_text(prompt = 'Root name for output files',val=fileroot, group = axis_ID) $
		 else fileroot = get_text(prompt = 'Root name for output files',val=fileroot)
	fileroot = fileroot + '_'
	ev = tmp_ls.x   ; set Energy-axis of linescan data


; -------------- define coefficents at energies of images -------------
IF NOT keyword_set(coeff) THEN BEGIN
	t = ' '
    t = dialog_message('Read fit parameter file ?', /question)
	if t(0) EQ 'Yes' then begin
	    par_file = pickfile2(title = ' Select fit parameter file', $
	           filter='*.par', /LPATH, DEFPATH=defpath )
	    if strlen(par_file(0)) NE 0 then coeff = par_file
	endif
ENDIF
IF keyword_set(coeff) THEN BEGIN
; ------- read spectral file info from a file
	check = findfile(coeff)
	if strlen(check(0)) NE 0 then pars = ax_par_load(coeff) else goto, get_user_info
	if n_tags(pars) EQ 0 then begin
		goto, get_user_info
	endif else begin
		ncomp = pars.n
		comp_names = pars.names
		comp_files = pars.files
	endelse

; ------ if fails then ask user to locate reference files one-by-one
ENDIF ELSE BEGIN
 get_user_info:
	if keyword_set(axis) then ncomp = get_num(prompt = 'Number of components',val=ncomp, group=axis_ID) $
	   else ncomp = get_num(prompt = 'Number of components',val=ncomp)
	if ncomp LE 0 then return
	comp_names = strarr(ncomp)
	comp_files = strarr(ncomp)
; ------- read coefficent information from spectra files
	for i = 0, ncomp-1 do begin
	 	text = 'Spectral file for component ' + strcompress(string(fix(i+1)))
		axis_log, text
		comp_files(i) = pickfile2(title = ' Spectrum of component '+ strtrim(string(i),2), $
	           filter='*.txt', /LPATH, DEFPATH=defpath )
	    if comp_files(i) EQ '' then return
		tmp = spc_load(file = comp_files(i))
;; ------ define name - keep short as used in map name
		text = 'short name of component ' + strcompress(string(fix(i+1)))
		if keyword_set(axis) then comp_names(i) = get_text(prompt = text, val = tmp.dl, group = axis_ID) $
		    else comp_names(i) = get_text(prompt = text, val = tmp.dl)
		comp_names(i) = strtrim(string(comp_names(i)),2)
	endfor

; ------ save coefficient array for later use (optional)
	par_file = pickfile2(filter='*.par', /write, path = DefPath, $
	           title = 'Name of fit parameter file')
	if strlen(par_file(0)) NE 0 then ax_par_save,par_file,ncomp,comp_names,comp_files
ENDELSE

; read in files and interpolate to the energy scale of the image set
	model = fltarr(ncomp, n_elements(ev))
	CurBuf = 0
	for i = 0, ncomp-1 do begin
		tmp = spc_load(file=comp_files(i))
; ----- interpolate to the same energies as the image sequence
		model(i,*)=interpol(tmp.d,tmp.x,ev)
		tmp = {t:'1d', d: model(i,*), x: ev, xl: 'X', dn:model(i,*), dl: 'interp '+ tmp.dl}
		CurBuf = i+1
		if  n_tags(tmp) NE 0 AND keyword_set(AXIS) THEN BEGIN
			HANDLE_VALUE, Data(Curbuf), tmp, /SET
			PlotBuf, Curbuf
		endif
	endfor


; -------START OF Linescan fitting procedure ---------
;	print, 'Fitting spectrum at each pixel to ', strcompress(string(ncomp)),' components'
	n_E = n_elements(ev)

; groom reference spectra to have zero background (but they are
; assumed be normalized to each other in some way)
;	if  keyword_set(AXIS) THEN BEGIN
;		nzero = get_num(Prompt='#bkground points for model spectra',Val = 0, group = axis_id)
;		nzero_img = get_num(Prompt='#bkground points for linescan spectra?', Val=0, group = axis_id)
;	endif else begin
;		nzero = get_num(Prompt='#bkground points for model spectra',Val = 0)
;		nzero_img = get_num(Prompt='#bkground points for linescan spectra?', Val=0)
;	endelse
;	Widget_control, /hourglass
;	if nzero GT 0 then begin
;		for i = 0, ncomp-1 do begin
;			tmp = reform(model(i,*))
;			groom1, tmp, bak, inte, nzero
;			model(i,*) = tmp
;			print,'reference spectrum ', strcompress(string(i+1)),' : subtract', bak,'; integrated intensity = ',inte
;		endfor
;	ENDIF

; switch plotting depending on keyword AXIS
	if NOT keyword_set(axis) then begin
	; plot reference signals
;		nwin,0
		window, 0, xsize=300, ysize=300
		wset,0
		plot, ev,model(0,*)
		for i = 1, ncomp-1 do oplot,ev,model(i,*)
;		oplot, ev, 10*(line/max(line))
	endif

; ------------ set up (x,y) scales from bsif_com information
	n_cols = n_elements(tmp_ls.x)
	n_rows = n_elements(tmp_ls.y)
	x_step = tmp_ls.x(1) - tmp_ls.x(0)	; assumes linear E-scal which is case when
	y_step = tmp_ls.y(1) - tmp_ls.y(0)  ; linescan converted to OD and written as AXB
	img_x = tmp_ls.x
	img_y = tmp_ls.y
; --------- set-up output arrays
	const=fltarr(n_rows)
	comp = fltarr(ncomp, n_rows)
	chisq=const
	bk=const
	int=const
	linear=const
	residual = float(tmp_ls.d)
	if NOT keyword_set(axis) then begin
		window, 1, xpos=0, ypos=200, xsize=n_cols, ysize=n_rows
		window, 2, xpos=0, ypos=400, xsize=n_cols, ysize=n_rows
		if ncomp EQ 3 then window, 3, xpos=0, ypos=600, xsize=n_cols, ysize=n_rows
		window, 4, xpos=0, ypos=0, xsize=n_cols, ysize=n_rows
		loadct,0
	endif

	; set up regression fit params: x is array of model spectra ; -- removed(& line)
	;  unweighted (weightsi=1, use /relative_weights keyword)
	if ncomp GE 1 then tmp0 = transpose(reform(model(0,*)))
	if ncomp GE 2 then tmp1 = transpose(reform(model(1,*)))
	if ncomp GE 3 then tmp2 = transpose(reform(model(2,*)))
	if ncomp GE 4 then tmp3 = transpose(reform(model(3,*)))
	if ncomp GE 5 then tmp4 = transpose(reform(model(4,*)))
	if ncomp EQ 1 then x=[tmp0]  ;,transpose(line)]
	if ncomp EQ 2 then x=[tmp0,tmp1]   ; ,transpose(line)]
	if ncomp EQ 3 then x=[tmp0,tmp1,tmp2] ;,transpose(line)]
	if ncomp EQ 4 then x=[tmp0,tmp1,tmp2,tmp3]  ;,transpose(line)]
	if ncomp EQ 5 then x=[tmp0,tmp1,tmp2,tmp3,tmp4]  ;,transpose(line)]
	wts = make_array(n_E,value=1.0)
	bak = 1.  & inte = 1.
	count = 0
	for i=0,n_rows-1 do begin
  		spectrum=reform(tmp_ls.d(*,i))
;		if nzero_img GT 0 then groom1, spectrum, bak, inte, nzero_img
		fit=regress(x, spectrum, wts, yfit, c, z, z, z, z, chi2, status, /relative_weight)
		comp(*,i)=fit
		const(i)=c
		chisq(i)=chi2
		residual(*,i) = tmp_ls.d(*,i) - yfit
; ------- user feedback ... working ....
		count = count +1
		if count EQ 10 then begin
			text = 'Line_fit: point '+ strcompress(string(i)) +' of ' + strcompress(string(fix(n_rows)))
			axis_log, text
			count = 0
		endif
	endfor

; finished the fit - store and display residuals
	if keyword_set(axis) then begin
; -----------  store results in buffers and thumbnails
		for i = 0, ncomp-1 do begin
			tmp = {t:'1d', d: reform(comp(i,*)), x: img_y, xl: 'position (um)', dl: comp_names(i)}
			bad_index = where(finite(tmp.d) EQ 0, count)
			if count GT 0 then begin
				tmp.d(bad_index) = 0.
				axis_log, 'Line_fit: WARNING - NaN or Inf pixels set to 0.'
			endif
			if ncomp EQ 4 then offset = 5 else offset = 4
			CurBuf = i + offset
			HANDLE_VALUE, Data(Curbuf), tmp, /SET
			Plotbuf, CurBuf
			file = fileroot + comp_names(i) + '.txt'
			test = spc_save(tmp,file=file,  /noLbl)
		endfor
; ---------- store residual map in buffer 7
		tmp = {t:'2d', d: residual, x: img_x, y: img_y, xl: 'E (eV)',  yl: 'position',dl: 'residual'}
		bad_index = where(finite(tmp.d) EQ 0, count)
		if count GT 0 then begin
			tmp.d(bad_index) = 0.
			axis_log, 'Line_fit: WARNING - NaN or Inf pixels set to 0.'
		endif
		file = fileroot + 'res.axb'
		test = axb_save(tmp,file=file, /noLbl)
		CurBuf = 7
		HANDLE_VALUE, Data(CurBuf), tmp, /SET
		PlotBuf,CurBuf

; -------- constant term in buffer 8 --------
		tmp = {t:'1d', d: const, x: img_y, xl: 'position (um)', dl: 'constant '}
		bad_index = where(finite(tmp.d) EQ 0, count)
		if count GT 0 then begin
			tmp.d(bad_index) = 0.
			axis_log, 'Line_fit: WARNING - NaN or Inf pixels set to 0.'
		endif
		const_file = fileroot+'const.txt'
		test = spc_save(tmp,file=const_file, /noLbl)

; -------- chi square in buffer 9 -------
		tmp = {t:'1d', d: chisq, x: img_y, xl: 'position (um)', dl: 'chi-square '}
		bad_index = where(finite(tmp.d) EQ 0, count)
		if count GT 0 then begin
			tmp.d(bad_index) = 0.
			axis_log, 'Line_fit: WARNING - NaN or Inf pixels set to 0.'
		endif
		chi_file = fileroot + 'chi.txt'
		test = spc_save(tmp,file=chi_file, /noLbl)

; recall the constant and chi files to fill buffers 8 and 9 properly (seems like a kluge !!)
		tmp=spc_load(file=const_file)
		CurBuf = 8
		HANDLE_VALUE, Data(CurBuf), tmp, /SET
		PlotBuf, CurBuf

		CurBuf = 9
		tmp=spc_load(file=chi_file)
		HANDLE_VALUE, Data(CurBuf), tmp, /SET
		PlotBuf, CurBuf


	endif

endif
END

;_________________________________________

; ANCIENT procedure by rk 4/99 to groom one spectrum.
;  Outputs bkd-subtracted spectrum, background, and integral of line
pro groom1,spectrum,back,integral, nzero

; The first nzero+1 points are averaged to get the background
;nzero=50
back = total(spectrum(0:nzero))/(nzero+1)
spectrum = spectrum-back
integral=total(spectrum)
return
end