; $Id: dcs_save_scattering.pro,v 1.19 2005/07/15 16:32:09 jrdcopley Exp $
; Written by J.R.D. Copley.

;************************************************************************************************
;###############################################################################
;
; LICENSE:
;  The software in this file is written by an employee of
;  National Institute of Standards and Technology
;  as part of the DAVE software project.
;
;  The DAVE software package is not subject to copyright protection
;  and is in the public domain. It should be considered as an
;  experimental neutron scattering data reduction, visualization, and
;  analysis system. As such, the authors assume no responsibility
;  whatsoever for its use, and make no guarantees, expressed or
;  implied, about its quality, reliability, or any other
;  characteristic. The use of certain trade names or commercial
;  products does not imply any endorsement of a particular product,
;  nor does it imply that the named product is necessarily the best
;  product for the stated purpose. We would appreciate acknowledgment
;  if the DAVE software is used or if the code in this file is
;  included in another product.
;
;###############################################################################
pro dcs_save_scattering_strategy_quit,event
;************************************************************************************************
;
compile_opt strictarr
;
widget_control,event.top,/destroy
return
end


;************************************************************************************************
pro dcs_save_scattering_strategy_cleanup,tlb
;************************************************************************************************
;
compile_opt strictarr
;
widget_control,tlb,get_uvalue=pState
;
ptr_free,pState
;
end


;************************************************************************************************
pro dcs_save_scattering_strategy_handler,event
;************************************************************************************************
;
compile_opt strictarr
;
widget_control,event.top,get_uvalue=pState
nfiles=(*pState).nfiles
tempswidget=(*pState).tempswidget
;
save_params=*(*pState).save_paramsPtr
temperatures=*(*pState).temperaturesPtr; (floats)
;
evid=event.id

if (evid eq (*pState).done) then begin
	dcs_save_scattering_strategy_quit,event
	return
endif
;
if (evid eq (*pState).cancel) then begin
	save_params[3]="Cancel"
	*(*pState).save_paramsPtr=save_params
	widget_control,event.top,set_uvalue=pState
	dcs_save_scattering_strategy_quit,event
	return
endif
;
widget_type=widget_info(evid,/type)
;
; widget_type=1 for a button widget.
if (widget_type eq 1) then begin
	evsel=event.select
	if (evsel eq 0) then begin
		case evid of
			(*pState).compress: save_params[2]="NOCOMPRESS"
			(*pState).deteffchoice: *(*pState).dodeteffPtr=0
			(*pState).zeroenergychoice: *(*pState).dozeroenergyPtr=0
			else:
		endcase
	endif
	if (evsel eq 1) then begin
		case evid of
			(*pState).dodt: save_params[0]="DODT"
			(*pState).dode: save_params[0]="DODE"
			(*pState).sqom: save_params[0]="SQOM"
			(*pState).ssqw: save_params[0]="SSQW"
			(*pState).kiqw: save_params[0]="KIQW"
			(*pState).gdos: save_params[0]="GDOS"
			(*pState).vdos: save_params[0]="VDOS"
			(*pState).sgdos: save_params[0]="SGDOS"
			(*pState).save_dave: save_params[1]="SAVE_DAVE"
			(*pState).octaveasc: save_params[1]="OCTAVEASC"
			(*pState).octavebin: save_params[1]="OCTAVEBIN"
			(*pState).idlformat: save_params[1]="IDLFORMAT"
			(*pState).treatform: save_params[1]="TREATFORM"
			(*pState).specops: begin
				operation="unknown"
				dcs_davePtr_ops,operation=operation,$
					davePtr=(*pState).davePtr,$
					base=event.top,$
					message_str=(*pState).message_str
				*(*pState).operationPtr=operation
			end
			(*pState).compress: save_params[2]="COMPRESS"
			(*pState).deteffchoice: *(*pState).dodeteffPtr=1
			(*pState).zeroenergychoice: *(*pState).dozeroenergyPtr=1
			else:
		endcase
		case evid of
			(*pState).dodt: widget_control,(*pState).tempvalsbase,sensitive=0
			(*pState).dode: widget_control,(*pState).tempvalsbase,sensitive=0
			(*pState).sqom: widget_control,(*pState).tempvalsbase,sensitive=0
			(*pState).ssqw: widget_control,(*pState).tempvalsbase,sensitive=1
			(*pState).kiqw: widget_control,(*pState).tempvalsbase,sensitive=1
			(*pState).gdos: widget_control,(*pState).tempvalsbase,sensitive=1
			(*pState).vdos: widget_control,(*pState).tempvalsbase,sensitive=1
			(*pState).sgdos: widget_control,(*pState).tempvalsbase,sensitive=1
			else:
		endcase
	endif
endif
;
for k=0,nfiles-1 do begin
	if (evid eq (*pState).tempswidget[k]) then begin
		if (event.type eq 0) then widget_control,evid,get_value=temperatura
		temperatura=float(stregex(temperatura,'[0-9.]+',/extract)); string to float
		temperatures[k]=temperatura
		temperatura="T = "+$
			string(stregex(temperatura,'[0-9.]+',/extract),format='(f7.2)')+"K"; float to string
		widget_control,evid,set_value=temperatura
	endif
endfor
;
*(*pState).save_paramsPtr=save_params
*(*pState).temperaturesPtr=temperatures; floats
widget_control,event.top,set_uvalue=pState
;
end


;************************************************************************************************
pro dcs_save_scattering_strategy,dcs_base,$
	save_paramsPtr,$
	temperaturesPtr,$
	dodeteffPtr,$
	dozeroenergyPtr,$
	operationPtr,$
	davePtr=davePtr,$
	message_str=message_str
;************************************************************************************************
;
compile_opt strictarr
;
temperatures=*temperaturesPtr; (floats)
nfiles=n_elements(temperatures)
samswidget=intarr(nfiles)
tempswidget=intarr(nfiles)
temperaturestrings=string(temperatures,format='(f7.2)'); floats to strings
;
tlb=widget_base(title='Save Strategy',/row,group_leader=dcs_base,/modal)
	leftcol=widget_base(tlb,/col)
		inchoice=widget_base(leftcol,/col,/exclusive,/frame)
			dodt=widget_button(inchoice,value="Save as ''d2sigma/dOmegadt''")
			dodE=widget_button(inchoice,value="Save as ''d2sigma/dOmegadEf''")
			sqom=widget_button(inchoice,value="Save as ''S(Q,omega)''")
			ssqw=widget_button(inchoice,value="Save as ''Symmetrized S(Q,omega)''")
			kiqw=widget_button(inchoice,value="Save as ''Chi'"+"'(Q,omega)''")
			gdos=widget_button(inchoice,value="Save as ''g(omega)''")
;			if (not lmgr(/runtime)) then $
				vdos=widget_button(inchoice,value="Save as ''vdos(omega)''") ;else vdos=0l
			sgdos=widget_button(inchoice,value="Save as summed ''g(omega)''")
		outchoice=widget_base(leftcol,/col,/exclusive,/frame)
			save_dave=widget_button(outchoice,value="Save in DAVE format")
			octaveasc=widget_button(outchoice,value="Save in Octave ascii format")
			octavebin=widget_button(outchoice,value="Save in Octave binary format")
;			idlformat=widget_button(outchoice,value="Save in internal IDL format")
;			treatform=widget_button(outchoice,value="Save in TREAT output format")
		compchoice=widget_base(leftcol,/col,/nonexclusive,/frame)
			compress=widget_button(compchoice,value="Compress the output data?")
	midcol=widget_base(tlb,/col)
		upsec=widget_base(midcol,/col,/frame)
			void=widget_label(upsec,value="Check the box below to correct for")
			void=widget_label(upsec,value="energy-dependent detector efficiency")
			deteffbase=widget_base(upsec,/col,/nonexclusive)
				deteffchoice=widget_button(deteffbase,value="Apply the correction")
		midsec=widget_base(midcol,/col,/frame)
			void=widget_label(midsec,value="Check the box below to correct for")
			void=widget_label(midsec,value="zero energy transfer channel")
			zeroenergybase=widget_base(midsec,/col,/nonexclusive)
				zeroenergychoice=widget_button(zeroenergybase,value="Apply the correction")
		lowsec=widget_base(midcol,/col,/frame)
			samtemptit1=widget_label(lowsec,value="Sample temperatures")
			samtemptit2=widget_label(lowsec,value="(needed to compute")
			samtemptit3=widget_label(lowsec,value="detailed balance factors)")
			samtempbase=widget_base(lowsec,/row)
				samvalsbase=widget_base(samtempbase,/col)
					for k=0,nfiles-1 do $
						samswidget[k]=widget_text(samvalsbase,xsize=8,value="Series"+strcompress(k+1,/remove)+":")
				tempvalsbase=widget_base(samtempbase,/col,sensitive=0)
					for k=0,nfiles-1 do $
						tempswidget[k]=widget_text(tempvalsbase,xsize=11,$
							value="T ="+strcompress(temperaturestrings[k])+"K",/editable)
;		specialops=widget_base(midcol,/col,/nonexclusive,/frame)
;			specops=widget_button(specialops,value="Special operations?")
		done=widget_button(midcol,value="Proceed")
		cancel=widget_button(midcol,value="Cancel")
;
dave_position_tlb,tlb
widget_control,tlb,/realize
;
widget_control,sqom,set_button=1
widget_control,save_dave,set_button=1
widget_control,deteffchoice,set_button=1
;
*save_paramsPtr=["SQOM","SAVE_DAVE","NOCOMPRESS","Proceed"]
*dodeteffPtr=1
;
state={$
	dodt:dodt,$
	dode:dode,$
	sqom:sqom,$
	ssqw:ssqw,$
	kiqw:kiqw,$
	gdos:gdos,$
	vdos:vdos,$
	sgdos:sgdos,$
	tempswidget:tempswidget,$
	done:done,$
	cancel:cancel,$
	save_paramsPtr:save_paramsPtr,$
	temperaturesPtr:temperaturesPtr,$; (floats)
	dodeteffPtr:dodeteffPtr,$
	dozeroenergyPtr:dozeroenergyPtr,$
	operationPtr:operationPtr,$
	nfiles:nfiles,$
	save_dave:save_dave,$
	octaveasc:octaveasc,$
	octavebin:octavebin,$
	idlformat:-1,$;idlformat,$
	treatform:-1,$;treatform,$
	compress:compress,$
	deteffchoice:deteffchoice,$
	zeroenergychoice:zeroenergychoice,$
	specops:-1,$;specops,$
	tempvalsbase:tempvalsbase,$
	davePtr:davePtr,$
	message_str:message_str}
;
pState=ptr_new(state)
widget_control,tlb,set_uvalue=pState
;
xmanager,'dcs_save_scattering_strategy',tlb,$
	event_handler='dcs_save_scattering_strategy_handler',$
	cleanup='dcs_save_scattering_strategy_cleanup'
;
end


;************************************************************************************************
pro dcs_save_scattering,dcs_base,samfiles,temperatures,logentry=logentry,message_str=message_str
;************************************************************************************************
;
compile_opt strictarr
;
restore,filename=samfiles[0]
;
save_paramsPtr=ptr_new(strarr(4))
temperaturesPtr=ptr_new(temperatures); (floats)
dodeteffPtr=ptr_new(1)
dozeroenergyPtr=ptr_new(0)
operationPtr=ptr_new("")
;
dcs_save_scattering_strategy,dcs_base,$
	save_paramsPtr,$
	temperaturesPtr,$
	dodeteffPtr,$
	dozeroenergyPtr,$
	operationPtr,$
	davePtr=davePtr,$
	message_str=message_str
;
save_params=*save_paramsPtr
temperatures=*temperaturesPtr
dodeteff=*dodeteffPtr
dozeroenergy=*dozeroenergyPtr
operation=*operationPtr
;
ptr_free,save_paramsPtr
ptr_free,temperaturesPtr
ptr_free,dodeteffPtr
ptr_free,dozeroenergyPtr
ptr_free,operationPtr
;
; Values of kT in meV.
ktemperatures=float(stregex(temperatures,'[0-9.]+',/extract))*!dcs_boltzmann/!dcs_meVcon
;
if (save_params[3] eq "Cancel") then return
;
if (dodeteff) then begin
	message="Performing detector efficiency corrections."
endif else begin
	message="SKIPPING detector efficiency corrections."
endelse
logentry=logentry+string(10B)+message
print,message
;
n_samfiles=n_elements(samfiles)
;
lines=strarr(2*n_samfiles+5)
lines[0:1]=["Please wait while the calculations proceed.",$
						"---------------------------------------------"]
bases=intarr(2*n_samfiles)
;
ngtchan0=(*(*(*davePtr).dataStrPtr).specificPtr).ngtchan0
;
; Uncorrected times of flight and time channel boundaries
times_uncorr=(*(*(*(*davePtr).dataStrPtr).specificPtr).time_propsPtr)[0:ngtchan0]
tstart_uncorr=times_uncorr[0:ngtchan0-1]
tstop_uncorr=times_uncorr[1:ngtchan0]
;
e0=!dcs_hsq2mn/(*(*(*davePtr).dataStrPtr).specificPtr).ch_wl^2
;
tsd0=(*(*(*davePtr).dataStrPtr).specificPtr).ch_wl/!dcs_hom*!dcs_dsd
;
for k=0,n_samfiles-1 do begin
;
; If requested, determine zero energy channel time correction.
	if (dozeroenergy) then begin
		peakfrac=0.5
		summed_int=total((*(*(*davePtr).dataStrPtr).commonStr.histPtr).qty,2)
		dcs_getcentroid,summed_int,peakfrac,tstart_uncorr,tstop_uncorr,tcentroid,error
		tcorr=tsd0-tcentroid
	endif else begin
		tcorr=0.0
	endelse
;
; Corrected times and energies.
	times=times_uncorr+tcorr
	energies=e0-!dcs_hsq2mn*!dcs_dsd^2/!dcs_hom^2/times^2
;
; Boundaries of corrected time and energy channels.
	tstart=times[0:ngtchan0-1]
	tstop=times[1:ngtchan0]
	etstart=energies[0:ngtchan0-1]
	etstop=energies[1:ngtchan0]
;
; Mean times of flight from sample to detector
	tmean=0.5*(tstart+tstop)
;
; Mean energy transfers
	etmean=0.5*(etstart+etstop)
;
; Mean scattering angles and elastic Q values
	angle_props=*(*(*(*davePtr).dataStrPtr).specificPtr).angle_propsPtr
	Q_props=*(*(*(*davePtr).dataStrPtr).specificPtr).Q_propsPtr
	tthmean=angle_props[*,0]
	qelmean=Q_props[*,0]
	ntth=n_elements(tthmean)
;
; Values of 2*cos(2theta) and sin(2theta)
	twocostth=2.0*cos(tthmean*!dtor)
	sintth=sin(tthmean*!dtor)
;
	ef=e0-etmean
	root=sqrt(e0*ef)
;
; Ratios (k_i/k_f)
	kikfratio=root/ef
;
; Inverse of detector efficiency as function of Ef.
	if (dodeteff) then begin
		det_thick=1.0; assumed detector thickness in cm.
		sigma_acoeff=0.4466; (from Jeremy in cm-1)
		sigma_s=0.0008; (from Jeremy in cm-1)
		wl0=(*(*(*davePtr).dataStrPtr).specificPtr).ch_wl
		deteff_inv=1.0/(1.0-exp(-det_thick*(sigma_s+sigma_acoeff*wl0*kikfratio)))
	endif else begin
		deteff_inv=fltarr(ngtchan0)+1
	endelse
;
	unitvectth=fltarr(ntth)+1
;
; Create 2-d array of inelastic Q values
; Next line creates matrices with dimensions [ngtchan0,ntth].
	hsqqsq=dcs_operate((e0+ef)#unitvectth,root#twocostth,"-")
	conversion=2*!pi/sqrt(!dcs_hsq2mn)
	qinel=sqrt(hsqqsq)*conversion
;
; Create 1-d array of (maximum Q^4 - minimum Q^4) values
	tthmin=10.0
	tthmax=140.0
	q4maxmin=(e0+ef-2.0*root*cos(tthmax*!dtor))^2 - (e0+ef-2.0*root*cos(tthmin*!dtor))^2
;
	compressed=save_params[2] eq "COMPRESS"
;
	out_dataPtr=ptr_new({description:save_params,tmean:tmean,etmean:etmean,$
		tthmean:tthmean,qelmean:qelmean,qinel:qinel,$
		qty:fltarr(ngtchan0,ntth),err:fltarr(ngtchan0,ntth)})
;
	filename=samfiles[k]
	lines[k*2+3]="Reading input file "+filename
	dcs_message_display,lines,base,group_leader=dcs_base,alignment=0
	bases[k*2]=base
	free_all_ptrs,davePtr
	restore,filename=filename
;
; Values of exp(-E/2kT)
	dbal=exp(-etmean/(2.0*ktemperatures[k]))
;
; Lines containing matrix multiplication (#) create matrices with
; dimensions [ngtchan0,ntth].
	case save_params[0] of
;	  "DODT": factor=twid_inv#unitvectth
;	  "DODE": factor=(twid_inv*kikfratio^3)#unitvectth
;	  "SQOM": factor=(twid_inv*kikfratio^4)#unitvectth
;	  "SSQW": factor=(twid_inv*kikfratio^4*dbal)#unitvectth
;	  "KIQW": factor=(twid_inv*kikfratio^4*(1-dbal^2))#unitvectth
;	  "GDOS": factor=((twid_inv*kikfratio^4*(1-dbal^2)*etmean)#unitvectth)/hsqqsq
;	  "VDOS": factor=(twid_inv*kikfratio^3*(1-dbal^2)*etmean/q4maxmin)#sintth
;		"SGDOS": factor=((twid_inv*kikfratio^4*(1-dbal^2)*etmean)#unitvectth)/hsqqsq
	  "DODT": factor=deteff_inv#unitvectth
	  "DODE": factor=(deteff_inv*kikfratio^3)#unitvectth
	  "SQOM": factor=(deteff_inv*kikfratio^4)#unitvectth
	  "SSQW": factor=(deteff_inv*kikfratio^4*dbal)#unitvectth
	  "KIQW": factor=(deteff_inv*kikfratio^4*(1-dbal^2))#unitvectth
	  "GDOS": factor=((deteff_inv*kikfratio^4*(1-dbal^2)*etmean)#unitvectth)/hsqqsq
	  "VDOS": factor=(deteff_inv*kikfratio^3*(1-dbal^2)*etmean/q4maxmin)#sintth
	  "SGDOS": factor=((deteff_inv*kikfratio^4*(1-dbal^2)*etmean)#unitvectth)/hsqqsq
	  else: if (!debug) then print,"Error in dcs_save_scattering: STOP"
	endcase
;
	case save_params[0] of
	  "DODT": what="d2sigma/dOmegadt"
	  "DODE": what="d2sigma/dOmegadE"
	  "SQOM": what="S(Q,omega)"
	  "SSQW": what="''Symmetrized'' S(Q,omega)"
	  "KIQW": what="Chi'"+"'(Q,omega)"
	  "GDOS": what="g(omega)"
	  "VDOS": what="vdos(omega)"
	  "SGDOS": what="Summed g(omega)"
  else:
	endcase
;
	(*out_dataPtr).qty=factor * $
		(*(*(*davePtr).dataStrPtr).commonStr.histPtr).qty
	(*out_dataPtr).err=factor * $
		(*(*(*davePtr).dataStrPtr).commonStr.histPtr).err
;
	logentry=logentry+string(10B)+"Sample run number"+strcompress(k+1)+":"
;
	case save_params[1] of
		"SAVE_DAVE": ext=".dave"
		"OCTAVEASC": ext=".ASC"
		"OCTAVEBIN": ext=".BIN"
		"IDLFORMAT": ext=".IDL"
		"TREATFORM": ext=".TRE"
	endcase
	newfilename=strmid(filename,0,strlen(filename)-8)+save_params[0]+ext
	if (compressed) then newfilename=newfilename+".gz"
;
	lines[k*2+4]="Writing output file "+newfilename
	dcs_message_display,lines,base,group_leader=dcs_base,alignment=0
	bases[k*2+1]=base
	logentry=logentry+string(10B)+$
		what+" is to be written to file "+newfilename+"."
;
	if (dozeroenergy) then begin
		message="Zero energy channel correction,"+strcompress(string(tcorr))+" us."
	endif else begin
		message="NO zero energy channel correction."
	endelse
	logentry=logentry+string(10B)+message
;
	bigfloat=exp(80.0)
	indices=where(1-finite((*out_dataPtr).qty),countq)
	if (countq gt 0) then begin
		(*out_dataPtr).qty[indices]=bigfloat
	endif
	indices=where(1-finite((*out_dataPtr).err),counte)
	if (counte gt 0) then begin
		(*out_dataPtr).err[indices]=bigfloat
	endif
	if (countq gt 0 or counte gt 0) then begin
		res=dialog_message(strcompress(countq)+" 'infinite' elements of qty and"+$
			strcompress(counte)+"' infinite' elements of err have been corrected")
	endif
;
	case save_params[1] of
		"SAVE_DAVE": begin
			(*(*(*davePtr).dataStrPtr).commonStr.histPtr).qty=(*out_dataPtr).qty
			(*(*(*davePtr).dataStrPtr).commonStr.histPtr).err=(*out_dataPtr).err
			if (save_params[0] ne "DODT") then begin
				(*(*(*davePtr).dataStrPtr).commonStr.histPtr).x=energies
				(*(*davePtr).dataStrPtr).commonStr.xunits="energy:meV"
				(*(*davePtr).dataStrPtr).commonStr.xlabel="Energy transfer"
				(*(*davePtr).dataStrPtr).commonStr.histlabel=what
			endif else begin
				(*(*(*davePtr).dataStrPtr).commonStr.histPtr).x=times
			endelse
			if (save_params[0] eq "SGDOS") then begin
				(*(*(*davePtr).dataStrPtr).commonStr.histPtr).qty=$
					total((*(*(*davePtr).dataStrPtr).commonStr.histPtr).qty,2)
				(*(*(*davePtr).dataStrPtr).commonStr.histPtr).err=$
					sqrt(total((*(*(*davePtr).dataStrPtr).commonStr.histPtr).err^2,2))
				(*(*(*davePtr).dataStrPtr).commonStr.histPtr).y=0.0
			endif
			if (operation eq "ytoQofE") then begin
				dcs_davePtr_ops,davePtr=davePtr,operation=operation,base=dcs_base
				logentry=logentry+string(10B)+$
				"Transforming Q to Q(E) along k0."
			endif
			if (operation eq "ytoQcompofE") then begin
				dcs_davePtr_ops,davePtr=davePtr,operation=operation,base=dcs_base
				logentry=logentry+string(10B)+$
				"Transforming Q to Q component along k0."
			endif
			(*(*(*davePtr).dataStrPtr).commonStr.treatmentPtr)=strsplit(logentry,string(10B),/extract)
			save,davePtr,filename=newfilename
		end
		"OCTAVEASC": dcs_write_asciiOctavefile,$
			newfilename,datafile=*out_dataPtr,compress=compressed
		"OCTAVEBIN": dcs_write_binaryOctavefile,$
			newfilename,datafile=*out_dataPtr,compress=compressed
		"IDLFORMAT":	save,filename=newfilename,*out_dataPtr,compress=compressed
		"TREATFORM": begin
			openw,lun,newfilename,/get_lun,compress=compressed
			for ich=0,ngtchan0-1 do begin
				printf,lun,etmean[ich],ntth
				for ig=0,ntth-1 do printf,lun,qinel[ich,ig],$
					(*out_dataPtr).qty[ich,ig],(*out_dataPtr).err[ich,ig]
			endfor
			free_lun,lun
		end
	endcase
;
endfor
;
;wait,3
for k=2*n_samfiles-1,0,-1 do widget_control,bases[k],/destroy
;
ptr_free,out_dataPtr
;
free_all_ptrs,davePtr
;
end
