; Copyright(c) 2008 A.P. Hitchcock  All rights reserved
;+
;NAME: READ_FTIR
;
;LAST CHANGED: ----------------------- 28-Dec-08 (aph)
;
; PURPOSE:
;   This prodecure reads an FTIR map (written by Thermo Omnic program)
;  consisting of IR spectra at an array of points
;  The data are converted to an aXis binary stack (*.ncb)
;
; CATEGORY:
;	utility; stand-alone or used through aXis2000
;
; CALLING SEQUENCE:
;READ_FTIR, file=file, filter=filter, xcol=xcol, ycol=ycol, $
;         multi=multi, skip=skip, int_factor=int_factor, typ=typ, NOF=nof, _extra=e
;
; INPUTS: only by keyword
;
; KEYWORDS:
; FILE		filename (*.map - array dimensions)
; FILTER  	filter extension (*.map)
; SKIP		number of header lines to skip  (if multi-column)
; MULTI		returns number of columns -1 = number of spectra
; NOF		no filter (*.*)
; VERBOSE	if set, print information as extra data
; _EXTRA	pass through keywords
;
; OUTPUTS:
;	Result is  an aXis2000 stack (x,y,wavenumber) (*.nc, *.dat files)
;   and a text file, (*.par) listing the parameters of the FTIR image acquisition
;
; COMMON BLOCKS:
;	AXIS_COM	        standard set of aXis2000 common blocks
;   STACK_PROCESS_COM   common for stacks
;   VOLUME_DATA         data array of stack
;   BSIF_COM            common for stack acquisition parameters
;
; MODIFICATION HISTORY:
; (18-dec-08 aph) first wrote a code for set of *.csv data sent by Paul Dumas (Soleil)
;            IRMAP_LOAD  - kept in case sets of csv are encountered from FTIR sources
; (26-Dec-08 aph) generated READ_FTIR code from nicolet.m  Matlab file (issues: IDL arrys start at 0, Matlab at 1)
; (28-Dec-08 aph) add choice of wavenmumber range to extract image(s)

FUNCTION READ_FTIR, file=file, filter=filter, verbose=verbose, $
           skip=skip, NOF=nof, int_factor=int_factor, multi=multi, _extra=e
@axis_com
@stack_process_com
COMMON volume_data, image_stack
@bsif_com
on_error,2

s=0
fltr='*.map'
if keyword_set(nof) then fltr='*'
if keyword_set(filter) then fltr=filter
if not keyword_set(file) then $
	file=pickfile2(/read, title='Read FTIR map file', filter='*.map', /LPATH)
if strlen(file) LE 0 THEN return, s    ; bail-out if no filename
t = ax_name(file)
path = t(0)
root = t(1)
openr, lun, file, /get_lun

axis_log, 'READ_FTIR: Thermo *map (v7) -> aXis2000 stack & par files'

; --- determine size & read-in all file
t =fstat(lun)
map_length = t.size/4
;print, 'length of file (in 32-bit reals = float*32)',map_length

mf=fltarr(map_length)
readu, lun, mf

; ------- skip 93*4 8-bit unsigned integers (uint8 in Matlab, byte in IDL)
POINT_LUN, lun, 372	; position map file to start of info block
tmp=ulong(1)
readu, lun, tmp
indicblock = (tmp-204)/4
UWN = mf(indicblock+3);   ; upper wavenumber limit
LWN = mf(indicblock+4);   ; lower wavenumber limit
BGG = mf(indicblock+13);  ; gain background
IDF = mf(indicblock+17);  ; identifier for start indices of spectra
LAS = mf(indicblock+19);  ; Laser wavenumber
APT = mf(indicblock+22);  ; aperture size spectrometer
HPF = mf(indicblock+39);  ; high pass filter
LPF = mf(indicblock+40);  ; low pass filter
VEL = mf(indicblock+46);  ; mirror velocity
if keyword_set(verbose) then begin
	t_start = systime(/seconds)
	print, ' '
	print, 'FTIR map being converted from file ', root+'.map' ,' STARTED: ', systime()
	print, '------------------------------'
	print, 'pointer= ', tmp
	print, 'indicblock = ', indicblock
	print, 'upper wavenumber= ', uwn
	print, 'lower wavenumber= ', lwn
	print, 'gain background= ', bgg
	print, 'spectral start index= ', idf
	print, 'Laser wavenumber= ', las
	print, 'aperture size= ', apt
	print, 'high pass filter= ', hpf
	print, 'low pass filter= ', lpf
	print, 'mirror velocity= ', vel
endif

POINT_LUN, lun, indicblock*4	; position map file to start of info block
mint= ulonarr(13)
readu, lun, mint
NOD = Mint(0);             ;  NOD - number of data points
NMP = Mint(6);             ;  NMP - number of measurement points
IPP = Mint(7);             ;  IPP - interferogram peak position
NSS = Mint(8);             ;  NSS - number of sample scans
NPF = Mint(10);            ;  NPF - number of points for FFT
NGS = Mint(12);            ;  NGS - number of background scans
WVS  = (UWN-LWN)/(NOD-1);    ;  WVS - wavenumber step
if keyword_set(verbose) then begin
	print, ' '
	print, 'parameters of each FTIR spectrum'
	print, '--------------------------------'
	print, '# of data points = ', nod
	print, '#  of measurement points = ', nmp
	print, '# of sample scans = ', nss
	print, '# of points for FFT = ', npf
	print, '# of background scans = ', ngs
	print, 'interferogram peak position = ', ipp
	WVS  = (UWN-LWN)/(NOD-1)
	print, 'wavenumber step = ', wvs
endif

; ------- get information about map -----------
POINT_LUN, lun, (indicblock+249)*4+2
mapprm = fltarr(6)
readu, lun, mapprm
Ystart = mapprm(0);
Ystop  = mapprm(1);
STY    = mapprm(2);           ; STY - step size in y direction
Xstart = mapprm(3);
Xstop  = mapprm(4);
STX    = mapprm(5);           ; STX - step size in x direction
SZY    = Ystop-Ystart;        ; SZY - map size in y-direction
SZX    = Xstop-Xstart;        ; SZX - map size in x-direction
YDI    = round(SZY/STY+1);    ; YDI - number of spectra in y-direction
XDI    = round(SZX/STX+1);    ; XDI - number of spectra in x-direction
XxY    = XDI * YDI;           ; XxY - number of spectra (equals xdim  * ydim)
if keyword_set(verbose) then begin
	print, ' '
	print, 'parameters of map'
	print, '-----------------'
	print, 'X-step = ', stx
	print, 'Y-step = ', sty
	print, 'map size = ', szx, ' by ', szy
	print, 'map dimensions = ', xdi, ' by ', ydi
	print, 'total number of spectra = ', xxy
endif

; ------- locate start of spectral data block
t=where(mf EQ idf, count)
if count LE 0 then begin
	print, 'cannot find start of spectral block = ', idf
endif

prog = 0
spcstart = t(1)+31
ctr = lonarr(2)
npts = nod-1				; # of points per spectrum
c = fltarr(npts, xxy)
for i=0, xxy-1 do begin
      ctr(0) = spcstart-NOD + (i+1)*(NOD+25)-1
      ctr(1) = ctr(0)+ NOD-2
;      print, ' spectrum ', i, ' from ', ctr(0),' to ', ctr(1)
      c(*,i) = mf(ctr(0):ctr(1))
; special patch for Michael C. Martin: checking whether the spectrum is a zero line spectrum
     if total(C(*,i)) EQ 0 then begin
      	tr = randomn(0.5, NOD)
      	C(*,i) = tr*0.0001
      endif
      progr   = (i+1)/XxY
      if (progr - prog) GT 0.01 then prog = progr
endfor

close, lun & free_lun, lun

; prepare * write aXsi2000 stack file
c = transpose(c)
; -------- Create aXis2000 stack
ncols = xdi & nrows = ydi
image_stack=reform(c, ncols, nrows, npts)

ev = lwn + wvs*findgen(npts)

filename_ev_msec_list = strarr(npts)
for i = 0, npts-1 do begin
	filename_ev_msec_list(i) = root + ' ' + string(ev(i))
endfor
x_start = Xstart  & x_stop = Xstop & x_step = stx
y_start = ystart  & y_stop = ystop & y_step = sty

; --------- ask user what file name they want - default name is 'root'
sname=pickfile2(/write, file= root +'.ncb', filter='*.ncb', title = 'name of binary stack file')
stack_wb, sname
s=1

;  ---------- write information file to accompany stack
; "*.par" file contains map info variables
t=ax_name(sname)
info_name =  t(0)+t(1)+'.par'
info_name = pickfile2(/write, file= info_name , filter='*.par', title = 'name of parameter file')
openw, lun, info_name, /get_lun
printf, lun, 'FTIR map converted from file ', root+'.map' ,' on ', systime()
printf, lun, ' '
printf, lun, 'Parameters reported in file'
printf, lun, '---------------------------'
printf, lun, UWN, ' ; upper wavenumber limit'
printf, lun,  LWN, ' ;  lower wavenumber limit'
printf, lun,  BGG, ' ;  gain background'
printf, lun,  IDF, ' ;  identifier for start indices of spectra'
printf, lun,  LAS, ' ;  Laser wavenumber'
printf, lun,  APT, ' ;  aperture size spectrometer'
printf, lun,  HPF, ' ;  high pass filter'
printf, lun,  LPF, ' ;  low pass filter'
printf, lun,  VEL, ' ;  mirror velocity'
printf, lun,  NOD, ' ;  number of data points'
printf, lun,  NMP, ' ;  number of measurement points'
printf, lun,  IPP, ' ;  interferogram peak position'
printf, lun,  NSS, ' ;  number of sample scans'
printf, lun,  NPF, ' ;  number of points for FFT'
printf, lun,  NGS, ' ;  number of background scans'
printf, lun,  WVS, ' ;  wavenumber step'
printf, lun,  Xstart, ' ; start X-position'
printf, lun,  Ystart, ' ; start Y-position'
printf, lun,  Ystop,  ' ; end   X-position'
printf, lun,  Xstop,  ' ; end   X-position'
printf, lun,  STY, ' ;  step size in y direction'
printf, lun,  STX, ' ;  step size in x direction'
printf, lun,  SZY, ' ;  map size in y-direction'
printf, lun,  SZX, ' ;  map size in x-direction'
printf, lun,  YDI, ' ;  number of spectra in y-direction'
printf, lun,  XDI, ' ;  number of spectra in x-direction'
printf, lun,  XxY, ' ;  number of spectra (equals xdi  * ydi)'
;
close, lun & free_lun, lun
t=ax_name(info_name)
axis_log, 'Parameters in ' + t(1)+ '.' + t(2)

if keyword_set(verbose) then begin
      print, 'FTIR map conversion END: ', systime()
      print, 'conversion time (sec) = ', systime(/seconds) - t_start
endif


; determine if there is an active widget, if so assume it is aXis2000
axis_on = widget_info(/active)
; generate average image and return as a 2d image
min_wnm = 1600
if axis_on NE 0 then begin
	 min_wnm = get_num(prompt='image sum: lower wavenumber', val=min_wnm, group=axis_ID)
endif else min_wnm = get_num(prompt='image sum: lower wavenumber ', val=min_wnm)
max_wnm = 1720
if axis_on NE 0 then begin
	 max_wnm = get_num(prompt='image sum: upper wavenumber', val=max_wnm, group=axis_ID)
endif else max_wnm = get_num(prompt='image sum: upper wavenumber ', val=max_wnm)
if max_wnm-min_wnm LE ev(1)-ev(0) then begin
	ind_wnm=where(ev GE min_wnm, count)				; SINGLE IMAGE
	if count GT 0 then begin
		img=image_stack(*,*,ind_wnm(0))
		e_wnm = ev(ind_wnm(0))
	endif else begin
		 img=image_stack(*,*,255)				; in case selected wavenumber is out of range
		 e_wnm = ev(255)
	endelse
	xl = 'um Image at ' + string(e_wnm, format='(I5)') + ' cm-1'
endif else begin
	ind_wnm=where(ev GE min_wnm AND ev LE max_wnm, count)
	if count GT 0 then begin						;SUM OF IMAGES
		t=image_stack(*,*,ind_wnm)
		e_wnm = mean([min_wnm,max_wnm])
	endif else begin								; ALL IMAGES
		t=image_stack
		e_wnm = 0.001
		ind_wnm=n_elements(ev)
	endelse
	img= total(t,3, /nan)
	xl = 'um    Avg from ' + string(ev(ind_wnm(0)), format='(I5)') + ' to '
	xl = xl + string(ev(ind_wnm(n_elements(ind_wnm)-1)),format='(I5)') + ' cm-1'
endelse
x= xstart+findgen(ncols)*x_step
y= ystart+findgen(nrows)*y_step

s = {t:'2d', x:x, y:y, E: e_wnm, d:img, xl:xl, yl:'  um', dl: 'average image of '+ root}
; ---- put spectrum in buffer 0
if axis_on NE 0 then begin
	d = fltarr(npts)
	for i = 0, npts-1 do d(i) = total(image_stack(*,*,i))

	xl = 'cm-1   Average over whole image'
	yl = 'FTIR spectrum'
	dl = 'FTIR spectrum of ' + root
	tmp = {t:'1d', x:ev, d:d, dn:d, xl: xl, yl:yl, dl: dl}
	HANDLE_VALUE, Data(0), tmp, /set
	Label(0) = tmp.dl
	PlotBuf,CurBuf
endif

return, s

end
