;$Id: integrate_2d_points.pro,v 1.2 2005/02/09 14:02:26 azuah Exp $
; Scheme 1.
; If there are 2 * 2 points, each point gets
; 1/4 of the area of the rectangle.
; The total volume becomes,
;
;	Area = ( R2 - R0 ) X ( R1 - R3 ) / 2,
;	Volume = ( Z1 + Z2 + Z3 + Z4 ) * Area / 4.
;
; Now consider the case where there are more than 2 * 2 points.
;
;
;
;###############################################################################
; 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.
;
;###############################################################################
pro integrate_2d_points_test

	xv = [ 0., 1., 2., 3., 4. ]
	yv = xv
	v = dblarr( 5 ) + 1d

	x = xv # v
	y = v # yv

	z = dblarr( size( x, /dimensions ) ) + 1d
	e = z

	print, integrate_2d_points( z, e, x, y )

end

function integrate_2d_points_mode1, z, e, x, y, index=index

	z_dims = size(z, /dimensions)
	; Fill in the area array.
	x0 = x[0:z_dims[0]-2, 0:z_dims[1]-2]
	x1 = x[1:*          , 0:z_dims[1]-2]
	x2 = x[1:*          , 1:*          ]
	x3 = x[0:z_dims[0]-2, 1:*          ]
	y0 = y[0:z_dims[0]-2, 0:z_dims[1]-2]
	y1 = y[1:*          , 0:z_dims[1]-2]
	y2 = y[1:*          , 1:*          ]
	y3 = y[0:z_dims[0]-2, 1:*          ]

	rect = ( x2 - x0 ) * ( y3 - y1 ) - ( x3 - x1 ) * ( y2 - y0 )
	; rect = rect / 2d

	; Prepare the area array that will be assigned to each point.
	area0 = dblarr( z_dims )
	area1 = dblarr( z_dims )
	area2 = dblarr( z_dims )
	area3 = dblarr( z_dims )

	area0[0:z_dims[0]-2, 0:z_dims[1]-2] = rect
	area1[1:z_dims[0]-1, 0:z_dims[1]-2] = rect
	area2[1:z_dims[0]-1, 1:z_dims[1]-1] = rect
	area3[0:z_dims[0]-2, 1:z_dims[1]-1] = rect

	area = area0 + area1 + area2 + area3

	; Mask
	mask = dblarr( z_dims )
	if n_elements( index ) ne 0 then begin
		mask[ index ] = 1d
	endif else begin
		mask = 1d
	endelse

	volume = total( area * z * mask )
	error = sqrt( total( area^2 * e^2 * mask ) )

	return, [ volume, error ] / 8d

end

;*******************************************************************************
;23456789 123456789 123456789 123465798 123456789 123456789 123465789 123456789
;       10        20        30        40        50        60        70        80
;*******************************************************************************
; Mode 0 does not calculate the volume exactly, but rather in a naive fasion
; assuming the points are distributed randomly in a RECTANGULAR area.
;
;

function integrate_2d_points_mode0, z, e, x, y, index = index

	z_dims = size(z, /dimensions)
	; Mask
	mask = dblarr( z_dims )
	if n_elements( index ) ne 0 then begin
		mask[ index ] = 1d
	endif else begin
		mask = 1d
		index = lindgen( n_elements( z ) )
	endelse

	xmax = max( x[ index ], min = xmin )
	ymax = max( y[ index ], min = ymin )
	area = ( xmax - xmin ) * ( ymax - ymin )

	volume = total( z * mask ) * area
	error = sqrt( total( e^2 * mask ) ) * area

	return, [ volume, error ] / 8d

end

function integrate_2d_points, z, e, x, y $
	, index = index $
	, mode = mode

	if n_elements( mode ) eq 0 then mode = 1
	case mode of
		0: return, integrate_2d_points_mode0( z, e, x, y, index=index )
		else: return, integrate_2d_points_mode1( z, e, x, y, index=index )
	endcase

end