; Copyright (c) 1998, ALS SpectroMicroscopy Facility, jd
; Copyright (c) 1998-2009 A.P. Hitchcock  All rights reserved
;+
;NAME:
;		LOADSPEM
;
;LAST CHANGED: ----------------------------------- 	01-Aug-09
;
; PURPOSE:
;	This procedure reads a file in XPD Labview format
;  written by early versions (pre-2000) of the ALS BL 7.0 SPEM
;
; CATEGORY:
;	Image display.
;
; CALLING SEQUENCE:
;	Result = LOADSPEM(fname, PK=pknum, NORM=normnum, AVG=avgnum, $
;                  PLOT=pl, GLITCH=ngli, SDEV=nsdev, QUIET=quiet)
;
; INPUTS:
; FNAME - filename

; KEYWORDS:
;	PK
;	NORM
;	AVG
;	PLOT	 - display (splot2d) image (first image if multiple)
;	GLITCH	 - run manual glitch threshold window
;	SDEV=n   - run auto glitch filter for data outside n std deviations
;	QUIET
;
; OUTPUTS:
;	No explicit outputs.   A new plot window is created if necessary
;
; COMMON BLOCKS:
;	AXIS_COM	standard set of common blocks
;
; SIDE EFFECTS:
;	A window is created if plot is selected.
;
; MODIFICATION HISTORY:
; (07-May-99 jdd) version taken over from Jonathan Denlinger
; (15nov00 cgz) Replaced call of DIMEN function by explicit (SIZE(var))(index+1)
; (15nov00 cgz) Migrated to RD_SPEM.PRO since it was only called from with RD_SPEM.PRO
;			LOADSPEM.PRO is now obsolete and can be eliminated
; (30-nov-00 cgz) Replaced call of nwin with explicit window statement
; (01-Aug-09 aph) AXIS standard header added
; (11-Jan-10 aph) commented out all lines calling GLITCH - can not find in package
;                 most likely in one of Eli rotenbergs' combined utility files
;-

function loadspem, fname, PK=pknum, NORM=normnum, AVG=avgnum, $
                  PLOT=pl, GLITCH=ngli, SDEV=nsdev, QUIET=quiet

on_error,2
@axis_com
  ;if (keyword_set(fname) eq 0) then fname=pickfile2()

if not keyword_set(quiet) then quiet=0
f=gen_read(fname, Q=quiet)
fileshort=strmid(f.fname,strlen(f.fname)-9, strlen(f.fname))

xx=reform(f.indvars(0,*))
;ncol=dimen(f.indvars,1)
ncol=(size(f.indvars))(2)

;determine repeat in 1st column  (slow loop)
xx0=xx(0)
nx=1
for i=1,ncol-1 do begin
   if (xx(i) eq xx0) then begin
      nx=i
      goto, break
   endif
endfor

break:
ny=floor(ncol/nx)
nc=nx*ny
;nd=dimen(f.singledata,0)
nd=(size(f.singledata))(1)
print, ncol, nc, nx, 'x', ny, ': nd=', nd

xaxis=xx(0:nx-1)
yy=reform(transpose(reform(f.indvars(1,0:nc-1),nx, ny)))
yaxis=yy(0:ny-1)

;image(s)
;if n_elements(pknum) eq 0 then pknum=0
if n_elements(avgnum) gt 1 then begin   ;average over indices in array
   navg=n_elements(avgnum)
   print, 'Avg Peak over ', navg, ' indices'
   img=fltarr(nx, ny)
   for i=0,navg-1 do img=img+reform(f.singledata(avgnum(i),0:nc-1),nx, ny)
   img=img/navg
   npk=1
endif else begin
   if n_elements(pknum) ge 1 then begin   ;peak number(s) specified
      npk=n_elements(pknum)
      img=fltarr(nx, ny, npk)
      for i=0,npk-1 do img(*,*,i)=reform(f.singledata(pknum(i),0:nc-1),nx, ny)
   endif else begin           ;get all images
      npk=f.summary[4]    ;or n_elements( f.singlenames )
      img=fltarr(nx, ny, npk)
      for i=0,npk-1 do img(*,*,i)=reform(f.singledata(i,0:nc-1),nx, ny)
   endelse
   img=reform(img)
endelse

;  s={d:dat, x:xx, y:yy, xl:f.indnames(0), yl:f.indnames(1), dl:fileshort}
;print, keyword_set(norm), n_elements(norm)
norm=0
if n_elements(normnum) gt 0 then begin
	norm=1
	if n_elements(normnum) gt 1 then begin
	   ;average over indices in array
	   nnorm=n_elements(normnum)
	   print, 'Avg Norm over ', npk, ' indices'
	   nrm=fltarr(nx, ny)
	   for i=0,nnorm-1 do nrm=nrm+reform(f.singledata(normnum(i),0:nc-1),nx, ny)
	   nrm=nrm/nnorm
	endif else begin
	   normnum=max([min([normnum,nd-1]),0])
	   nrm=reform(f.singledata(normnum,0:nc-1),nx, ny)
	endelse
endif

if (npk gt 1) then begin
   if norm then $
      s={d:img, dn:nrm, x:xaxis, y:yaxis, z:indgen(npk), xl:'x (um)', yl:'y (um)', zl:'BE', sl:fileshort, dl:''} $
   else $
      s={d:img, x:xaxis, y:yaxis, z:indgen(npk), xl:'x (um)', yl:'y (um)', zl:'BE', sl:fileshort, dl:''}

endif else begin
   if norm then $
      s={d:img, dn:nrm, x:xaxis, y:yaxis, xl:'x (um)', yl:'y (um)',sl:fileshort, dl:''} $
   else $
      s={d:img, x:xaxis, y:yaxis, xl:'x (um)', yl:'y (um)', sl:fileshort, dl:''}
endelse

;  if keyword_set(pknum) eq 0 then pknum=0
;;  if keyword_set(nps) eq 0 then nps=max([dimen(xx,0),dimen(yy,0)])
;  if keyword_set(nps) eq 0 then nps=max([(size(xx))(1),(size(yy))(1)])
;  if keyword_set(nm) then $
;    zz=reform(f.singledata(pknum,*)/f.norm1dim(pknum,*)) $
;  else $
;    zz=reform(f.singledata(pknum,*))
;  triangulate,xx,yy,tri
;  ans={d:trigrid(xx,yy,zz,tri,abs([(min(xx)-max(xx))/(nps-1), (min(yy)-max(yy))/(nps-1)])),$
;       x:xx, y:yy, xl:f.indnames(0), yl:f.indnames(1)}

; ------------ 11-Jan-2010 -----COMMENTED OUT as do not have GLITCH procedure
; if keyword_set(ngli) then GLITCH, s, ngli
; if keyword_set(nsdev) then GLITCH, s, SDEV=nsdev

if keyword_set(pl) then begin
   log=0
   if pl eq 2 then log=1
;   if !d.window eq -1 or !d.window gt 31 then nwin, 0
   if !d.window eq -1 or !d.window gt 31 then window, 0, xsize=300, ysize=300
;   if !d.window eq 20 and keyword_set(ngli) then nwin		; 20=glitch filter window
   if !d.window eq 20 and keyword_set(ngli) then $
         window, !d.window+1<0, xsize=300, ysize=300		; 20=glitch filter window
		; cgz - don't think there is any possibility that !d.window+1 could be less than zero
		;		i.e., no window can ever be anything other than -1, 0, 1, 2, ...
   splot2d, s, NORM=norm, LOG=log, _extra=e, charsize=1
endif
return, s
end
