; $Id: dcs_read_binaryoctavefile.pro,v 1.22 2005/02/14 17:50:33 azuah 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_read_binaryOctavefile_done_reading,unit,ivar,ivaru,variables_read,unknowns,$
	davePtr=davePtr,fname=fname
;************************************************************************************************	;
;
compile_opt strictarr
;
if (ivar ge 0) then variables_read=variables_read[0:ivar] else variables_read=["None"]
if (ivaru ge 0) then unknowns=unknowns[0:ivaru] else begin
	unknowns=["None"]
endelse
;
free_lun,unit
;
; Populate (*histPtr).x and (*histPtr).y.
ntchan=(*(*(*davePtr).dataStrPtr).specificPtr).ntchan
nhedet=(*(*(*davePtr).dataStrPtr).specificPtr).nhedet
(*(*(*davePtr).dataStrPtr).commonStr.histPtr).x=indgen(ntchan)+1
(*(*(*davePtr).dataStrPtr).commonStr.histPtr).y=indgen(nhedet)+1
;
; With respect to the variables ncycles, duration_csum, he3_csum and high_csum, there are
;	three possible situations:
; 1.	If the raw data file was created before 10/31/03, ncycles is unknown. The variable ngood
;			is determined from temp_control except that in case B (when temperatures were not read)
;			it is set equal to 1. Then ncycles is set equal to ngood.
;			The csum variables are then set to average values determined from ngood.
;	2.	If the raw data file was created after 10/31/03 and before 5/4/04 it does not contain
;			ncycles nor the csum variables but it does contain runinfo so ncycles is set to runinfo[0].
; 		The number of good cycles, ngood, is determined from temp_control except that in case B
;			it is unknown so it is set equal to ncycles.
;			The csum variables are set to average values determined from ngood.
;	3.	If the raw data file was created after 5/4/04 it contains ncycles, duration_csum,
;			he3_csum and high_csum.
;
caseB=(*(*(*davePtr).dataStrPtr).specificPtr).temp_control[0] eq 300.0 and $
	(*(*(*davePtr).dataStrPtr).specificPtr).temp_sample[0] eq 300.0 and $
	(*(*(*davePtr).dataStrPtr).specificPtr).temp_setpoint eq 300.0
;
situation=0
if (fname lt '20031103_003') then situation=1  															;<----mod 9/9/04 JRDC
if (fname ge '20031103_003' and fname lt '20040504_001') then situation=2   ;<----mod 9/9/04 JRDC
if (fname ge '20040504_001') then situation=3
;
; For case B (when temperatures are not read, i.e. when rectemp=0), temp_control and temp_sample
; are scalars but duration_csum, he3_csum and high_csum all contain ngood (or 18*ngood)
;	useful numbers.
;
; For cases A, C and D, temp_control, temp_sample, duration_csum, he3_csum and high_csum
; all contain ngood (or 18*ngood) useful numbers.
;
; The following explanation of cases A, C and D is from dcs_runsinfo.
;
;
; temp_sample and temp_control are initialized in dcs_definedavestructure as vectors of length nvals
; with value -1.0, then in read_binaryoctavefile they become vectors of length nvals, value -1.1
; They are then modified according to the content of the raw data file.	In the comments that follow
;	#cycles is the number of accepted (good) cycles and the descriptions refer to the content of the
;	raw data file.
;
; Case A:	If rectemp=1 but the temperature cannot be read temp_sample and temp_control (in the raw data
;					file) are vectors of length #cycles, value -1.0 (thus temp_sample and temp_control in this
;					procedure, if there were 3 good cycles, would be [-1.0,-1.0,-1.0,-1.1,-1.1,-1.1,-1.1,-1.1,-1.1 ...])
; Case B:	If rectemp=0, temp_sample,temp_control and temp_setpoint (in the raw data file) are scalars,
;					value "300.0", i.e. temp_sample and temp_control in this procedure, would be [300.0,-1.1,-1.1,-1.1,
;					-1.1,-1.1,-1.1,-1.1,-1.1,-1.1,-1.1 ...])
; Case C:	Normally temp_sample and temp_control (in the raw data file) are vectors of length #cycles
;					with positive numbers.
;	Case D:	If the temperature controller were set to <none> temp_sample and temp_control (in the raw data
;					file) would be vectors of length #cycles, value 0.0 (thus temp_sample and temp_control in this
;					procedure, if there were 3 good cycles, would be [0.0,0.0,0.0,-1.1,-1.1,-1.1,-1.1,-1.1,-1.1 ...])
;
;
; Test for creation before 10/31/03.
if (situation eq 1) then begin
	if (caseB) then begin
		ngood=1
	endif else begin
		res=where((*(*(*davePtr).dataStrPtr).specificPtr).temp_control ge 0.0,ngood)
		res=where((*(*(*davePtr).dataStrPtr).specificPtr).temp_control ge -1.05,ngood)   ; <---- this line added 9/3/04
	endelse
	ncycles=ngood
endif
;
; Test for creation after 10/31/03 and before 5/4/04.
if (situation eq 2) then begin
	ncycles=(*(*(*davePtr).dataStrPtr).specificPtr).runinfo[0]
	if (caseB) then begin
		ngood=ncycles
	endif else begin
		res=where((*(*(*davePtr).dataStrPtr).specificPtr).temp_control ge 0.0,ngood)
		res=where((*(*(*davePtr).dataStrPtr).specificPtr).temp_control ge -1.05,ngood)   ; <---- this line added 9/3/04
	endelse
endif
;
; Test for creation before 5/4/04.
if (situation eq 1 or situation eq 2) then begin
	(*(*(*davePtr).dataStrPtr).specificPtr).ncycles=ncycles
	ngood=float(ngood)
	(*(*(*davePtr).dataStrPtr).specificPtr).duration_csum=-1.0
	(*(*(*davePtr).dataStrPtr).specificPtr).he3_csum=-1.0
	(*(*(*davePtr).dataStrPtr).specificPtr).high_csum=-1.0
	(*(*(*davePtr).dataStrPtr).specificPtr).duration_csum[0:ngood-1]=$
		((*(*(*davePtr).dataStrPtr).specificPtr).duration)/ngood
	(*(*(*davePtr).dataStrPtr).specificPtr).he3_csum[0:ngood-1]=$
		total((*(*(*davePtr).dataStrPtr).specificPtr).detsum[0:912])/ngood
	(*(*(*davePtr).dataStrPtr).specificPtr).high_csum[*,0:ngood-1]=$
		((*(*(*davePtr).dataStrPtr).specificPtr).detsum[913:930]/ngood)#(fltarr(ngood)+1)
endif
;
;
;;;; Temporary code to construct artifically varying quantities for test purposes.
;;;doit=0
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;if (doit) then begin
;;;	mult=1+0.07*(-1)^(indgen(ncycles))
;;;	mult[0]=1
;;;	(*(*(*davePtr).dataStrPtr).specificPtr).duration_csum=(*(*(*davePtr).dataStrPtr).specificPtr).duration_csum*mult
;;;	mult=1-0.05*(-1)^(indgen(ncycles))
;;;	mult[0]=1
;;;	(*(*(*davePtr).dataStrPtr).specificPtr).he3_csum=(*(*(*davePtr).dataStrPtr).specificPtr).he3_csum*mult
;;;	mult=1-0.025*(-1)^(indgen(ncycles))
;;;	mult[0]=1
;;;	mult=(fltarr(18)+1)#mult
;;;	(*(*(*davePtr).dataStrPtr).specificPtr).high_csum=(*(*(*davePtr).dataStrPtr).specificPtr).high_csum*mult
;;;endif
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;
;
; Initialize normalization_sum and fbnormalization_sum.
(*(*(*davePtr).dataStrPtr).specificPtr).normalization_sum=fltarr(nhedet)
(*(*(*davePtr).dataStrPtr).specificPtr).fbnormalization_sum=fltarr(nhedet)
;
; Define time channel boundaries in units of time sample-detector.
twid=(*(*(*davePtr).dataStrPtr).specificPtr).tchanlook[0,*]*0.05
tsdmin=(*(*(*davePtr).dataStrPtr).specificPtr).tsdmin
times=tsdmin+[0,reform(total(twid,/cumulative))]
free_ptr,(*(*(*davePtr).dataStrPtr).specificPtr).time_propsPtr
(*(*(*davePtr).dataStrPtr).specificPtr).time_propsPtr=ptr_new(times)
;
; Define time channel boundaries in units of energy transfer.
conversion=!dcs_hsq2mn/(!dcs_hom/!dcs_dsd)^2
wl0=(*(*(*davePtr).dataStrPtr).specificPtr).ch_wl
if (wl0 eq 0.0) then begin
	energies=-conversion/times^2
endif else begin
	energies=!dcs_hsq2mn/wl0^2-conversion/times^2
endelse
free_ptr,(*(*(*davePtr).dataStrPtr).specificPtr).energy_propsPtr
(*(*(*davePtr).dataStrPtr).specificPtr).energy_propsPtr=ptr_new(energies)
;
; Define detector angles (as a vector).
dcs_getangles,angles,nhedet=nhedet
free_ptr,(*(*(*davePtr).dataStrPtr).specificPtr).angle_propsPtr
(*(*(*davePtr).dataStrPtr).specificPtr).angle_propsPtr=ptr_new(angles)
;
; Define elastic Q values (as a vector).
free_ptr,(*(*(*davePtr).dataStrPtr).specificPtr).Q_propsPtr
Qs=4.0*!dpi*sin(0.5*angles*!dtor)/(*(*(*davePtr).dataStrPtr).specificPtr).ch_wl
(*(*(*davePtr).dataStrPtr).specificPtr).Q_propsPtr=ptr_new(Qs)
;
return
end


;************************************************************************************************
pro dcs_read_binaryOctavefile,filepath,nofirstline,omithisto,variables_read,unknowns,error,$
	davePtr=davePtr,histodata_floating=histodata_floating
;************************************************************************************************
;
compile_opt strictarr
;
; This is a procedure that reads a binary Octave file.
; If omithisto=0, nothing is omitted, everything is read.
; If omithisto=1, histohigh is read but histodata is omitted.
; If omithisto=2, histohigh and histodata are both omitted.
; If omithisto=3, histohigh and histodata are both omitted, and there is no call to
; 	dcs_read_binaryOctavefile_done_reading.
;
on_ioerror, no_file
;
fname=file_basename(filepath)
;
if (n_elements(histodata_floating) eq 0) then histodata_floating=0
;
; Determine whether or not file is compressed.
ending=strmid(filepath,2,3,/reverse_offset)
compressed=(ending eq ".gz")
;
; Read header to determine byte ordering ("endianness") of input file.
openr,unit,filepath,compress=compressed,/get_lun
header=bytarr(11)
readu,unit,header
free_lun,unit
;
; When opening the file take into account the endianness of the computer that is **reading** the file.
; Goran Gasparovic pointed out that my original coding (commented out) was incorrect.
case string(header[9]) of
;	"B": openr,unit,filepath,/swap_endian,compress=compressed,/get_lun
;	"L": openr,unit,filepath,compress=compressed,/get_lun
	"B": openr,unit,filepath,/swap_if_little_endian,compress=compressed,/get_lun
	"L": openr,unit,filepath,/swap_if_big_endian,compress=compressed,/get_lun
	else: begin
		error="Cannot determine byte ordering of file named "+filepath
		return
	end
endcase
;
readu,unit,header
;
; Define the types of some variables
nlgth=0l & dldocgf=bytarr(5) & dt=0b & padding=0b & scalar=0.0d
n1=0l & n2=0l & nstrgs=0l & nchars=0l
;
variables_read=strarr(100)
unknowns=strarr(100)
error=""
ivar=-1
ivaru=-1
;
ptr_free,(*(*davePtr).dataStrPtr).commonStr.histPtr
ntchan=(*(*(*davePtr).dataStrPtr).specificPtr).ntchan
nhedet=(*(*(*davePtr).dataStrPtr).specificPtr).nhedet
(*(*davePtr).dataStrPtr).commonStr.histPtr = $
	ptr_new({qty:fltarr(ntchan,nhedet),err:fltarr(ntchan,nhedet),$
	x:fltarr(ntchan),y:fltarr(nhedet)})
;
ptr_free,(*(*davePtr).dataStrPtr).commonStr.treatmentPtr
(*(*davePtr).dataStrPtr).commonStr.treatmentPtr = ptr_new("")
;
nvread=0
;
while (not eof(unit)) do begin
	readu,unit,nlgth; read name_length
	name=bytarr(nlgth) & readu,unit,name & name=string(name); read name
	if ((omithisto eq 2 and name eq "histohigh") or $
		(omithisto eq 1 and name eq "histodata")) then begin
		dcs_read_binaryOctavefile_done_reading,unit,ivar,ivaru,$
			variables_read,unknowns,davePtr=davePtr,fname=fname
		return
	endif
	if (omithisto eq 3 and name eq "histohigh") then begin
		free_lun,unit
		return
	endif
	;
	readu,unit,dldocgf; read 5 bytes of doc_length, doc, and global flag
	readu,unit,dt; read data type
	case dt of
		1: begin; read scalar
			readu,unit,padding
			readu,unit,scalar
				case name of
					"badbrd": (*(*(*davePtr).dataStrPtr).specificPtr).badbrd=float(scalar)
					"baddet": (*(*(*davePtr).dataStrPtr).specificPtr).baddet=float(scalar)
					"ch_ms": (*(*(*davePtr).dataStrPtr).specificPtr).ch_ms=float(scalar)
					"ch_res": (*(*(*davePtr).dataStrPtr).specificPtr).ch_res=fix(scalar)
					"ch_srdenom": (*(*(*davePtr).dataStrPtr).specificPtr).ch_srdenom=fix(scalar)
					"ch_srmode": (*(*(*davePtr).dataStrPtr).specificPtr).ch_srmode=fix(scalar)
					"ch_wl": (*(*(*davePtr).dataStrPtr).specificPtr).ch_wl=float(scalar)
					"coll_amp": (*(*(*davePtr).dataStrPtr).specificPtr).coll_amp=float(scalar)
					"coll_mean": (*(*(*davePtr).dataStrPtr).specificPtr).coll_mean=float(scalar)
					"coll_osc": (*(*(*davePtr).dataStrPtr).specificPtr).coll_osc=fix(scalar)
					"datamax": (*(*(*davePtr).dataStrPtr).specificPtr).datamax=float(scalar)
					"det_dis": (*(*(*davePtr).dataStrPtr).specificPtr).det_dis=float(scalar)
					"duration": (*(*(*davePtr).dataStrPtr).specificPtr).duration=float(scalar)
					"highmax": (*(*(*davePtr).dataStrPtr).specificPtr).highmax=float(scalar)
					"ncycles": (*(*(*davePtr).dataStrPtr).specificPtr).ncycles=long(scalar)
					"nframes": (*(*(*davePtr).dataStrPtr).specificPtr).nframes=long(scalar)
					"shutter_stat": (*(*(*davePtr).dataStrPtr).specificPtr).shutter_stat=fix(scalar); added 10/31/03
					"temp_setpoint": (*(*(*davePtr).dataStrPtr).specificPtr).temp_setpoint=float(scalar)
					"tsdmin": (*(*(*davePtr).dataStrPtr).specificPtr).tsdmin=float(scalar)
					else: begin & ivaru=ivaru+1 & unknowns[ivaru]=name & end
				endcase
				ivar=ivar+1
				variables_read[ivar]=name
		end
		2: begin; read matrix
			readu,unit,n1,n2
			nhi=n1>n2
			nlo=n1<n2
			if (nlo eq 1) then data=dblarr(nhi) else data=dblarr(n1,n2)
			readu,unit,padding
			readu,unit,data
				case name of
					"ch_bid" : (*(*(*davePtr).dataStrPtr).specificPtr).ch_bid=float(data)
					"ch_delay": (*(*(*davePtr).dataStrPtr).specificPtr).ch_delay=float(data)
					"ch_dis": (*(*(*davePtr).dataStrPtr).specificPtr).ch_dis=float(data)
					"ch_input": (*(*(*davePtr).dataStrPtr).specificPtr).ch_input=float(data)
					"ch_phase": (*(*(*davePtr).dataStrPtr).specificPtr).ch_phase=float(data)
					"ch_slots": (*(*(*davePtr).dataStrPtr).specificPtr).ch_slots=float(data)
					"detsum": (*(*(*davePtr).dataStrPtr).specificPtr).detsum=ulong(data)
					"duration_csum": begin
						(*(*(*davePtr).dataStrPtr).specificPtr).duration_csum=-1.1
						(*(*(*davePtr).dataStrPtr).specificPtr).duration_csum=float(data)
					end
					"fc_dis": (*(*(*davePtr).dataStrPtr).specificPtr).fc_dis=float(data)
					"grandsum": (*(*(*davePtr).dataStrPtr).specificPtr).grandsum=ulong(data)
					"he3_csum": begin
						(*(*(*davePtr).dataStrPtr).specificPtr).he3_csum=-1.1
						(*(*(*davePtr).dataStrPtr).specificPtr).he3_csum=float(data)
					end
					"high_csum": begin
						(*(*(*davePtr).dataStrPtr).specificPtr).high_csum=-1.1
						(*(*(*davePtr).dataStrPtr).specificPtr).high_csum=float(data)
					end
					"histodata": begin
						if (histodata_floating) then begin
							(*(*(*davePtr).dataStrPtr).commonStr.histPtr).qty=float(data)
							(*(*(*davePtr).dataStrPtr).commonStr.histPtr).err=sqrt(float(data))
						endif else begin
							(*(*(*davePtr).dataStrPtr).commonStr.histPtr).qty=ulong(data) ; CHANGED 4/15/02
							(*(*(*davePtr).dataStrPtr).commonStr.histPtr).err=sqrt(ulong(data)) ; CHANGED 4/15/02
						endelse
					end
					"histohigh": (*(*(*davePtr).dataStrPtr).specificPtr).histohigh=ulong(data) ; CHANGED 4/15/02
					"motor_pos": (*(*(*davePtr).dataStrPtr).specificPtr).motor_pos=float(data)
					"repeats": (*(*(*davePtr).dataStrPtr).specificPtr).repeats=ulong(data)
					"resets": (*(*(*davePtr).dataStrPtr).specificPtr).resets=ulong(data)
					"runinfo": (*(*(*davePtr).dataStrPtr).specificPtr).runinfo=ulong(data); added 10/31/03
					"tchanlook": (*(*(*davePtr).dataStrPtr).specificPtr).tchanlook=ulong(data)
					"temp_control": begin
						(*(*(*davePtr).dataStrPtr).specificPtr).temp_control=-1.1
						(*(*(*davePtr).dataStrPtr).specificPtr).temp_control=float(data)
					end
					"temp_sample": begin
						(*(*(*davePtr).dataStrPtr).specificPtr).temp_sample=-1.1
						(*(*(*davePtr).dataStrPtr).specificPtr).temp_sample=float(data)
					end
					"timsum": (*(*(*davePtr).dataStrPtr).specificPtr).timsum=ulong(data)
					"totals": (*(*(*davePtr).dataStrPtr).specificPtr).totals=ulong(data)
					else: begin & ivaru=ivaru+1 & unknowns[ivaru]=name & end
				endcase
				ivar=ivar+1
				variables_read[ivar]=name
		end
		5: begin; read string
			readu,unit,nchars; read number of characters in string element
			if (nchars gt 0) then begin
				str=bytarr(nchars) & readu,unit,str; read string element
				str=string(str)
				case name of
					"command": (*(*(*davePtr).dataStrPtr).specificPtr).command=str
					"comments": (*(*(*davePtr).dataStrPtr).specificPtr).comments=str
					"sample_desc": (*(*(*davePtr).dataStrPtr).specificPtr).sample_desc=str
					"start_date": (*(*(*davePtr).dataStrPtr).specificPtr).start_date=str
					"startchoice": (*(*(*davePtr).dataStrPtr).specificPtr).startchoice=str
					"stop_date": (*(*(*davePtr).dataStrPtr).specificPtr).stop_date=str
					"user": (*(*(*davePtr).dataStrPtr).specificPtr).user=str; added 10/31/03
					else: begin & ivaru=ivaru+1 & unknowns[ivaru]=name & end
				endcase
				ivar=ivar+1
				variables_read[ivar]=name
			endif
		end
		7: begin; read string array
			readu,unit,nstrgs; read number of strings in array
			strarray=strarr(nstrgs)
			for j=1,nstrgs do begin
				readu,unit,nchars; read number of characters in string element
				str=bytarr(nchars) & readu,unit,str; read string element
				strarray[j-1]=strjoin(string(str)); store, joined, in strarray
			endfor
			case name of
				"highsource": (*(*(*davePtr).dataStrPtr).specificPtr).highsource=strarray
				else: begin & ivaru=ivaru+1 & unknowns[ivaru]=name & end
			endcase
			ivar=ivar+1
			variables_read[ivar]=name
		end
		else: error="N.B. An unsupported variable type has been encountered: "+strcompress(dt)
	endcase
;
	nvread=nvread+1
;
endwhile
;
;print,nvread,' parameters in ',fname
dcs_read_binaryOctavefile_done_reading,unit,ivar,ivaru,$
	variables_read,unknowns,davePtr=davePtr,fname=fname
return
;
no_file: error="There is no file named "+filepath
end