;$Id: odave__define.pro,v 1.20 2008/01/25 16:08:41 azuah Exp $
;###############################################################################
; 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 of if the code in this file is
;  included in another product.
;
;###############################################################################
;*******************************************************************************
;
; oDave
;
; 	Uses:
;
;		Pointer_recursive_cleanup
;
; I/O
;
;	1 = self -> read_dave()
;	1 = self -> write_dave()
;
;	1 = self -> read_spe()
;	1 = self -> write_spe()
;
;	1 = self -> read_dcs()
;
;	1 = self -> read_asd()
;
;	1 = self -> read_txt()
;	1 = self -> write_txt()
;
;	1 = self -> read_self()
;	1 = self -> write_self()
;
;	self -> print
;	self -> plot
;	self -> view
;
; Processing
;
;	o = self -> rebin()
;	o = self -> cut()
;	self -> transpose
;
; Arithmetic
;
;	o = self -> add()
;	o = self -> sub()
;	o = self -> mul()
;	o = self -> div()
;	[qty,err] = self -> sum()
;	[qty,err] = self -> integrate_2d()
;
; Basic 2
;
;	1 = self -> compare()
;
; Basic 1
;
;	1 = self -> valid()
;	1 = self -> set()
;	1 = self -> get()
;	1 = self -> get_range()
;	self -> empty
;
; Internal
;
;	1 = self -> _self2dave
;
; Lifecycle
;
;	o = obj_new('oDave')
;	o = obj_new('oDave', ...=... )
;	obj_destroy, o
;
;*******************************************************************************

@lib_file

;*******************************************************************************
;23456789 123456789 123456789 123465798 123456789 123456789 123465789 123456789
;       10        20        30        40        50        60        70        80
;*******************************************************************************
;
; I/O Methods
;
;*******************************************************************************
;
; .DAVE files contain 'davePtr'.
;


;*******************************************************************************
function oDave::read_dave, file = file

	on_ioerror, io_error

	if n_elements( file ) ne 1 then begin
		file = dialog_pickfile( /read $
			, title = 'Please choose a DAVE file to read.' $
			, filter = [ '*.dave' ] $
			)
	endif

	if not ( file_test( file ) )[0] then begin
		message, 'Unable to comply.', /continue
		message, 'The file does not exist.', /continue
		return, 0
	endif

	restore, filename = file[0]

	if not self -> _dave2self( davePtr = davePtr ) then begin
		message, 'Unable to convert DAVE to oDAVE.', /continue
		return, 0
	endif

	Pointer_recursive_cleanup, davePtr

	return, 1

	io_error:
		message, 'An I/O error occurred.', /continue
		return, 0
end
;*******************************************************************************
function oDave::write_dave, file = file

	on_ioerror, io_error
	if n_elements( file ) ne 1 then begin
		message, 'Using default file name: oDave.dave', /continue
		file = 'oDave.dave'
	endif

	if not self -> _self2dave( davePtr = davePtr ) then begin
		message, '_self2dave failed.', /continue
		return, 0
	endif

	save, davePtr, filename=file, /compress

	Pointer_recursive_cleanup, davePtr

	return, 1

	io_error:
		message, 'An I/O error occurred while writing DAVE.', /continue
		return, 0

end
;*******************************************************************************
function oDave::read_spe, file = file

	oDave = file_read_spe_one( file = file )
	if obj_valid( oDave ) then begin

		return_value = self -> set( oDave )
		obj_destroy, oDave
		return, return_value

	endif else begin
		message, 'Nonvalid object', /continue
		return, 0
	endelse

end
;*******************************************************************************
function oDave::write_spe, file = file, classic = classic

	on_ioerror, io_error

	if n_elements( file ) ne 1 then begin
		file = dialog_pickfile( /write $
			, title = 'Please choose a SPE file to write.' $
			, filter = [ '*.spe' ] $
			)
	endif

	s = size( *self.qty )
	nx = s[1]
	ny = s[2]

	openw, lun, file, /get_lun
	if not keyword_set( classic ) then begin
		printf,lun, ' Title '
	endif
	printf,lun, nx, ny
	printf,lun, self.xlabel
	printf,lun, (*self.x)[*, 0]
	printf,lun, self.ylabel
	printf,lun, (*self.y)[0, *]

	for i=0,nx - 1 do begin

		printf,lun, self.histlabel
		printf,lun, reform( (*self.qty)[i,*] ), format='(4(e15.5))'
		printf,lun,'  Error'
		printf,lun, reform( (*self.err)[i,*] ), format='(4(e15.5))'

	endfor
	free_lun,lun

	return, 1

io_error:
	message, 'An I/O Error Occurred.', /continue
	return, 0

end


;*******************************************************************************
;
; Read a DCS data file in abs(angles) vs. energy.
;
function oDave::read_dcs_abs, file = files $
	; The desired angle range.
	, angle_range = angle_range

	if n_elements( files ) eq 0 then begin
		files = dialog_pickfile( /read, /multiple_files $
			, title = 'Please select a DCS file to read.' $
			, filter = ['*.dcs', '*.dcs.gz'] $
			)
	endif

	if total( file_test(files, /read) ) ne n_elements( files ) then begin
		message, 'The DCS file, ' + $
			files[0] + ' does not exists.', /continue
		return, 0
	endif

	; Number of useful time channels is fixed at this number.
	n_time = 1000

	; Read files.
	trifenestra_dcs_readfiles, files, d
	trifenestra_dcs_det, ang

	; Find the temperature.
	Tidx = where( d.temp_sample gt 0, nT )
	if (nT ge 1) then begin
		self.T = total( d.temp_sample[ Tidx ] ) / nT
	endif

	; Sort the angle in the ascending order.
	sort_key_ang = sort( abs( ang ) )
	if n_elements( angle_range ) eq 2 then begin
		idx = where( abs(ang[ sort_key_ang ]) le angle_range[1] $
			and abs(ang[ sort_key_ang ]) ge angle_range[0] )
		key_ang = sort_key_ang[ idx ]
	endif else begin
		key_ang = sort_key_ang
	endelse
	ang = abs( ang[ key_ang ] )
	n_ang = n_elements( ang )

	; Calculate time.
	DT = d.tchanlook[0, 0:n_time - 1 ]
		; Delta Time * 20.
	time2 = d.tsdmin + total( DT, /cumulative) * .05
	time1 =	time2 - DT * .05

	; Convert time to energy.
	Ei = 81.8 / d.ch_wl^2
	E1 = 81.8 / (3.956 * time1 / d.det_dis)^2
	E2 = 81.8 / (3.956 * time2 / d.det_dis)^2
	Ef = (E2 + E1 ) / 2d
	DE = E1 - E2
		; Delta Energy
	hw	=	Ei - Ef
	n_hw	=	n_time ; + 1

	; Prepare empty vectors.
	oneset = 1 + fltarr(n_ang)
	onepix = 1 + fltarr(n_hw)


	; Obtain monitor normalized scattering intensity.
	for dcs_k = 0, 16 do begin
		if strmid( d.highsource[dcs_k], 0, 3 ) eq 'FC1' then begin
			dcs_ibm = dcs_k
		endif
	endfor
	monitor = total( d.histohigh[0:999, dcs_ibm] ); / d.duration
	int	= double( transpose(d.histodata[0:(n_time - 1),*]) )

	; Phase Volume
	DP = oneset # (monitor * DE)

	; Arrange int in the angle ascending order.
	int = int[ key_ang, * ]

	; Merge the data points that have the same absolute angle.
	qty = int
	err = int
	vol = int

	angle = [ -1d ]
	n_abs_angles = 0
	beginning = 0
	for i = 1, n_ang - 1 do begin

		if i eq n_ang - 1 or ang[i] ne ang[i-1] then begin

			ending = i

			; Merge from beginning to ending.
			angle = [ angle, ang[beginning] ]
			int_vec = total( int[ beginning:ending, *], 1 )
			err_vec = sqrt( int_vec )
			DP_vec = DP[ beginning, * ] * ( ending - beginning + 1 )

			qty[ n_abs_angles, * ] = int_vec / DP_vec
			err[ n_abs_angles, * ] = err_vec / DP_vec
			vol[ n_abs_angles, * ] = DP_vec
			n_abs_angles = n_abs_angles + 1

			beginning = ending + 1

		endif

	endfor

 	return, self -> set( qty=qty[0:n_abs_angles - 1, *] $
 		, err=err[0:n_abs_angles - 1, *] $
 		, vol=vol[0:n_abs_angles - 1, *] $
 		, xdata=angle[1:n_abs_angles] $
 		, ydata=hw $
 		, xlabel='Angle' $
 		, ylabel='hw' $
 		, xunits='degree' $
 		, yunits='meV' $
 		, histlabel='I' $
 		, histunits='arb. units' $
 		, Ei = Ei $
 		)


end
;*******************************************************************************
;
; Read a DCS data file.
;
function oDave::read_dcs_raw $
	; INPUT
	, file = files $
	; OPTION
	, normalize=normalize $
	, old=old $
	; OUTPUT
	, Ei = Ei $
	, T = T $
	, monitor = monitor_sum $
	, duration = duration $
	, DT = DT

	n_time = 1000

	if n_elements( files ) eq 0 then begin
		files = dave_dialog_pickfile( /read, /multiple_files $
			, title = 'Please select DCS files to read.' $
			, filter = ['*.dcs', '*.dcs.gz'] $
			)
	endif

	if total( file_test(files, /read) ) ne n_elements( files ) then begin
		message, 'The DCS file, ' + $
			files[0] + ' does not exists.', /continue
		return, 0
	endif

	; Many system variables are defined here.
	davesysvars

	; !dave_defaults is defined in dave.pro.
	defsysv, '!dave_defaults', exists=exists
	if not exists then begin
		cd, current=directory
		pdd=ptr_new({datDir:directory,workDir:directory,PDFReader:''})
		defsysv,'!dave_defaults',pdd
	endif

	; DCS system variables are defined here.
	dcs_sysvars

	if files[0] eq '' then return, 0

	dave_makeDavePtr,davePtr=davePtr
	dcs_definedavestructure,davePtr=davePtr

	nofirstline=0
	omithisto=0
	printinfo=0

	o = obj_new('oDave')
	p = obj_new('oDave')

	monitor_sum = 0d
	monitor_old_sum = 0d
	duration = 0d
	Tvec = [0]
	for i = 0, n_elements( files ) - 1 do begin

		dcs_read_binaryOctavefile $
			,files[i], nofirstline, omithisto, variables_read $
			,unknowns,error,davePtr=davePtr

		if (strlen(error) gt 0) then begin
			result=dialog_message(["Error reading binary file.","Exiting ..."])
			obj_destroy, [ o, p ]
			return, 0
		endif

 		; DCS specific pointer.
		sPtr = (*(*davePtr).dataStrPtr).specificPtr
		; Obtain monitor normalized scattering intensity.
		for dcs_k = 0, 16 do begin
			if strmid( (*sPtr).highsource[dcs_k], 0, 3 ) eq 'FC1' then begin
				dcs_ibm = dcs_k
			endif
		endfor
		monitor_old = (*sPtr).detsum[916]
		monitor = total( (*sPtr).histohigh[0:999, dcs_ibm] ); / d.duration
		monitor_sum = monitor_sum + monitor
		monitor_old_sum = monitor_old_sum + monitor_old

		duration = duration + (*sPtr).duration

		; Find the temperature.
		Tidx = where( (*sPtr).temp_sample gt 0, nT )
		if (nT ge 1) then begin
			Tvec = [ Tvec, (*sPtr).temp_sample[ Tidx ] ]
		endif

		void = p -> _dave2self( davePtr = davePtr )
		o -> add, p

	endfor

	T = total( Tvec ) / double( n_elements( Tvec ) - 1 )

	o -> get, qty=qty
	obj_destroy, [o, p]

	dcs_getangles, angles, nhedet=913
	; angles = abs( angles )
	; merge_2d_counts, angles, transpose( qty ), ang, int, v = v
	int = transpose( qty )
	ang = angles
	v = dblarr( n_elements( ang ) ) + 1d


	; Calculate time.
	int = (transpose( int ))[0:1000-1, *]
	DT = transpose( (*sPtr).tchanlook[0, 0:1000 - 1 ] )
	time_d = total(  DT  , /cumulative )
	time = (*sPtr).tsdmin + [ time_d * 0.05 ]

	vx = dblarr( n_elements( time ) ) + 1d
	vy = dblarr( n_elements( ang ) ) + 1d

	; Convert time to energy.
	Ei = 81.81 / (*sPtr).ch_wl^2
	Ef = 81.81 / (3.956 * time / (*sPtr).det_dis)^2
	hw	=	Ei - Ef

	print, 'old/new = ', monitor_old_sum / monitor_sum
	if keyword_set( old ) then monitor_sum = monitor_old_sum

	if keyword_set( normalize ) then begin

		vol = monitor_sum * vx # v

	endif else begin

		size_int = size( int, /dimensions )
		vol = dblarr( size_int[0], size_int[1] ) + 1d

	endelse

	err = sqrt( int ) / vol
	int = int / vol

	print, 'max intensity is ', max( int )
	print, 'max error is ', max( err )

	return, self -> set( qty = int $
		, err = err $
		, vol = vol $
		, xdata = hw # vy $
		, xlabel = 'Energy (meV)' $
		, ydata = vx # ang $
		, ylabel = 'Q (' + labels(/inv_A) + ')' $
		)

end

;*******************************************************************************
function oDave::read_dcs_old_angle, file = files

	if n_elements( files ) eq 0 then begin
		files = dialog_pickfile( /read, /multiple_files $
			, title = 'Please select a DCS file to read.' $
			, filter = ['*.dcs', '*.dcs.gz'] $
			)
	endif

	if total( file_test(files, /read) ) ne n_elements( files ) then begin
		message, 'The DCS file does not exists.', /continue
		return, 0
	endif

	; Number of useful time channels is fixed at this number.
	n_time = 1000

	; Read files.
	trifenestra_dcs_readfiles, files, d
	trifenestra_dcs_det, ang
	n_ang = n_elements( ang )

	; Sort the angle in the ascending order.
	key_ang = sort( ang )
	ang = ang[ key_ang ]

	; Calculate time.
	DT = d.tchanlook[0, 0:n_time - 1 ]
		; Delta Time * 20.
	time2 = d.tsdmin + total( DT, /cumulative) * .05
	time1 =	time2 - DT * .05

	; Convert time to energy.
	Ei = 81.8 / d.ch_wl^2
	E1 = 81.8 / (3.956 * time1 / d.det_dis)^2
	E2 = 81.8 / (3.956 * time2 / d.det_dis)^2
	Ef = (E2 + E1 ) / 2d
	DE = E1 - E2
		; Delta Energy
	hw	=	Ei - Ef
	n_hw	=	n_time ; + 1

	; Prepare empty vectors.
	oneset = 1 + fltarr(n_ang)
	onepix = 1 + fltarr(n_hw)

	; Make hw & angle 2D arrays.
	hw = oneset # hw
	Ef = oneset # Ef
	angle = ang # onepix

	; Phase Volume
	DP = oneset # (DE)

	; Obtain monitor normalized scattering intensity.
	monitor = total( d.histohigh[0:d.ngtchan1-1, 3] ) / d.duration
	int	= double( transpose(d.histodata[0:(n_time - 1),*]) )
	error = sqrt(int)
	int = int / DP / monitor ;  d.detsum[916] ; / DE
	error = error / DP / monitor ; d.detsum[916] ; / DE
		; By dividing by DP, we obtain differential cross section.

	; Arrange int and error in the angle ascending order.
	int = int[ key_ang, * ]
	error = error[ key_ang, * ]

 	return, self -> set( qty=int $
 		, err=error $
 		, vol=DP $
 		, xdata=angle $
 		, ydata=hw $
 		, xlabel='Angle' $
 		, ylabel='hw' $
 		, xunits='degree' $
 		, yunits='meV' $
 		, histlabel='I' $
 		, histunits='arb. units' $
 		, Ei = Ei $
 		)
end
;*******************************************************************************
function oDave::read_dcs_old, file = files

	if not self -> read_dcs_old_angle( file = files ) then return, 0

	return, 1

end

;*******************************************************************************

function oDave::read_dcs_old2, file = files $
	, angle=angle, q=q $
	, raw=raw, normalize=normalize $
	, monitor=monitor, T = T, Ei = Ei, duration = duration, DT = DT

	if n_elements( angle ) eq 0 or keyword_set( q ) then q = 1
	if n_elements( raw ) eq 0 or keyword_set( normalize ) then normalize =1 $
		else normalize = 0

	if n_elements( files ) eq 0 then begin
		files = dialog_pickfile( /read $
			, /multiple_files $
			, title = 'Please select a DCS file to read.' $
			, filter = ['*.dcs', '*.dcs.gz'] $
			)
	endif

	if total( file_test(files, /read) ) ne n_elements( files ) then begin
		message, 'The DCS file does not exists.', /continue
		return, 0
	endif

	if not self -> read_dcs_raw( file = files $
		, normalize=normalize, Ei=Ei, monitor=monitor, T=T $
		, duration = duration, DT = DT) then	begin
			message, 'Error in reading the DCS file.', /continue
	endif


	; Make it I(Q, hw).
	self -> transpose
	self -> get, err = err
	print, 'max err is ', max(err)

	; By default, I(Q, hw) must be presented.
	if q then begin
		self -> angle2q, Ei = Ei
	endif

	return, 1

end

;function odave::read_dcs, file = file $
;	, DT = DT $
;	, wavelength = wavelength $
;	, _extra = _extra
;
;	s = read_dcs( file, _extra = _extra )
;
;	DT = s.DT
;	wavelength = s.wavelength
;
;	; Fill in m array.
;	dim_f = size( s.f, /dimensions )
;	m = replicate( s.monitor, dim_f[0], dim_f[1] )
;
;	return, self -> set( qty = s.f / s.monitor $
;		, err = sqrt( s.f ) / s.monitor $
;		, vol = m $
;		, xdata = s.q $
;		, ydata = s.hw )
;
;end

;*******************************************************************************
function oDave::read_self, file = file

	on_ioerror, io_error
	if not file_test(file, /read) then begin
		message, 'No file to read.', /continue
		return, 0
	endif

	self -> empty
	oTmp = self
	restore, filename = file
	obj_destroy, oTmp
	return, 1

io_error:
		message, 'I/O Error', /continue
		return, 0

end
;*******************************************************************************
function oDave::write_self, file = file

	on_ioerror, io_error

	save, self, /compress, filename = file
	return, 1

io_error:
	message, 'I/O Error', /continue
	return, 0

end
;*******************************************************************************
;23456789 123456789 123456789 123465798 123456789 123456789 123465789 123456789
;       10        20        30        40        50        60        70        80
;*******************************************************************************
;
; OTHER I/O METHODS
;
; self->plot_fy
; self->plot_fx
; self->print
; self->read_self( file = file )
; self->write_self( file = file )
;
;*******************************************************************************
pro oDave::print

	print, '-------------------- oDave'
	if not self -> valid() then return

	print, 'self.n_qty = ', self.n_qty
	print, 'self.histtype = ', self.histtype
	print, 'self.xtype = ', self.xtype
	print, 'self.ytype = ', self.ytype

	maxf = max( *self.qty , min = minf )
	print, 'Max = ', maxf, ' Min = ', minf

	max_iteration = min( [ 20, self.n_qty ] )
	for i = 0, max_iteration -1 do begin
		print, (*self.qty)[i], (*self.x)[i], (*self.y)[i]
	endfor

end
;*******************************************************************************
pro oDave::plot_fy, _extra=_extra

	sort_key = sort( *self.y )
	plot, (*self.y)[sort_key], (*self.qty)[sort_key], _extra=_extra
	errplot, (*self.y)[sort_key] $
		, (*self.qty)[sort_key] - (*self.err)[sort_key] $
		, (*self.qty)[sort_key] + (*self.err)[sort_key]


end
;*******************************************************************************
pro oDave::plot_fx, _extra=_extra

	sort_key = sort( *self.x )
	plot, (*self.x)[sort_key], (*self.qty)[sort_key], _extra=_extra
	errplot, (*self.x)[sort_key] $
		, (*self.qty)[sort_key] - (*self.err)[sort_key] $
		, (*self.qty)[sort_key] + (*self.err)[sort_key]

end
;*******************************************************************************
;
; VIEW method only displays the object if the data is in 2D.
;
pro oDave::view $

;	polygonal_plot, (*self.qty), (*self.x), (*self.y), _extra = _extra

	, color_table = color_table $

	, zrange = zrange $
	, xrange = xrange $
	, yrange = yrange $

	, ratio = ratio $
	, logscal = logscal $
	, white = white $
	, pos = pos $
	, print = print $
	, _extra = _extra

	; Set default parameter values.
	minx = min( *self.x, max = maxx )
	if n_elements( xrange ) ne 2 then xrange =[ minx, maxx ]

	miny = min( *self.y, max = maxy )
	if n_elements( yrange ) ne 2 then yrange= [ miny, maxy ]

	if n_elements( pos ) ne 4 then pos = [.2,.2,.85,.9]

	; Fill in the empty points.

	; if n_elements( white ) ne 1 then white = -1d30
	; white_idx = where( finite( *self.qty, /nan ) )
	; if white_idx[0] ne -1 then (*self.qty)[ white_idx ] = white
	; z_ = *self.qty
	; white_idx = where( z_ eq white )
	; if white_idx[0] ne -1 then z_[ white_idx ] = - white
	; minz = min( z_, max = maxz )
	minz = min( *self.qty, max = maxz )

	if n_elements( zrange ) ne 2 then zrange =  [ minz, maxz ]
	if n_elements( ratio ) eq 1 then ra = ratio else ra=1.

	plot, /nodata, findgen(1) $
  		, xstyle=1, ystyle=1, xrange=xrange, yrange=yrange	$
  		, pos=pos, _extra=_extra

 	vert_colors = bytscl( reform( (*self.qty), self.n_qty ) $
			, max = zrange[1] $
			, min = zrange[0] $
			, top = 254 )

	self.oPoly -> setproperty, vert_colors = vert_colors

	pos_device = convert_coord( $
				[ [ pos[0], pos[1] ], [ pos[2], pos[3] ] ]$
				, /normal, /to_device )

	self.oBuffer -> setProperty, dimensions = $
		[ pos_device[0,1] - pos_device[0,0] $
		, pos_device[1,1] - pos_device[1,0] ]

	if keyword_set( print ) then bkgd_color = 255 else bkgd_color = 0
	self.oView -> setproperty $
		, viewplane_rect = [ xrange[0], yrange[0] $
		, xrange[1] - xrange[0], yrange[1] - yrange[0] ] $
		, color = bkgd_color

	self.oBuffer -> draw, self.oView

	oImage = self.oBuffer -> read()
	oImage -> getproperty, data = data
	obj_destroy, oImage


;	device, decomposed = 0
	loadct, 39
	tv, data, pos[0], pos[1], /normal $
		, xsize = (pos[2] - pos[0]), ysize = (pos[3] - pos[1])

	loadct, 0
	plot $
		, title = self.histlabel $
		, xtitle = self.xlabel $
		, ytitle = self.ylabel $
		, /nodata , [0]	$
  		, xstyle = 1 $
  		, ystyle = 1 $
  		, /noerase $
  		, xrange = xrange $
  		, yrange = yrange $
  		, pos = pos $
  		, _extra=_extra

end
;*******************************************************************************
;23456789 123456789 123456789 123465798 123456789 123456789 123465789 123456789
;       10        20        30        40        50        60        70        80
;*******************************************************************************


;
; Data Format Conversion
;
;*******************************************************************************
;23456789 123456789 123456789 123456789 123456789 123456789 123456789 123456789
;       10        20        30        40        50        60        70        80
;*******************************************************************************
; cgrain.pro
;
;  smooth w/ consideration to empty elements.
;
pro oDave::cgrain_old, nxa,nya,white=white
	; INPUT
	; z = a 2D array
	; [nxa, nya] = size of the small rectangle that is used for rebinning
	; white = the value for the empty elements.
	;	If this is set, those empty elements are ignored.

z = self.qty

nx	=	n_elements( (*z)[*,0] )
	; Size along x direction
ny	=	n_elements( (*z)[0,*] )
	; Size along y direction
zave	=	(*z)
	; Initialize the output matrix.

if n_elements(white) ne 0 then $ ; If a value is designated for some empty points...
	for i=0, nx-1 do begin
 		for j=0,ny-1 do $
 			if (*z)[i,j] ne white then begin
  				nav	=	0 ; Number of contributing elements
  				tav	=	0 ; Total value
  				; The following iav and jav are indices
  				;  for a small rectangle surrounding [i, j]
  				for iav=i-nxa/2,i+nxa/2 do begin
   					for jav=j-nya/2,j+nya/2 do begin
    					if (iav ge 0 and iav le nx-1 and jav ge 0 and jav le ny-1) then $
							if z(iav,jav) ne white then  begin
								tav	=	tav+(*z)(iav,jav)
								nav	=	nav+1
    						endif
   					endfor
  				endfor
  				; Calculate the average and store it.
  				zave(i,j)	=	tav/nav
 			endif ; It was 'end' before.
 		; endfor j
	endfor $ ; i, It was 'end' before.
else $
	for i=0,nx-1 do begin
 		for j=0,ny-1 do begin
  			nav	=	0
  			tav	=	0
  			for iav=i-nxa/2,i+nxa/2 do begin
   				for jav=j-nya/2,j+nya/2 do begin
    				if (iav ge 0 and iav le nx-1 and jav ge 0 and jav le ny-1) then begin
						tav	= tav+(*z)(iav,jav)
						nav	= nav+1
    				endif
   				endfor
  			endfor
  			zave(i,j)	=	tav/nav
 		endfor
	endfor

	ptr_free, z
	self.qty = ptr_new( zave, /no_copy )

end
;*******************************************************************************
;23456789 123456789 123456789 123465798 123456789 123456789 123465789 123456789
;       10        20        30        40        50        60        70        80
;*******************************************************************************
;
; Rebin Methods
;
; 'rebin' employs a few steps.
; The first step is to properly divide the points in the designated range
; (XRange, YRange) into pre-sized (Dimension) array, p_array.
; Those points whose indices appear in an element of p_array are merged
; later on.
; Note that the resultant array has become 2 dimensional.
;
;
pro oDave::rebin $
	, XRange = XRange $
	, YRange = YRange $
	, Dimension = Dimension

	; Find the default values when no input is given.
	if n_elements( XRange ) ne 2 then begin
		xmax = max( *self.x , min = xmin )
		XRange = [ xmin, xmax ]
	endif
	if n_elements( YRange ) ne 2 then begin
		ymax = max( *self.y, min = ymin )
		YRange = [ ymin, ymax ]
	endif
	if n_elements( Dimension ) ne 2 then begin
		Dimension = [ 1, 1 ]
	endif

	; Create empty array of indices.
	p_array = ptrarr( Dimension[0], Dimension[1] )
	for i = 0, Dimension[0] - 1 do begin
	for j = 0, Dimension[1] - 1 do begin
		p_array[i, j] = ptr_new([-1])
	endfor
	endfor

	; Distribute points on the p_array.
	x_length = XRange[1] - XRange[0]
	y_length = YRange[1] - YRange[0]
	for i = 0l, self.n_qty - 1 do begin

		; Calculate x index for the point.
		if (*self.x)[i] eq XRange[1] then begin
			ix = Dimension[0] - 1
		endif else begin
			fx = ( (*self.x)[i] - XRange[0] ) / x_length * Dimension[0]
			ix = fix( fx )
		endelse

		; Calculate the y index for the point.
		if (*self.y)[i] eq YRange[1] then begin
			iy = Dimension[1] - 1
		endif else begin
			fy = ( (*self.y)[i] - YRange[0] ) / y_length * Dimension[1]
			iy = fix( fy )
		endelse

		; Assign the point to the calculated indices.
		if ix ge 0 and ix lt Dimension[0] $
		and iy ge 0 and iy lt Dimension[1] then begin
			new_idx = [ *p_array[ix, iy], i ]
			ptr_free, p_array[ix, iy]
			p_array[ix, iy] = ptr_new(new_idx, /no_copy)
		endif

	endfor

	; Create empty arrays to store the results.
	f = dblarr( Dimension[0], Dimension[1] )
	e = f
	p = f
	x = f
	y = f

	; Reduce the points.
	for i = 0, Dimension[0] - 1 do begin
	for j = 0, Dimension[1] - 1 do begin

		x[ i, j ] = x_length * ( i  + .5d ) / Dimension[0] + XRange[0]
		y[ i, j ] = y_length * ( j  + .5d ) / Dimension[1] + YRange[0]

		n_points = n_elements( *p_array[i, j] )
		if n_points le 1 then begin
				f[i, j] = !values.d_nan
				e[i, j] = !values.d_nan
				p[i, j] = 0d
				continue
		endif
			; Get rid of [ -1 ]'s.
			n_points = n_points - 1
			new_idx = (*p_array[i, j])[1:*]
			ptr_free, p_array[i, j]
			p_array[i, j] = ptr_new(new_idx, /no_copy)

			; Create sub arrays.
			p_ = (*self.vol)[ *p_array[i, j] ]
			f_ = (*self.qty)[ *p_array[i, j] ]
			e_ = (*self.err)[ *p_array[i, j] ]

			; Calculate the new values.
			new_p = total( p_ )
			new_f = total( p_ * f_ ) / new_p
			new_e = sqrt( total( p_^2 * f_^2 ) ) / new_p

			f[ i, j ] = new_f
			e[ i, j ] = new_e
			p[ i, j ] = new_p
	endfor
	endfor

	; Clean up
	ptr_free, p_array

 	void = self -> set( qty = f $
 		, err = e $
 		, vol = p $
 		, xdata = x $
 		, ydata = y $
 		)

end
;*******************************************************************************

;*******************************************************************************
pro oDave::rebin_grid $
	; INPUT
	, Group = Group

	; Set default values.
	if n_elements( Group ) ne 2 then Group = [ 2, 2 ]
	if Group[0] le 0 then Group[0] = 1
	if Group[1] le 0 then Group[1] = 1


	; Check the dimension of the data array.
	size_f = size( *self.qty )
	if size_f[0] ne 2 then begin
		message, 'Unable to comply. F array is not 2D.', /continue
		return
	endif

	; Calculate the size of the resultant arrays.
	dimension_x = size_f[1] / Group[0]
	dimension_y = size_f[2] / Group[1]

	; Create empty arrays.
	f = dblarr( dimension_x, dimension_y )
	e = f
	p = f
	x = f
	y = f

	; Start the big loop.
	for i = 0, dimension_x - 1 do begin
	for j = 0, dimension_y - 1 do begin

		; Find the indices of those elements associated with i & j.
  		ix = Group[0] * i
  		iy = Group[1] * j

		; Create sub arrays.
		p_ = (*self.vol)[ ix:ix+Group[0]-1, iy:iy+Group[1]-1 ]
		f_ = (*self.qty)[ ix:ix+Group[0]-1, iy:iy+Group[1]-1 ]
		e_ = (*self.err)[ ix:ix+Group[0]-1, iy:iy+Group[1]-1 ]

		; Calculate the new values.
		new_p = total( p_ )
		new_f = total( p_ * f_ ) / new_p
		new_e = sqrt( total( p_^2 * e_^2 ) ) / new_p
		new_x = (*self.x)[ix, iy] + (*self.x)[ix+Group[0]-1,iy+Group[1]-1]
		new_x = new_x / 2.
		new_y = (*self.y)[ix, iy] + (*self.y)[ix+Group[0]-1,iy+Group[1]-1]
		new_y = new_y / 2.

	   	f[ i, j ] = new_f
	   	e[ i, j ] = new_e
	   	p[ i, j ] = new_p
	   	x[ i, j ] = new_x
	   	y[ i, j ] = new_y

	endfor
	endfor

 	void = self -> set( qty = f $
 		, err = e $
 		, vol = p $
 		, xdata = x $
 		, ydata = y $
 		)

end

pro oDave::rebin_old, dx,dy

	x = self.x
	y = self.y
	z = self.qty
	dz = self.err

;
; input
; x,y : one-dimensional vectors
; z,dz : two-D arrays
; dx,dy : required steps of x- and y-grids
;
; output
; xr,yr : resulting x and y variables with step dx and dy, respectively
; zr, dzr : averaged z and dz (error for z) for the new (xr,yr) points

; Calculate the number of elements for the output vectors.
	nxr	=	fix((max(*x)-min(*x))/dx)
	nyr	=	fix((max(*y)-min(*y))/dy)

; Generate the output vectors.
	xr	=	findgen(nxr)*dx+min(*x)+dx/2.
	yr	=	findgen(nyr)*dy+min(*y)+dy/2.

; print, 'dx = ', dx, ' dy = ', dy
; help, yr
; print, yr


	; Create empty output arrays.
	zr	=	dblarr(nxr,nyr)
	dzr	=	dblarr(nxr,nyr)

	; Create empty intermediate arrays.
	ny	=	n_elements((*y)[0, *])
	zrt	=	fltarr(nxr,ny)
	dzrt	=	fltarr(nxr,ny)

	; Create an empty array that will store number of indices
	;	used for averaging.
	facx	=	fltarr(nxr)

; Sum over x direction.
for i=0,nxr-1 do begin

	; Indices that fall into the right x range.
	lx	=	where( (*x)[*,0] ge xr(i)-dx/2. and (*x)[*,0] lt xr(i)+dx/2.,nx)
	facx(i)	=	nx

	zrt[i,*]	=	total((*z)([lx],*),1)
	dzrt[i,*]	=	total((*dz)([lx],*)^2,1)

end

; Sum over y direction and calculate average.
for j=0,nyr-1 do begin

	; Indices that fall into the right y range.
	ly	=	where((*y)[0,*] ge yr(j)-dy/2. and (*y)[0,*] lt yr(j)+dy/2.,ny)
	if ny eq -1 then begin
		message, 'increase energy step', /continue
		return
	endif

 	if ny eq 1 then begin
		zr(*,j)	=	zrt(*,ly)/facx/ny
		dzr(*,j)=	sqrt(dzrt(*,ly))/facx/ny
	endif

	if ny gt 0 then begin
		zr(*,j)	=	total( zrt(*,ly),2 )/facx/ny
		dzr(*,j)=	sqrt(total( dzrt(*,ly),2 ))/facx/ny
	endif

endfor

	vol = zr
	vol = vol - vol + 1d

	void = self -> set( $
		 qty = zr $
		, err = dzr $
		, vol = vol $
		, xdata = xr $
		, ydata = yr $
	)

end
;*******************************************************************************
;23456789 123456789 123456789 123465798 123456789 123456789 123465789 123456789
;       10        20        30        40        50        60        70        80


function oDave::cut $
	, XRange = XRange, YRange = YRange
	; RETURN
	;	-1 : if self is not valid.
	;	ooDave

	; Self not valid? Then return -1.
	if not self->valid() then return, -1

	; If XRange is not given, use the total x range as default.
	if n_elements( XRange ) ne 2 then begin
		xmin = min( *self.x, max = xmax )
		XRange = [ xmin, xmax ]
	endif

	; If YRange is not given, use the total y range as default.
	if n_elements( YRange ) ne 2 then begin
		ymin = min( *self.y, max = ymax )
		YRange = [ ymin, ymax ]
	endif

	; Select those that fall into the given range.
	idx = where( *self.y ge YRange[0] and *self.y le YRange[1] $
		and *self.x ge XRange[0] and *self.x le XRange[1] )

	if idx[0] eq -1 then return, -1

	f = (*self.qty)[idx]
	e = (*self.err)[idx]
	p = (*self.vol)[idx]
	x = (*self.x)[idx]
	y = (*self.y)[idx]

	return, obj_new( obj_class(self), qty=f, err=e, vol=p, xdata=x, ydata=y )

end

;*******************************************************************************
;23456789 123456789 123456789 123465798 123456789 123456789 123465789 123456789
;       10        20        30        40        50        60        70        80
;*******************************************************************************
; Reduce_Dimension
;
;	For a 2D data set, cut those portion that belongs to the given range,
;	and merge those points that are in the same row or column.
;	Note that they do not necessarily have the same y- or x- coordinate
;	even if they belong to the same row or column.
;
pro oDave::reduce_dimension, xrange=xrange, yrange=yrange

	qty = *self.qty
	vol = *self.vol
	err = *self.err


	if n_elements( xrange ) eq 2 then begin

		idx = where( *self.x ge xrange[0] and *self.x le xrange[1], /complement )
		qty[idx] = 0d
		vol[idx] = 0d

		vol_vec = total( vol, 1 )
		qty_vec = total( vol * qty, 1 ) / vol_vec
		err_vec = sqrt( total( (vol * err)^2, 1 ) ) / vol_vec
		; x_vec =


	endif

	if n_elements( yrange ) ne 2 then begin

		idx = where( *self.y ge yrange[0] and *self.y le yrange[1], /complement )
		qty[idx] = 0
		x_vec = total( qty, 2 )

	endif

end
;*******************************************************************************
;23456789 123456789 123456789 123465798 123456789 123456789 123465789 123456789
;       10        20        30        40        50        60        70        80
;*******************************************************************************
pro oDave::transpose

	void = self->set( qty = transpose( *self.qty ) $
		, err = transpose( *self.err ) $
		, vol = transpose( *self.vol ) $
		, xdata = transpose( *self.y ) $
		, ydata = transpose( *self.x ) $
		, xlabel = self.ylabel $
		, ylabel = self.xlabel $
  		, xunits = self.yunits $
  		, yunits = self.xunits $
		)

end
;*******************************************************************************
;23456789 123456789 123456789 123465798 123456789 123456789 123465789 123456789
;       10        20        30        40        50        60        70        80
;*******************************************************************************
;
;
; self -> angle_powder_correction
;	For powder samples, the intensity near the direct beam (angle = 0)
;	is larger than expected because of the finite and uniform detector
;	size.
;
pro oDave::powder_correction $
	, transpose = transpose $
	, _extra = _extra

	if keyword_set( transpose ) then begin
		angle = self.y
		hw = self.x
	endif else begin
		angle = self.x
		hw = self.y
	endelse

	*self.qty = *self.qty * abs( sin( *angle * !dtor ) )
	*self.err = *self.err * abs( sin( *angle * !dtor ) )
	; Phase volume does not need to be changed
	; because the calculated assumed non-powder sample.

end

pro oDave::angle2q $
	, Ei = Ei $
	, transpose = transpose $
	, kikf = kikf $
	, _extra = _extra

	if n_elements( Ei ) eq 0 then Ei = self.Ei

	if keyword_set( transpose ) then begin
		angle = *self.y
		q = self.y
		hw = self.x
	endif else begin
		angle = *self.x
		q = self.x
		hw = self.y
	endelse

	Ef = Ei - *hw
	*q = sqrt( $
		( Ei + Ef - 2 * cos(angle * !dtor) * sqrt(Ei * Ef) )$
		/ 2.0717 $
		)  * angle / abs( angle )

	; jacobian = abs( *q * 2.0717 / sqrt( Ei * Ef ) / sin( angle * !dtor ) )
	jacobian = 1d


	i_factor = 1d
	if keyword_set( kikf ) then begin
		i_factor = sqrt( Ei / Ef )
	endif

	; Conversion to wavevector space.
	*self.qty = *self.qty * jacobian * i_factor
	*self.err = *self.err * jacobian * i_factor
	*self.vol = *self.vol / jacobian

	; Q needs to be X. -> must be fixed.
	key = sort( (*q)[ *, 0 ] )
	*self.x = (*self.x)[ key, * ]
	*self.qty = (*self.qty)[ key, * ]
	*self.err = (*self.err)[ key, * ]
	*self.vol = (*self.vol)[ key, * ]
	*self.y = (*self.y)[ key, * ]

	if obj_valid( self.oPoly ) then obj_destroy, self.oPoly
	void = self -> _self_check()



	; Conversion to wavevector space.
;	int = int * p_factor
;	error = error * p_factor
;	DP = DP / p_factor

end
;*******************************************************************************
;23456789 123456789 123456789 123465798 123456789 123456789 123465789 123456789
;       10        20        30        40        50        60        70        80
;*******************************************************************************
pro odave::t3_background, yrange = yrange, DT = DT

	nx = n_elements( (*self.y)[*,0] )
	average = dblarr( nx )

	for i = 0, nx - 1 do begin

		idx = where( (*self.y)[i,*] le yrange[0] or (*self.y)[i,*] ge yrange[1] )
		duration = total( DT[ idx ] )
		intensity = total( (*self.qty)[i,idx], 2 )
		average[i] = double( intensity ) / duration

	endfor

	(*self.qty) = average # DT

end

;*******************************************************************************
;23456789 123456789 123456789 123465798 123456789 123456789 123465789 123456789
;       10        20        30        40        50        60        70        80
;*******************************************************************************
pro odave::t3_sub, yrange = yrange, DT = DT

	nx = n_elements( (*self.y)[*,0] )
	average = dblarr( nx )

	for i = 0, nx - 1 do begin

		idx = where( (*self.y)[i,*] le yrange[0] or (*self.y)[i,*] ge yrange[1] )
		duration = total( DT[ idx ] )
		intensity = total( (*self.qty)[i,idx], 2 )
		average[i] = double( intensity ) / duration

		; (*self.qty)[i,*] = (*self.qty)[i,*] - average[i] * DT

	endfor

	(*self.qty) = (*self.qty) - average # DT


end
;*******************************************************************************
;23456789 123456789 123456789 123465798 123456789 123456789 123465789 123456789
;       10        20        30        40        50        60        70        80
;*******************************************************************************
;
; ARITHMETIC
;
;	self->mul, obj
;	self->add, obj
;	self->sub, obj
;
; Note on dimensions:
;	Self's dimension prevails in all cases.
;

;*******************************************************************************
pro oDave::add $
	, obj

	if obj_valid( obj ) then begin
	if obj_class( obj ) eq obj_class( self ) then begin

		if self.n_qty eq 0 then begin

			void = self -> set( obj )

		endif

		*self.qty = *self.qty + *obj.qty
		*self.err = sqrt( (*self.err)^2 + (*obj.err)^2 )

		return

	endif
	endif

	*self.qty = *self.qty + obj
	return

end
;*******************************************************************************
pro oDave::sub $
	; INPUT
	, obj
		; A number, an array, or an object

	if obj_valid( obj ) then begin
	if obj_class( obj ) eq obj_class( self ) then begin
	if obj -> valid() then begin

		*self.qty = *self.qty - *obj.qty
		*self.err = sqrt( (*self.err)^2 + (*obj.err)^2 )
		return

	endif
	endif else return
	endif

	*self.qty = *self.qty - obj
	return

end

pro oDave::mul $
	, obj

	if obj_valid( obj ) then begin
	if obj_class( obj ) eq obj_class( self ) then begin

		message, 'Multiplication Not Implemented.', /continue
		return

	endif
	endif

	*self.qty = *self.qty * obj
	*self.err = *self.err * obj

	return

end


;*******************************************************************************
;23456789 123456789 123456789 123465798 123456789 123456789 123465789 123456789
;       10        20        30        40        50        60        70        80
;*******************************************************************************
;
; BASIC METHODS 2
;
; COMPARE
;	- incomplete
;
function oDave::compare, oDave

	if not obj_valid( oDave ) then begin
		message, 'Input is not an object.', /continue
		return, 0
	endif

	if obj_class( oDave ) ne obj_class( self ) then begin
		message, 'Input is not the same object as self.', /continue
		return, 0
	endif

	size_self = size( *self.qty, /dimensions )
	size_arg = size( *oDave.qty, /dimensions )

	n_size_self = n_elements( size_self )
	if n_size_self ne n_elements( size_arg ) then begin
		message, 'Input has different dimensionality.', /continue
		return, 0
	endif

	if total( size_self eq size_arg ) ne n_size_self then begin
		message, 'Input has different size.', /continue
		return, 0
	endif

	return, 1

end
;*******************************************************************************
;
;
function oDave::sum, xrange=xrange, yrange=yrange, index=index

	if n_elements( xrange ) ne 2 then begin
		maxx = max( *self.x, min = minx )
		xrange = [ minx, maxx ]
	endif

	if n_elements( yrange ) ne 2 then begin
		maxy = max( *self.y, min = miny )
		yrange = [ miny, maxy ]
	endif

	if n_elements( index ) eq 0 then index = lindgen( n_elements( *self.qty ) )

	idx = where( $
		(*self.x)[ index ] ge xrange[0] and $
		(*self.x)[ index ] le xrange[1] and $
		(*self.y)[ index ] ge yrange[0] and $
		(*self.y)[ index ] le yrange[1] )

	sum = total( (*self.qty)[ idx ] )
	error = sqrt( total( (*self.err)[ idx ]^2 ) )

	return, [ sum, error ]

end
;*******************************************************************************
;
;
function oDave::integrate_2d, xrange = xrange, yrange = yrange


	idx = where( *self.x ge xrange[0] and *self.x le xrange[1] and $
		*self.y ge yrange[0] and *self.y le yrange[1] )

	int = integrate_2d_points( *self.qty, *self.err, *self.x, *self.y $
		, index = idx )

	return, int

end
;*******************************************************************************
;
; From time to time, we need to integrate values over x or y axis or both.
; How do we do that? It is done through rebinning and multiplying the range
; of integration.
;
function oDave::integrate_2d_old $
	, xrange = xrange $
	, yrange = yrange $
	, _extra = _extra

	; Assuming POINTS.
	idx = where( $
		*self.x ge xrange[0] and *self.x le xrange[1] $
		and *self.y ge yrange[0] and *self.y le yrange[1] $
		and *self.vol ne 0d )

	; If there are no points in the region, return [ NaN, NaN, 0d ]
	if idx[0] eq -1 then goto, return_null

	; Rebin the points into a point.
	vol = total( (*self.vol)[idx] )
	qty = total( (*self.vol)[idx] * (*self.qty)[idx] ) / vol
	err = sqrt( $
		total( ( (*self.vol)[idx] * (*self.err)[idx] )^2) $
		) / vol

	; Calculate Area
	area = ( xrange[1] - xrange[0] ) * ( yrange[1] - yrange[0] )

	; Return the integrated value.
	return, [ qty, err ] * area

return_null:

	return, [ !values.d_nan, !values.d_nan ]

end
;*******************************************************************************
;
; NOTE
; If not valid, the pointers must be free.
;
function oDave::valid

	if self.n_qty eq 0 then return, 0
	return, 1

end
;*******************************************************************************
;
; SET method
;
; If any keyword argument is present, self value is set with it.
; After all keywords are processed, consistency is checked.
; Returns 1 if all is O.K., 0 otherwise.
;
function oDave::set $
	, oDave $
	, qty=qty, err=err, vol=vol, xdata=x, ydata=y $
	, xlabel=xlabel, ylabel=ylabel, histlabel=histlabel $
	, xunits=xunits, yunits=yunits, histunits=histunits $
	, instrument=instrument $
	, treatment=treatment $
	, Ei = Ei $
	, _extra = _extra

	if obj_valid( oDave ) eq 1 then begin

		if not oDave -> valid() then begin
			message, 'The given oDave is not valid.', /continue
			self -> empty
			return, 0
		endif

		oDave -> get $
			, qty = qty, err = err, vol = vol, xdata = x, ydata = y $
			, histlabel = histlabel, xlabel = xlabel, ylabel = ylabel $
			, histunits = histunits, xunits = xunits, yunits = yunits $
			, instrument = instrument, Ei = Ei

	endif

	if n_elements( Ei ) eq 1 then self.Ei = Ei

	n_x = n_elements( x )
	n_y = n_elements( y )
	n_qty = n_elements( qty )
	n_err = n_elements( err )
	n_vol = n_elements( vol )

	if n_qty ne 0 then begin
		if ptr_valid( self.qty ) then ptr_free, self.qty
		self.qty = ptr_new( qty )
		self.n_qty = n_qty
		if obj_valid( self.oPoly ) then obj_destroy, self.oPoly
	endif

	if n_x ne 0 then begin
		if n_x eq 1 then x = [ x ]
		if ptr_valid( self.x ) then ptr_free, self.x
		self.x = ptr_new( x )
		self.n_x = (size( x ))[1]
		if obj_valid( self.oPoly ) then obj_destroy, self.oPoly
	endif

	if n_y ne 0 then begin
		if n_y eq 1 then y = [ y ]
		if ptr_valid( self.y ) then ptr_free, self.y
		self.y = ptr_new( y )
		size_y = size( y )
		self.n_y = (size_y)[ size_y[0] < 2 ]
		if obj_valid( self.oPoly ) then obj_destroy, self.oPoly
	endif

	if n_err ne 0 then begin
  		if ptr_valid( self.err ) then ptr_free, self.err
  		self.err = ptr_new( err )
  		self.n_err = n_err
  	endif

	if n_vol ne 0 then begin
  		if ptr_valid( self.vol ) then ptr_free, self.vol
  		self.vol = ptr_new( vol )
  		self.n_vol = n_vol
  	endif

	if n_elements(xlabel) eq 1 then	self.xlabel = xlabel
	if n_elements(ylabel) eq 1 then	self.ylabel = ylabel
	if n_elements(histlabel) eq 1 then	self.histlabel = histlabel
	if n_elements(xunits) eq 1 then	self.xunits = xunits
	if n_elements(yunits) eq 1 then	self.yunits = yunits
	if n_elements(histunits) eq 1 then	self.histunits = histunits
	if n_elements(instrument) eq 1 then self.instrument = instrument
	if n_elements(treatment) eq 1 then self.treatment = treatment


	return, self -> _self_check()

end
;*******************************************************************************
pro oDave::get $
	, dimension = dimension $
	, n=n $
	, qty=qty, err=err, vol=vol, xdata=x, ydata=y $
	, xlabel = xlabel, ylabel = ylabel, histlabel = histlabel $
	, xunits = xunits, yunits = yunits, histunits = histunits $
	, xrange = xrange, yrange = yrange, frange = frange $
	, T = T, Ei = Ei, instrument=instrument $
	, treatment=treatment

	if not self->valid() then return

	n = self.n_qty

	if ptr_valid( self.qty ) then qty = *self.qty
	if ptr_valid( self.err ) then err = *self.err
	if ptr_valid( self.vol ) then vol = *self.vol
	if ptr_valid( self.x ) then x = *self.x
	if ptr_valid( self.y ) then y = *self.y

	T = self.T
	Ei = self.Ei

	xlabel = self.xlabel
	ylabel = self.ylabel
	histlabel = self.histlabel
	xunits = self.xunits
	yunits = self.yunits
	histunits = self.histunits
	instrument = self.instrument

	minx = min( *self.x, max = maxx )
	xrange = [minx, maxx]
	miny = min( *self.y, max = maxy )
	yrange = [miny, maxy]
	minf = min( *self.qty, max = maxf )
	frange = [minf, maxf]

	return

end
;*******************************************************************************
function oDave::get_dimensions

	if ptr_valid( self.qty ) then return, size( *self.qty, /dimensions )

end
;*******************************************************************************
pro oDave::get_label $
	, xlabel=xlabel, ylabel=ylabel, histlabel=histlabel

	xlabel = self.xlabel
	ylabel = self.ylabel
	histlabel = self.histlabel

end
;*******************************************************************************
pro oDave::get_range $
	, xrange=xrange, yrange = yrange, frange=frange

	if not self->valid() then return

	xmax = max( *self.x, min = xmin )
	xrange = [ xmin, xmax ]

	ymax = max( *self.y, min = ymin )
	yrange = [ ymin, ymax ]

	; Find the range only from finite values.
	fmax = max( (*self.qty)[ where( finite( *self.qty ) ) ], min = fmin )
	frange = [ fmin, fmax ]

	return

end
;*******************************************************************************
;
; Empties the object if full.
;
pro oDave::empty

	if ptr_valid( self.qty ) then ptr_free, self.qty
	if ptr_valid( self.err ) then ptr_free, self.err
	if ptr_valid( self.vol ) then ptr_free, self.vol
	if ptr_valid( self.x ) then ptr_free, self.x
	if ptr_valid( self.y ) then ptr_free, self.y

	if obj_valid( self.oBuffer ) then obj_destroy, self.oBuffer
	if obj_valid( self.oView ) then obj_destroy, self.oView
	if obj_valid( self.oPalette ) then obj_destroy, self.oPalette

	self.n_qty = 0
	self.n_err = 0
	self.n_vol = 0
	self.n_x = 0
	self.n_y = 0


end
;*******************************************************************************
;23456789 123456789 123456789 123465798 123456789 123456789 123465789 123456789
;       10        20        30        40        50        60        70        80
;*******************************************************************************
;
; INTERNAL METHODS
;
; _DISTRIBUTE
;
;	distributes points into several bins defined by
;
function oDave::_distribute $
	, x, y

end
; _MERGE
;
;	merges points given in the one dimensional 'index' array.
;	'vol' is the fraction of the phase volume involved from each point.
;	If not given, 'vol' is set to *self.vol, i.e., the whole point participates.
;
function oDave::_merge $
	; INPUT
	, index $
	, vol $
	, _extra=_extra

	; Initialize.
	n_index = n_elements( index )
	n_vol = n_elements( vol )
	if n_vol eq 0 then begin
		vol = (*self.vol)[ index ]
		n_vol = n_index
	endif

	; 'fraction' must have the same size as the 'index'.
	if n_fraction ne n_index then begin
		message, 'n_fraction is not identical as n_index.', /continue
		return, 0
	endif

	f = total( vol * *self.qty[ index ] )
	e = sqrt( total( ( vol * *self.err[ index ] )^2 ) )
	v = total( vol )

	return, [ f, e, v ]

end
;*******************************************************************************

function oDave::_dave2self $
	; INPUT
	, davePtr = davePtr

	if not dave_valid( davePtr ) then begin
		message, 'davePtr failed dave_valid check.', /continue
		return, 0
	endif else begin
		message, 'davePtr passed rigorous test.', /continue
	endelse

	self -> empty

	dataStrPtr = (*davePtr).dataStrPtr
	histPtr = (*dataStrPtr).commonStr.histPtr

 	return, self -> set( qty = (*histPtr).qty $
 		, err = (*histPtr).err $
 		, xdata = (*histPtr).x $
 		, ydata = (*histPtr).y $
 		, instrument = (*dataStrPtr).commonStr.instrument $
		, xlabel = (*dataStrPtr).commonStr.xlabel $
		, ylabel = (*dataStrPtr).commonStr.ylabel $
		, histLabel = (*dataStrPtr).commonStr.histLabel $
		, xunits = (*dataStrPtr).commonStr.xunits $
		, yunits = (*dataStrPtr).commonStr.yunits $
		, histUnits = (*dataStrPtr).commonStr.histUnits $
		, treatment = *((*dataStrPtr).commonStr.treatmentPtr) $
		)

end
;*******************************************************************************
function oDave::_self2dave $
	; INPUT, OPTIONAL
	, dummy = dummy $
		; Creates a dummy DAVE structure for future use.
	; OUTPUT
	, davePtr = davePtr


	; Put data into a DAVE structure.
	if keyword_set( dummy ) then begin
		histStr = { qty:[0,1], err:[1,1], x:[0,1], y:[0,1] }
	endif else begin
		if not self->valid() then return, 0
		histStr = {qty:*self.qty, err:*self.err, x:*self.x, y:*self.y}
	endelse

	if not dave_valid( davePtr ) then dave_makedaveptr, daveptr = daveptr
	*((*(*davePtr).dataStrPtr).commonStr.histPtr) = histStr
	return, 1

	commonStr = $
		{ instrument:self.instrument $
		, histPtr:ptr_new(histStr, /no_copy) $
		, xlabel: self.xlabel $
		, ylabel: self.ylabel $
		, histlabel: self.histlabel $
		, xunits: self.xunits $
		, yunits: self.yunits $
		, histUnits: self.histunits $
		, xtype: self.xtype $
		, ytype: self.ytype $
		, treatmentPtr:ptr_new( self.treatment ) $
		, commonStrPtr:ptr_new() $
		}

	dataStr = $
		{ commonStr:commonStr $
		, specificPtr:ptr_new( /allocate_heap ) $
		}

	state = $
		{ dataStrPtr:ptr_new( dataStr, /no_copy ) $
		}

	davePtr =  ptr_new( state, /no_copy )

	return, 1


end


;*******************************************************************************

;*******************************************************************************
function oDave::_self_check

	; If there is no data, the object is going to be emptied.
	; However, it is not an error.
	if self.n_qty eq 0 then begin
		self -> empty
		return, 1
	endif

	; If we have come this far, self.n_qty is not zero.
	; Error and phase volume is not necessary.
	switch 1 of
		self.n_x eq 0:
		self.n_y eq 0: begin
			err_message = 'Either X, Y, or QTY is empty.'
			goto, input_error
		end
	endswitch

	; Change X & Y into 2D if necessary.
	; This replaces the commented out part later.
	if (size( *self.qty ))[0] eq 2 then begin
		void = array_oned2twod( *self.x, *self.y, x, y )
		*self.x = x
		*self.y = y
		self.n_x = (size( *self.x ))[1]
		self.n_y = (size( *self.y ))[2]
	endif

	size_qty = size( *self.qty )
	size_x = size( *self.x )
	size_y = size( *self.y )
	self.histtype = size_qty[0]

	; If exist, err and vol must be of the same dimension and size as qty.
	if self.n_err ne 0 then begin
		size_err = size( *self.err )
		if 	size_qty[0] ne size_err[0] then goto, input_error
		if total( size_qty[ 0:size_qty[0] ] eq size_err[ 0:size_err[0] ] ) $
			ne size_qty[0] + 1 then begin
				err_message = 'ERR has different size than QTY.'
				goto, input_error
		endif
	endif else begin
		err_array = *self.qty - *self.qty + 1d
		self.err = ptr_new( err_array, /no_copy )
	endelse

	if self.n_vol ne 0 then begin
		size_vol = size( *self.vol )
		if	size_qty[0] ne size_vol[0] then goto, input_error
		if total( size_qty[ 0:size_qty[0] ] eq size_vol[ 0:size_vol[0] ] ) $
			ne size_qty[0] + 1 then begin
				err_message = 'VOL has different size than QTY.'
				goto, input_error
		endif
	endif else begin
		vol_array = *self.qty - *self.qty + 1d
		self.vol = ptr_new( vol_array, /no_copy )
	endelse



	; Fix X only if X is 1D and QTY is 2D.
;	if size_x[0] eq 1 and size_qty[0] eq 2 $
;		and size_y[0] ge 1 and size_y[0] le 2 then begin
;		y_vec = dblarr( size_y[ size_y[0] ] ) + 1d
;		*self.x = *self.x # y_vec
;		size_x = size( *self.x )
;	endif

	; Fix Y only if Y is 1D and QTY is 2D.
;	if size_y[0] eq 1 and size_qty[0] eq 2 $
;		and size_x[0] ge 1 and size_x[0] le 2 then begin
;		x_vec = dblarr( size_x[ 1 ] ) + 1d
;		*self.y = x_vec # *self.y
;		size_y = size( *self.y )
;	endif


	; Determine the data types.
	; X and Y must be of the same dimension as qty.
	; However, they may be bigger than qty in one dimension by 1.

	if size_x[0] ne size_qty[0] then begin
		err_message = 'X dimension does not match that of QTY.'
		goto, input_error
	endif

	if size_x[0] eq 0 then begin
		self.xtype = 'POINTS'
	endif else begin
		x_dim = size_x[1]
		qty_x_dim = size_qty[1]
		case (x_dim - qty_x_dim) of
			0: self.xtype = 'POINTS'
			1: self.xtype = 'HISTOGRAM'
			else: begin
				err_message = 'X size does not match that of QTY.'
				print, x_dim, qty_x_dim
				goto, input_error
			end
		endcase
	endelse

	if size_y[0] ne size_qty[0] then begin
		err_message = 'Y dimension does not match that of QTY.'
		goto, input_error
	endif

	if size_y[0] eq 0 then begin
		self.ytype = 'POINTS'
	endif else begin
		y_dim = ( size_y[0] ge 2 ) ? size_y[2] : 1
		qty_y_dim = ( size_qty[0] ge 2 ) ? size_qty[2] : 1

		case (y_dim - qty_y_dim) of
			0: self.ytype = 'POINTS'
			1: self.ytype = 'HISTOGRAM'
			else: begin
				err_message = 'Y size does not match that of QTY.'
				goto, input_error
			end
		endcase
	endelse

	; X and Y should have the same size.
	if total( size_x[ 0:size_x[0] ] eq size_y[ 0:size_y[0] ] ) $
			ne size_qty[0] + 1 then begin
				print, size_x
				print, size_y
				print, size_qty
				print, total(size_x[ 0:size_x[0] ] eq size_y[ 0:size_y[0] ] )
				err_message = 'X has different size than Y.'
				goto, input_error
	endif

	; Check the sizes of the remaining dimesions.
	if size_qty[0] gt 2 then begin
		x_eq_total = total( size_qty[3:size_qty[0]] eq size_x[2:size_x[0]] )
		y_eq_total = total( size_qty[3:size_qty[0]] eq size_y[2:size_y[0]] )
		if (x_eq_total ne size_qty[0] - 2) $
			or (y_eq_total ne size_qty[0] - 2) then begin
			err_message = 'The remaining dimensions do not match that of QTY.'
			goto, input_error
		endif
	endif

	; **************************************************************************
	; Prepare for the plot.

	; 1. Palette object
	if not obj_valid( self.oPalette ) then begin
		; device, get_decomposed = decomposed
		; device, decomposed = 0
		; loadct, 39
		loadct, 39
		tvlct, r, g, b, /get
		self.oPalette = obj_new('IDLgrPalette', r, g, b )
	endif

	; 2. Buffer object
	if not obj_valid( self.oBuffer ) then begin
		self.oBuffer = obj_new('IDLgrBuffer' $
			, color_model = 1, palette = self.oPalette )
	endif

	; 3. View object
	if not obj_valid( self.oView ) then begin
		self.oView = obj_new('IDLgrView', color = $
			( (!d.name eq 'PS') ? 255 : 0 ) )
			; Background color = 0
	endif

	; 4. Model object
	if not obj_valid( self.oModel ) then begin
		self.oModel = obj_new('IDLgrModel')
		self.oView -> add, self.oModel
	endif

	; 5. Polygon object
	if obj_valid( self.oPoly ) then obj_destroy, self.oPoly
;	if not obj_valid( self.oPoly ) then begin

		dim = size( *self.qty, /dimensions )
		if n_elements( dim ) ne 2 then return, 1

		case self.xtype of
			'POINTS': begin
				x = *self.x
				y = *self.y
			end
			'HISTOGRAM': begin
				x = (*self.x)[0:self.n_x-2,*] + (*self.x)[1:self.n_x-1,*]
				y = (*self.y)[0:self.n_x-2,*] + (*self.y)[1:self.n_x-1,*]
				x = x / 2d
				y = y / 2d
			end
		endcase

		case self.ytype of
			'POINTS':
			'HISTOGRAM': begin
				x = x[*,0:self.n_y-2] + x[*,1:self.n_y-1]
				y = y[*,0:self.n_y-2] + y[*,1:self.n_y-1]
				x = x / 2d
				y = y / 2d
			end
		endcase


		nx = n_elements( *self.qty )
		verts = dblarr( 3, nx )
		verts[0, *] = x
		verts[1, *] = y

		;vert_colors = bytscl( reform( z, nx ) $
		;	, max = zrange[1] $
		;	, min = zrange[0] $
		;	, top = 255 )

		conn = lonarr( 5 * ( dim[0] - 1l) * (dim[1] - 1l) )
		conn_idx = lindgen( ( dim[0] - 1l ) * (dim[1] - 1l) ) * 5l
		vert_idx = $
			reform( ( lindgen( dim[0], dim[1] ) )[0:dim[0]-2l,0:dim[1]-2l]  $
				, (dim[0]-1l)*(dim[1]-1l) )

		conn[ conn_idx ] = 4l
		conn[ conn_idx + 1 ] = 	0l + vert_idx
		conn[ conn_idx + 2 ] =  1l + vert_idx
		conn[ conn_idx + 3 ] =  1l + vert_idx + dim[0]
		conn[ conn_idx + 4 ] =  0l + vert_idx + dim[0]

		self.oPoly = obj_new('IDLgrPolygon' $
			, data = verts, polygons = conn $
			, shading = 1 )

		self.oModel -> add, self.oPoly

	; end

	return, 1

input_error:

	self -> empty

	message, '_self_check failed.', /continue
	if n_elements( err_message ) eq 1 then message, err_message, /continue
	message, 'Returning 0.', /continue
	return, 0

end
;*******************************************************************************
;23456789 123456789 123456789 123465798 123456789 123456789 123465789 123456789
;       10        20        30        40        50        60        70        80
;*******************************************************************************
function oDave::oPoly

	dim = size( *self.qty, /dimensions )
	if n_elements( dim ) ne 2 then return, 0

	case self.xtype of
		'POINTS': begin
			x = *self.x
			y = *self.y
		end
		'HISTOGRAM': begin
			x = (*self.x)[0:self.n_x-2,*] + (*self.x)[1:self.n_x-1,*]
			y = (*self.y)[0:self.n_x-2,*] + (*self.y)[1:self.n_x-1,*]
			x = x / 2d
			y = y / 2d
		end
	endcase

	case self.ytype of
		'POINTS':
		'HISTOGRAM': begin
			x = x[*,0:self.n_y-2] + x[*,1:self.n_y-1]
			y = y[*,0:self.n_y-2] + y[*,1:self.n_y-1]
			x = x / 2d
			y = y / 2d
		end
	endcase


	nx = n_elements( *self.qty )
	verts = dblarr( 3, nx )
	verts[0, *] = x
	verts[1, *] = y

	;vert_colors = bytscl( reform( z, nx ) $
	;	, max = zrange[1] $
	;	, min = zrange[0] $
	;	, top = 255 )

	conn = lonarr( 5 * ( dim[0] - 1l) * (dim[1] - 1l) )
	conn_idx = lindgen( ( dim[0] - 1l ) * (dim[1] - 1l) ) * 5l
	vert_idx = $
		reform( ( lindgen( dim[0], dim[1] ) )[0:dim[0]-2l,0:dim[1]-2l]  $
			, (dim[0]-1l)*(dim[1]-1l) )

	conn[ conn_idx ] = 4l
	conn[ conn_idx + 1 ] = 	0l + vert_idx
	conn[ conn_idx + 2 ] =  1l + vert_idx
	conn[ conn_idx + 3 ] =  1l + vert_idx + dim[0]
	conn[ conn_idx + 4 ] =  0l + vert_idx + dim[0]

	return, obj_new('IDLgrPolygon' $
		, data = verts, polygons = conn $
		, shading = 1 )

end
;*******************************************************************************
;23456789 123456789 123456789 123465798 123456789 123456789 123465789 123456789
;       10        20        30        40        50        60        70        80
;*******************************************************************************
;*******************************************************************************
;23456789 123456789 123456789 123465798 123456789 123456789 123465789 123456789
;       10        20        30        40        50        60        70        80
;*******************************************************************************
;
; DEFINITION & LIFECYCLE METHODS
;
;*******************************************************************************
function oDave::init, _extra = _extra

	self.n_qty = 0

	cd, current = current

	self.instrument = 'No instrument'
	self.treatment = ''
	self.xlabel = 'X axis'
	self.xunits = 'arb. units'
	self.xtype = 'POINTS'
	self.ylabel = 'Y axis'
	self.yunits = 'arb. units'
	self.ytype = 'POINTS'
	self.histlabel = 'Intensity'
	self.histunits = 'arb. units'
	self.histtype = 2

	return, self->set( _extra = _extra )

end
;*******************************************************************************
pro oDave::cleanup

	self->empty

end
;*******************************************************************************
pro oDave__define

	void = { oDave $

		; Size of the data
		, n_qty:0l $
		, n_err:0l $
		, n_vol:0l $
		, n_x:0l $
		, n_y:0l $

		; Data themselves
		, qty:ptr_new() $		; Functional Value
		, err:ptr_new() $		; Error
		, vol:ptr_new() $		; Phase Volume
		, x:ptr_new() $			; Horizontal Axis
		, y:ptr_new() $			; Vertical Axis

		, Ei:0d	$
		, T:0d $

		; Labels and Units
		, instrument:'' $
		, treatment:'' $

		, xlabel:'' $
		, xunits:'' $
		, xtype:'' $

		, ylabel:'' $
		, yunits:'' $
		, ytype:'' $

		, histlabel:'' $
		, histunits:'' $
		, histtype:0 $

		; Special Physical Quantities

		; Plot related objects
		, oBuffer:obj_new() $
		, oView:obj_new() $
		, oModel:obj_new() $
		, oPoly:obj_new() $
		, oPalette:obj_new() $

		}

end
;*******************************************************************************
