; Copyright (c) 1998-2008 A.P. Hitchcock  All rights reserved
;+
;NAME:
;	  AX_POL_MAP
;
;LAST CHANGED: ----------------------------------- 18-sep-08
;
;PURPOSE:
;	This procedure fits a stack of angle dependent images
; to I(angle)=constant + amplitude*cos^2(angle-angle_Ref)
; with provision to force the reference angle within unser-set limits
;
;CATEGORY:
;	AXIS:  polarization stack analysis (operates stand-alone as well)
;
;CALLING SEQUENCE:
; for procedures:
;	AX_POL_MAP(/AXIS,/verbose)
;
;CALLED FROM AXIS:
;	->Stacks->maps~polarization fit
;
;INPUTS:
;	All input parameters are obtained through user dialogs.
;
;KEYWORDS:
;	AXIS - if set, indicates called from AXIS
;	VERBOSE - if set, displays results as generated line-by-line (SLOW!)
;
;ROUTINES
;	AX_POL_FIT - main routine
;
;OUTPUTS:
;	If run from AXIS, the results are placed in AXIS buffers
; buffer 4 = contant
; buffer 5 = amplitude
; buffer 6 = reference angle
; buffer 7 - Standard deviation of the constant
; buffer 8 - Standard deviation of the amplitude
; buffer 9 - Standard deviation of the ref. angle
;
;COMMON BLOCKS:
;	@AXIS_COM	standard set of common blocks
;	@stack_process_com
;	Volume_data - STACK
; 	@BSIF_COM - nc data
;
;PROCEDURE:
;	Uses ax_pol_fit applied to the curve at each pixel
;
;MODIFICATION HISTORY:
; ***********
; (04-May-06 aph) first developed from stack_fit.pro
; (19-Feb-07 aph) make tolerance a parameter; report no-change parameters
; (18-mar-07 aph) add limits to ref_angle
; (17-Sep-08 aph) replace analcom with stack_process_com
;-


PRO ax_pol_map, axis=axis, verbose = verbose
@axis_com
@stack_process_com
COMMON volume_data, image_stack
@bsif_com

on_error,2


; ------------ read-in stack data -----------------------
file = pickfile2(/READ, FILTER='*.ncb', /LPATH)
if strlen(file) EQ 0 THEN RETURN  ; bail-out if no filename

stack_rb, file		; read in stack in image_stack (Volume_data), energies in ev (stack_process_com)
t=ax_name(file)
fileroot =  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 = t(0) + fileroot + '_'

; ------ select starting parameters -------------
t=size(image_stack)
nx = t(1)  & ny = t(2)
n_angles=t(3)
axis_log, 'nx =' + strtrim(string(fix(nx)),2)
axis_log, 'ny =' + strtrim(string(fix(ny)),2)
axis_log, 'n_angles =' + strtrim(string(fix(n_angles)),2)
param=fltarr(3)
param(0) = 0.5
param(1) = 0.5
param(2)= 90.
;---- check if wish to use a fit parameter file -------------------
if keyword_set(axis) then begin
	t = dialog_message('Display average for definition of starting parameters ?', /question)
	if t(0) EQ 'Yes' then begin
		avg = fltarr(n_angles)
		for i = 0, nx-1 do begin
			for j = 0, ny-1 do begin
				avg = avg + image_stack(i,j,*)
			endfor
		endfor
		avg = avg/(nx*ny)
		tmp={t:'1d', x: ev, d:avg, xl:'angle(deg)', dl:'avg of pol. dep. stack'}
	 	HANDLE_value, Data(0),tmp,/set
	   	Plotbuf,0
		param(0) = min(avg, max=max_avg)
		param(1) = max_avg-param(0)
		param(2)= 90.
	endif
endif

; -------- acquire starting values from user
if keyword_set(axis) then $
    param(0) = get_num(prompt='constant', val=param(0), group=axis_ID) $
    else param(0) = get_num(prompt='constant', val=param(0))
if keyword_set(axis) then $
	param(1) = get_num(prompt='amplitude', val=param(1), group=axis_ID)$
	else param(1) = get_num(prompt='amplitude', val=param(1))
if keyword_set(axis) then $
	param(2) = get_num(prompt='ref_angle', val=param(2), group=axis_ID) $
	else param(2) = get_num(prompt='ref_angle', val=param(2))

; ------- apply clip to ref_angle if desired ----
limits_on = dialog_message('Clip reference angle at limits ?', /question)
if limits_on(0) EQ 'Yes' then begin
	limits = fltarr(2)
	if keyword_set(axis) then $
		 limits(0) = get_num(prompt='lower limit', val = limits(0), group = axis_ID) $
		   else limits(0) = get_num(prompt='lower limit',  val = limits(0))
	if keyword_set(axis) then $
		limits(1) = get_num(prompt='upper limit',  val = limits(1), group = axis_ID) $
		   else limits(1) = get_num(prompt='upper limit',  val = limits(1))
	if limits(0) GT limits(1) then begin ; protect user
		t = limits(0)
		limits(0)=limits(1)
		limits(1) = t
	endif
	print, 'Force ref_angle between ', limits(0),' and ', limits(1)
endif

; ------- set tolerance for fits
tol = 1e-6
if keyword_set(axis) then $
	tol = get_num(prompt='tolerance', val=tol, group=axis_ID) $
	else tol = get_num(prompt='tolerance', val=tol)

; ---------- get max iterations
max_iter=200
if keyword_set(axis) then $
	max_iter = get_num(prompt='maximum # of iterations', val=max_iter, group=axis_ID) $
	else max_iter = get_num(prompt='maximum # of iterations', val=max_iter)

; ------ execute polarization fit on pixel-by-pixel basis
const_map = fltarr(nx,ny)
ampl_map = fltarr(nx,ny)
Refang_map = fltarr(nx,ny)
stdDev_const_map = fltarr(nx,ny)
stdDev_ampl_map = fltarr(nx,ny)
stdDev_refang_map = fltarr(nx,ny)
curv=fltarr(2,n_angles)
curv(0,*) = ev
n_same_angle = 0
n_same_amplitude = 0
n_same_constant = 0
n_successful = 0
n_divergent = 0
n_nonconvergent = 0
deviation = 0.01
WIDGET_CONTROL, /Hourglass
n_low=0 & n_hi=0
for j = 0, ny-1 do begin
	for i = 0, nx-1 do begin
		par=param
		curv(1,*) = image_stack(i,j,*)
		polfit = Ax_Pol_Fit(curv,param=par, sigma=sigma, max_iter=max_iter, status=status, tol=tol, /silent)
;		print, status
		const_map(i,j)=par(0)
		ampl_map(i,j)=par(1)
; --------- apply limits to ref_angle if requested
		if limits_on(0) EQ 'Yes' then begin
			if par(2) LT limits(0) then begin
				par(2)=limits(0)
				n_low = n_low+1
			endif
			if par(2) GT limits(1) then begin
				par(2)=limits(1)
				n_hi=n_hi+1
			endif
		endif
		Refang_map(i,j)=par(2)
		stddev_const_map(i,j)= total(sigma(0))/float(n_angles)
		stddev_ampl_map(i,j)= total(sigma(1))/float(n_angles)
		stddev_refang_map(i,j)= total(sigma(2))/float(n_angles)
		if status EQ 0 then n_successful = n_successful + 1
		if status EQ 1 then n_divergent = n_divergent + 1
		if status EQ 2 then n_nonconvergent = n_nonconvergent + 1
		if abs(par(0) - param(0)) LT deviation*param(0) $
			 then  n_same_constant = n_same_constant + 1
		if abs(par(1) - param(1)) LT deviation*param(1) $
			 then   n_same_amplitude = n_same_amplitude + 1
		if abs(par(2) - param(2)) LT deviation*param(2) $
			 then  n_same_angle = n_same_angle + 1
	endfor
; ------- user feedback ... working ....
;	if keyword_set(AXIS)THEN
	axis_log, 'line ' + strcompress(string(fix(j)))
endfor

; ----------- output section -------------

; How good was the fit ??

if limits_on(0) EQ 'Yes' then begin
	text = 'Ref. angle limits applied. # low = ' + strcompress(string(fix(n_low)))
	text = text + '. # high = ' +  strcompress(string(fix(n_hi)))
	axis_log, text
endif
text = 'unchanged constant' +  strcompress(string(fix(n_same_constant))) + ' of ' + strcompress(string(fix(nx*ny)))
axis_log, text
text = 'unchanged amplitude' + strcompress(string(fix(n_same_amplitude))) + ' of ' + strcompress(string(fix(nx*ny)))
axis_log, text
text = 'unchanged angle' +  strcompress(string(fix(n_same_angle))) + ' of ' + strcompress(string(fix(nx*ny)))
axis_log, text
axis_log, 'of ' + strcompress(string(fix(nx*ny))) + ' pixels fitted'
text = strcompress(string(fix(n_successful))) + ' were successful'
axis_log, text
text = strcompress(string(fix(n_divergent))) + ' were divergent'
axis_log, text
text = strcompress(string(fix(n_nonconvergent))) + ' did not converge'
axis_log, text

; ------------ set up (x,y) scales from bsif_com information
x_step = (x_stop - x_start)/n_cols
y_step = (y_stop - y_start)/n_rows
img_x = x_start + findgen(n_cols)*x_step
img_y = y_start + findgen(n_rows)*x_step


; finished the fit - store and display residuals
if keyword_set(axis) then begin
; ----------- replace NaN values with 0.
	tmp = {t:'2d', d: const_map, x: img_x ,y: img_y, E: 0.0, xl:'X',yl: 'Y',  dl: 'constant '}
	bad_index = where(finite(tmp.d) EQ 0, count)
	if count GT 0 then begin
		tmp.d(bad_index) = 0.
		print, 'ax_pol_fit: WARNING - NaN or Inf pixels set to 0.'
	endif
	HANDLE_VALUE, Data(4), tmp, /SET
	file = fileroot+'const.axb'
;	overwrite_all = 0   ; reset overwrite flag in case was set by earlier uses
	test = axb_save(tmp,file=file)
	PlotBuf,4
; ----- store and plot amplitudes -----------
; ----------- replace NaN values with 0.
	tmp = {t:'2d', d: ampl_map, x: img_x ,y: img_y, E: 0.0, xl:'X',yl: 'Y',  dl: 'pol. amplitude '}
	bad_index = where(finite(tmp.d) EQ 0, count)
	if count GT 0 then begin
		tmp.d(bad_index) = 0.
		print, 'polarization_fit: WARNING - NaN or Inf pixels set to 0.'
	endif
	HANDLE_VALUE, Data(5), tmp, /SET
	file = fileroot+'ampl.axb'
;	overwrite_all = 0   ; reset overwrite flag in case was set by earlier uses
	test = axb_save(tmp,file=file)
	PlotBuf,5
; ------------------ store reference angle maps ------------
; ----------- replace NaN values with 0.
	tmp = {t:'2d', d: refang_map, x: img_x ,y: img_y, E: 0.0, xl:'X',yl: 'Y',  dl: 'ref. angle '}
	bad_index = where(finite(tmp.d) EQ 0, count)
	if count GT 0 then begin
		tmp.d(bad_index) = 0.
		print, 'ax_pol_fit: WARNING - NaN or Inf pixels set to 0.'
	endif
	HANDLE_VALUE, Data(6), tmp, /SET
	file = fileroot+'angle.axb'
;	overwrite_all = 0   ; reset overwrite flag in case was set by earlier uses
	test = axb_save(tmp,file=file)
	PlotBuf,6
; -------- store errors
	tmp = {t:'2d', d: stddev_const_map, x: img_x ,y: img_y, E: 0.0,  xl: 'X', yl: 'Y',  dl: 'stddev_const '}
	bad_index = where(finite(tmp.d) EQ 0, count)
	if count GT 0 then begin
		tmp.d(bad_index) = 0.
		print, 'ax_pol_fit: WARNING - NaN or Inf pixels set to 0.'
	endif
	tmp = ax_order(tmp)
	HANDLE_VALUE, Data(7), tmp, /SET
	PlotBuf,7
	tmp = {t:'2d', d: stddev_ampl_map, x: img_x ,y: img_y, E: 0.0,  xl: 'X', yl: 'Y',  dl: 'stddev_ampl '}
	bad_index = where(finite(tmp.d) EQ 0, count)
	if count GT 0 then begin
		tmp.d(bad_index) = 0.
		print, 'ax_pol_fit: WARNING - NaN or Inf pixels set to 0.'
	endif
	tmp = ax_order(tmp)
	HANDLE_VALUE, Data(8), tmp, /SET
	PlotBuf,8
	tmp = {t:'2d', d: stddev_refang_map, x: img_x ,y: img_y, E: 0.0,  xl: 'X', yl: 'Y',  dl: 'stddev_refang '}
	bad_index = where(finite(tmp.d) EQ 0, count)
	if count GT 0 then begin
		tmp.d(bad_index) = 0.
		print, 'ax_pol_fit: WARNING - NaN or Inf pixels set to 0.'
	endif
	tmp = ax_order(tmp)
	HANDLE_VALUE, Data(9), tmp, /SET
	PlotBuf,9
	CurBuf = 9
ENDIF
END
