; Copyright (c) 1998-2017 A.P. Hitchcock  All rights reserved
;+
;NAME:
;  STACK_Illumination_Correct
;
;LAST CHANGED: ----------------------------------- 28 Sep 2017
;
; PURPOSE:
;	This function corrects for uneven vertical illumination on a line-by-line basis
; the illumination along columns as f(energy) is returned as an image
;
; CATEGORY:
;	stack processing; stand alone operation or from axis2000
;
; CALLING SEQUENCE:
;	RESULT = STACK_Illumination_Correct()
;
; INPUTS: none
;
; KEYWORDS: none
;
; OUTPUTS: 	routine stores the modified stack to disk
;
; COMMON BLOCKS:
; AXIS_COM	standard set of common blocks
; stack_process_com
; BSIF_com
; volume_data, image_stack
;
; RESTRICTIONS: none
;
; MODIFICATION HISTORY:
; (28-Sep-17 aph) first version -  needed for cryo-STXM in cases where OSA not scanned vertically
;-

FUNCTION stack_illumination_correct
@axis_com
@stack_process_com
COMMON volume_data, image_stack
@bsif_com

on_error,2

; ---- read in illumination stack (a truncated version of the full stack containing the illumination change
zero = 0
stack_ill   = pickfile2(/READ, title='select illumination stack', FILTER='*.ncb',/LPATH, DEFPATH=defpath)
if strlen(stack_ill) EQ 0 then return, zero
stack_rb, stack_ill
ill=image_stack
ts = size(ill)
ill_ncols=ts(1)
ill_nrows=ts(2)
ill_npts=ts(3)
print, 'Ill: ncols, nrows, npts ', ill_ncols, ill_nrows, ill_npts

; ---- read in stack to be corrected ----
stack_data = pickfile2(/READ, title='select stack to correct', FILTER='*.ncb',/LPATH, DEFPATH=defpath)
if strlen(stack_data) EQ 0 then return, zero
stack_rb, stack_data
stack=image_stack
ts = size(stack)
stack_ncols=ts(1)
stack_nrows=ts(2)
stack_npts=ts(3)
print, 'STACK: ncols, nrows, npts ', stack_ncols, stack_nrows, stack_npts

; ------ check for data file compatibility ---
if ill_nrows NE stack_nrows or ill_npts NE stack_npts then begin
	axis_log, '# of rows or # of energies of Illumination stack  differ from those of the transmission stack'
endif

; ----- generate illumination image-spectrum by averaging each row of the Io stack ---
Ill_image = fltarr(stack_npts, stack_nrows)
for i = 0, stack_npts - 1 do begin
	for j = 0, ill_nrows - 1 do begin
		t =  moment(ill(*,j,i))
;		print, j, '  ', t(0)
		Ill_image(i,j) = t(0)
	endfor
; ------ modify Ill_image to be centered about 1.0 at each energy
     t = moment(Ill_image(i,*))
     Ill_image(i,*) = Ill_image(i,*)/t(0)
endfor

print, 'Ill_image size: (n_energies, nrows) ', stack_npts, stack_nrows

; generate axis format image  of the Ill signal, normalized to centre about 1
y = findgen(stack_nrows)
nx = ill_nrows
x = ev
t=ax_name(stack_data)
text = t(1) + ' illumination vertical at each E'
tmp = {t:'2d', x:x, y:y, xl:'Energy (eV)', e: ev(0), yl:'rows', d: Ill_image , dl:text}
help, tmp, /struct

; ------- generate illumination corrected stack
image_stack = fltarr(stack_ncols, stack_nrows, stack_npts)
for k = 0, stack_npts -1 do begin
	for i = 0, stack_ncols - 1 do begin
		for j = 0, stack_nrows - 1 do image_stack(*,j,k)=stack(*,j,k)*Ill_image(k,j)
	endfor
endfor

; ----- write out converted stack -----------------
file = pickfile2(/READ, TITLE = 'Filename for illumination corrected stack', FILTER='*.ncb', /LPATH, DEFPATH=defpath)
if strlen(file) NE 0 then begin
	t = ax_name(file)
	filename = string(t(0)) + string(t(1)) + '.ncb'
	stack_wb, filename
endif

ax_wait,/off		; turn off hourglass (sometimes sticks on)
return, tmp

end

