; Copyright (c) 1998-2013 A.P. Hitchcock  All rights reserved
;+
;NAME:
;		XAS2.PRO
;
;LAST CHANGED: ----------------------------------- 	30-Aug-2013
;
; PURPOSE:
;	This function reads two Beamline 7 STXM X-ray absorption data files
;  and computes Beers- law absorption as ln(hole/DET)
;
; CATEGORY:
;	read in
;
; CALLING SEQUENCE:
;	Result = XAS2 (file, PLOT=pl, DEFPATH=defpath, OPT=opt, _extra=e )
;
; INPUTS:
;	FILE	filename1
;
; KEYWORDS:
;
; OPT = 'raw' (default) - do not use OSA normalisation or dark correct
; OPT = 'dark' - dark correct but no OSA normalisation
; OPT = 'OSA' - OSA normalise but no dark correction
; OPT = 'OSA-dark' - dark correct AND OSA normalize
;
;  ASSUMES identical sampling in both files !!!!!!!!!
;
; OUTPUTS:
;	No explicit outputs.   A new window is created if necessary
;
; COMMON BLOCKS: none
;
; MODIFICATION HISTORY:
; (20-Dec-97 aph) auto-detect version of file format
;                 all of 2-file readin functions put into one program with opt= options
; (30-Aug-13 aph) AXIS standard header added
;-

function xas2, file, PLOT=pl, DEFPATH=defpath, OPT=opt, _extra=e

on_error,2

fltr='*.xas'
if keyword_set(nof) then fltr='*'
file1=PICKFILE2(/Read, FILTER=fltr, /LPATH, DEFPATH=defpath)
file1=strlowcase(file1)
file2=PICKFILE2(/Read, FILTER=fltr, /LPATH, DEFPATH=defpath)
file2=strlowcase(file2)
; print, 'It from ', file1, '.    Io from ', file2

file1=strmid(file1,0,strlen(file1)-4)
file2=strmid(file2,0,strlen(file2)-4)
dat1=text_read1(file1+'.dat')
dat2=text_read1(file2+'.dat')
a1=text_read1(file1+'.xas')
a2=text_read1(file2+'.xas')

ics1=a1(0,*)   ; Energy scales for each file
ics2=a2(0,*)

; AUTO-detect file type (aph - 6-dec-97)       FILE # 1 (sample signal)
dat=dat1
a = a1
test = size(dat)
if test(1) EQ 3 then begin
	nul1=dat(0,1)	; (format A) pre-Nov-97 - 1-section file
	samp1=dat(2,0)
	dw = dat(2,1)/1e3  ; convert dwell in usec to msec
	b1 = a(1,*)/dw  & c1 = a(2,*)/dw
endif else begin		; post Nov-97 formats - multi-section files
	b = a(1,*)  & c = a(2,*)
; (format B)(Nov-20 to 1-Dec-97) - no dark channel info - deduce from file
    if dat(1,0) NE (test(2)-1) then begin  ;(format B)
        nseg = test(2)                 ; number of segments
        xtst = size(a)
        samp = total(dat(2,*))          ; total number of data points
        nul = (xtst(2) - samp)/2       ; number of dark points
        first = 0 & last = nseg - 1
    endif else begin
; (format C) (after 3-Dec-97) - 1st line of dat has dark channel info
		nseg = dat(1,0)               ;(format C)
		samp  = total(dat(2,1:*))
    	nul = dat(0,0)
    	first = 1 & last = nseg
    endelse
; common to both post_Nov-97 formats
    nb = fix(nul)
    for i = first,last do begin    ; correct for different dwell for each section
       dw = dat(3,i)/1e3
       nf = fix(nb + dat(2,i))
       print,'section ',i,' dwell (ms) = ',dw, ' npts = ',nf-nb
       b(nb:nf) = b(nb:nf)/dw
       c(nb:nf) = c(nb:nf)/dw
       nb = nf + 1
    endfor
; correct dark counts for dwell (NB can have error here for format B)
 dw = dat(3,0)/1e3
 ;help,b
 if nul GT 0 then begin           ; to handle 0 dark point data ALS Feb14/98
 	 b(0:nul-1) = b(0:nul-1)/dw
	 b(nul+samp:*)=b(nul+samp:*)/dw
	 c(0:nul-1) = c(0:nul-1)/dw
	 c(nul+samp:*) = c(nul+samp:*)/dw
 endif
; print, '#dark= ',fix(nul),' dwell dark (ms) ',dw,'   #data= ', samp
; load into arrays for file 1
	nul1 = nul  & samp1 = samp
	b1 = b  &  c1 = c
	endelse

; AUTO-detect file type (aph - 6-dec-97)       FILE # 2 (hole signal)
dat=dat2                             ; NB does not assume will use same format
a = a2
test = size(dat)
if test(1) EQ 3 then begin
	nul2=dat(0,1)	; (format A) pre-Nov-97 - 1-section file
	samp2=dat(2,0)
	dw = dat(2,1)/1e3  ; convert dwell in usec to msec
	b2 = a(1,*)/dw  & c2 = a(2,*)/dw
endif else begin		; post Nov-97 formats - multi-section files
	b = a(1,*)  & c = a(2,*)
; (format B)(Nov-20 to 1-Dec-97) - no dark channel info - deduce from file
    if dat(1,0) NE (test(2)-1) then begin  ;(format B)
        nseg = test(2)                 ; number of segments
        xtst = size(a)
        samp = total(dat(2,*))          ; total number of data points
        nul = (xtst(2) - samp)/2       ; number of dark points
        first = 0 & last = nseg - 1
    endif else begin
; (format C) (after 3-Dec-97) - 1st line of dat has dark channel info
		nseg = dat(1,0)               ;(format C)
		samp  = total(dat(2,1:*))
    	nul = dat(0,0)
    	first = 1 & last = nseg
    endelse
; common to both post_Nov-97 formats
    nb = fix(nul)
    for i = first,last do begin    ; correct for different dwell for each section
       dw = dat(3,i)/1e3
       nf = fix(nb + dat(2,i))
;       print,'section ',i,' dwell (ms) = ',dw, ' npts = ',nf-nb
       b(nb:nf) = b(nb:nf)/dw
       c(nb:nf) = c(nb:nf)/dw
       nb = nf + 1
    endfor
; correct dark counts for dwell (NB can have error here for format B)
 dw = dat(3,0)/1e3
 if nul GT 0 then begin
 	b(0:nul-1) = b(0:nul-1)/dw
 	b(nul+samp:*)=b(nul+samp:*)/dw
 	c(0:nul-1) = c(0:nul-1)/dw
 	c(nul+samp:*) = c(nul+samp:*)/dw
 endif
; print, '#dark= ',fix(nul),' dwell dark (ms) ',dw,'   #data= ', samp
; load into arrays for file 2
nul2 = nul  & samp2 = samp
b2 = b  &  c2 = c
endelse

; processing common to all 4 options  -------------
ics1 = ics1(nul1:nul1+samp1-1)
ics2 = ics2(nul2:nul2+samp2-1)
b1d=b1(nul1:nul1+samp1-1)
c1d=c1(nul1:nul1+samp1-1)
b2d=b2(nul2:nul2+samp2-1)
c2d=c2(nul2:nul2+samp2-1)
if nul1 GT 0 then begin
	b1i=b1(0:nul1-1)
	b1f=b1(nul1+samp1:*)
	c1i=c1(0:nul1-1)
	c1f=c1(nul1+samp1:*)
	m1i=median(b1i)
	m1f=median(b1f)
	n1i=median(c1i)
	n1f=median(c1f)
endif
if nul2 GT 0 then begin
	b2i=b2(0:nul2-1)
	b2f=b2(nul2+samp2:*)
	c2i=c2(0:nul2-1)
	c2f=c2(nul2+samp2:*)
	m2i=median(b2i)
	m2f=median(b2f)
	n2i=median(c2i)
	n2f=median(c2f)
endif
samp = samp1
ics = ics1
; ------ interpolate to finest mesh over common scale -------
; (aph - allows for analysis of files with unequal mesh size)
if samp1 NE samp2 OR ics1(0) NE ics2(0) OR ics1(samp1-1) NE ics2(samp2-1)then begin
  print, '        npts         begin          end '
  print, '1    ',samp1, ics1(0),ics1(samp1-1)
  print, '2    ',samp2, ics2(0),ics2(samp2-1)
  print, 'WARNING !!! : Files not on same energy scale - INTERPOL-ated to common scale'
  min1=min(ics1,max=max1)
  min2=min(ics2,max=max2)
  xmin=max(min1,min2)
  xmax=min(max1,max2)
  xt1=ics1(where(ics1 GE xmin AND ics1 LE xmax,count1))
  xt2=ics2(where(ics2 GE xmin AND ics2 LE xmax,count2))
  if (count1 NE 0 AND count2 NE 0) then begin            ; use finer mesh
	if count1 GT count2 then begin
		b2d = interpol(b2d(where(ics2 GE xmin AND ics2 LE xmax)),xt2,xt1)
		c2d = interpol(c2d(where(ics2 GE xmin AND ics2 LE xmax)),xt2,xt1)
		b1d = b1d(where(ics1 GE xmin AND ics1 LE xmax))
		c1d = c1d(where(ics1 GE xmin AND ics1 LE xmax))
		samp = count1
		ics = xt1
	endif else begin
		b1d = interpol(b1d(where(ics1 GE xmin AND ics1 LE xmax)),xt1,xt2)
		c1d = interpol(c1d(where(ics1 GE xmin AND ics1 LE xmax)),xt1,xt2)
		b2d = b2d(where(ics2 GE xmin AND ics2 LE xmax))
		c2d = c2d(where(ics2 GE xmin AND ics2 LE xmax))
		samp = count2
		ics = xt2
	endelse
  endif else begin
  	print, 'ERROR !!! - files do not share any common energies'
  	s = 0
  	return, s
  	endelse
endif

; OPTION-specific processing ------------------
if NOT keyword_set(opt) then OPT = 'raw' ; default
; print, 'OPT= ',OPT
; ************************* no OSA or dark corrections ****************
if OPT EQ 'raw' then begin
	uf=(b1d/b2d)
	fin=-alog(uf)
	t_tag = ' no osa, no dark'
	ytit = '-ln[IN /OUT] (no dark corr.)'
endif
; ********************** OSA normalisation (ONLY) ***************
if OPT EQ 'OSA' then begin
	b1d=b1d/c1d
	b2d=b2d/c2d
	uf=(b1d/b2d)
	fin=-alog(uf)
	t_tag = ' osa, no dark'
	ytit = '-ln[(IN/osa / OUT/osa)] (no dark corr.)'
endif
; ***********************************  dark-only (no OSA norm) ************************
if OPT EQ 'dark' then begin
	bia1=findgen(samp)*(m1f(0)-m1i(0))/samp + m1i(0)
	bia2=findgen(samp)*(m2f(0)-m2i(0))/samp + m2i(0)
	cia1=findgen(samp)*(n1f(0)-n1i(0))/samp + n1i(0)
	cia2=findgen(samp)*(n2f(0)-n2i(0))/samp + n2i(0)
	b1d = b1d-bia1
	b2d = b2d-bia2
	c1d = c1d-cia1
	c2d = c2d-cia2
	uf=(b1d(where(b2d NE 0))/b2d(where(b2d NE 0)))
	fin=-alog(uf)
	t_tag = ' no osa, dark'
	ytit = '-ln[IN/OUT] [dark corrected]'
endif
; ***********************************  OSA-dark ************************
if OPT EQ 'OSA-dark' then begin
	bia1=findgen(samp)*(m1f(0)-m1i(0))/samp + m1i(0)
	bia2=findgen(samp)*(m2f(0)-m2i(0))/samp + m2i(0)
	cia1=findgen(samp)*(n1f(0)-n1i(0))/samp + n1i(0)
	cia2=findgen(samp)*(n2f(0)-n2i(0))/samp + n2i(0)
	b1d = b1d-bia1
	b2d = b2d-bia2
	c1d = c1d-cia1
	c2d = c2d-cia2
	b1d=b1d/c1d
	b2d=b2d/c2d
	uf=(b1d(where(b2d NE 0))/b2d(where(b2d NE 0)))
	fin=-alog(uf)
	t_tag = ' osa & dark'
	ytit = 'ln[(IN/OSA)/(OUT/OSA)] [dark corrected]'
endif

file1=strmid(file1,strlen(file1)-8,strlen(file1))
file2=strmid(file2,strlen(file2)-8,strlen(file2))

s={t:'1d', d:fin, x:ics, xl:'Photon Energy',dn:c1, dl:+file1+'/'+file2+t_tag}

if keyword_set(pl)then  splot,s,ytitle= ytit,_extra=e, charsize=1

return, s
end




