; Copyright (c) 1998-2009  A.P. Hitchcock  All rights reserved
;+
;NAME:
;  STACK_DIFF
;
;LAST CHANGED: ----------------------------------- 12-Jan-09
;
; PURPOSE:
;	This procedure differentiates the spectrum at every pixel of a stack
;
; CATEGORY:
;	stack processing; stand alone operation
;
; CALLING SEQUENCE:
;	STACK_DIFF, FILE+FILE
;
; INPUTS: none required
;
; KEYWORDS:
; 	FILE	name of the output file
;
; OUTPUTS: 	differentiated stack
;
; COMMON BLOCKS:
; AXIS_COM	standard set of common blocks
; stack_process_com
; BSIF_com
; volume_data, image_stack
;
; SIDE EFFECTS: none
;
; RESTRICTIONS: none
;
; MODIFICATION HISTORY:
; (12-Jan-09 aph) first version
;-

PRO stack_diff, file=file
@axis_com
@stack_process_com
COMMON volume_data, image_stack
@bsif_com

on_error,2

fltr = '*.ncb'

IF not keyword_set(file) then begin
	file = pickfile2(/READ, TITLE = 'STACK 1', FILTER=fltr, /LPATH, DEFPATH=defpath)
	if strlen(file) LE 0 THEN return  ; bail-out if no filename
ENDIF

stack_rb, file
stack = image_stack
t1 = size(image_stack)
;print, 'stack-1: ', t1

; ask if want to smooth before differentiating
test=DIALOG_MESSAGE('Smooth before differentiating ?', /QUESTION,/DEFAULT_NO)
IF test EQ 'Yes' then begin
	if widget_info(/active) EQ 1 then $		; assume if a widget is going, it is axis
           smth = get_num(Prompt = 'n-point smooth. n= (3,5,..)', val = 3, group=AXIS_ID) $
      else smth = get_num(Prompt = 'n-point smooth. n= (3,5,..)', val = 3)
; smooth if 3 or more
	smth_width= [1,1,smth]   ;this should only smooth in energy, not (x,y) dimensions
	image_stack = smooth(image_stack,smth_width,/edge_truncate, /NaN)
ENDIF

ax_wait		; turn on hourglass- this could be slow

; differentiate
;(from spectra~differentiate)
;	for i = 0, n_elements(tmp.d)-2 do begin
;		dax = tmp.x(i+1) - tmp.x(i)
;		if dax NE 0 then tmp.d(i) = (tmp.d(i+1) - tmp.d(i))/dax
;	endfor
;	tmp.d(n_elements(tmp.d)-1) = tmp.d(n_elements(tmp.d)-2)
;	tmp.d = shift(tmp.d,1)
;	tmp.d(0) = tmp.d(1)

for k = 0, n_cols-1 do begin
	for j = 0, n_rows-1 do begin
		y=reform(image_stack(k,j,*))
		for i = 0, n_elements(ev)-2 do begin
			dax = ev(i+1) - ev(i)
			if dax NE 0 then y(i) = (y(i+1) - y(i))/dax
		endfor
		y(n_elements(ev)-1) = y(n_elements(ev)-2)
		y = shift(y,1)
		y(0) = y(1)
		image_stack(k,j,*) = y
	endfor
endfor
ax_wait,/off		; turn off hourglass (sometimes sticks on)

; write out result
file = pickfile2(/READ, TITLE = 'name of differentiated stack', FILTER=fltr, /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

return
end

