;+
; NAME:
;	flat_centerofmass
; PURPOSE:
;	Calculations related to the centroid of a flat density distribution
; CATEGORY:
;	gen/idl/toolbox
; CALLING SEQUENCE:
	FUNCTION flat_centerofmass, box, centroid	, $
		polar		= polar 		, $
		shift_origin= shift_origin	, $
		degrees 	= degrees

;	R = flat_centerofmass(box [, /polar, /degrees])
;			Returns centroid of 'box'
;	R = flat_centerofmass(box ,/shift_origin [, /polar, /degrees])
;			Returns coordinate of 'box' relative to its centroid
;	R = flat_centerofmass(box , centroid [, /polar, /degrees])
;			Positions box with its centroid matching 'centroid'
; INPUTS:
;	If /polar NOT set:
;	box			array[2,2]; type: integer
;					defines two corners of box in the form [ [x1,y1], [x2,y2] ]
;
;	If /polar set:
;	box			array[2,2]; type: float
;					limiting values in phase angle and radius of the wedge
;					in the form [[angle1,radius1],[angle2,radius2]].
;					Angle1 and angle2 are in radians between [-!pi,+!pi].
;					The wedge runs counterclockwise from 'angle1' to 'angle2'
;					over less than 180 degrees: either angle2 > angle1 with
;					angle2-angle1 < !pi or angle2 < angle1 with
;					angle2+2*!pi-angle1 < !pi. Always radius1 < radius2.
; OPTIONAL INPUT PARAMETERS:
;	centroid	array[2]; type: float or integer
;					position of centroid in rectangular (if /polar NOT set)
;					or polar (if /polar set) coordinates. If this argument is specified
;					then box is positioned with its centroid on this point.
;	/polar		if set, box and centroid is specified in polar coordinates
;					(default: rectangular coordinates)
;	/shift_origin
;				if set, then the centroid of 'box' is calculated. Then the coordinates
;					of 'box' relative to this center are returned.
;	/degrees	if set, in- and output angles are in degrees (default: radians)
; OUTPUTS:
;	R			array[2]
;				array[2,2]
; INCLUDE:
	@compile_opt.pro			; On error, return to caller
; CALLS:
;	InitVar, IsType, AngleRange, ToRadians
; PROCEDURE:
;	For a rectangular, constant density distribution the centroid is the
;	geometric center of the box (the mean of the x and y coordinates of the corners).
;
;	For a wedge shaped box the phase angle of the centroid is the geometric mean
;	of the phase angles of the corners: phi(centroid) = (phi1+phi2)/2.
;	For the radius, r(centroid)^2 = (r1^2+r2^)/2.
;
;	If /shift_origin is set then the return values are the coordinates
;	of 'box' relative to the centroid (i.e. the centroid coordinates are
;	subtracted from the box coordinates

;	If argument 'centroid' is specified then the centroid of 'box' is calculated.
;	Then the box is positioned on 'centroid' such that the box centroid coincides
;	with 'centroid'.
; MODIFICATION HISTORY:
;	MAR-2000, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;-

InitVar, polar		  , /key
InitVar, shift_origin , /key
InitVar, get_centroid , /key
add2centroid = IsType(centroid, /defined)

IF polar THEN BEGIN								; Polar coordinates

	CASE add2centroid OF
	0: BEGIN

		rpm	= ToRadians(degrees=degrees)

		; box specifies the phase angle between [-!pi,!pi]. The phase angle
		; covered runs from box[0,0] to box[0,1]. The phase angle for the centroid
		; is the mean of the two, but some care must be taken when the box
		; straddles phase angle !pi.

		phi = box[0,*]							; Phase angle limits
		phi += (2*!pi/rpm)*[0B, phi[1] LT phi[0]]; In case the box straddles !pi
		phi_centroid = 0.5*total(phi)			; Phase angle of centroid

		; The radial location of the centroid follow from an equal area argument:
		; the area between the lower radius and the centroid must be the same as
		; between the centroid and the outer radius.

		rad = box[1,*]							; Radial limits
		rad_centroid = sqrt(0.5*total(rad*rad))	; Radius of centroid

		CASE shift_origin OF
		0: result = [AngleRange(phi_centroid, /pi, degrees=degrees), rad_centroid]
		1: result = [phi-phi_centroid, rad-rad_centroid]
		ENDCASE

	END
	1: BEGIN

		; Call dr1 = box[1] and dr2 = box[3]
		; Dropping the box at radius r the radial boundaries are:
		;	r1 = r+dr1, r2 = r+dr2
		; The centroid of this box is located at s determined by:
		;	s^2 = (r1^2+r2^2)/2 = r^2+r*(dr1+dr2)+(dr1^2+dr2^2)/2
		; The left hand side is the value stored in centroid[1]
		; Solving this for r (taking the positive root) gives
		;	r = -(dr1+dr2)/2+sqrt(s^2-(dr1-dr2)^2/4)
		; Dropping the box at r will put the centroid at the required value s.

		rad = sqrt(centroid[1]^2-0.25*(box[1]-box[3])^2) - 0.5*(box[1]+box[3])

		result = [centroid[0],rad]#[1,1]+box
		result[0,*] = AngleRange(result[0,*], /pi, degrees=degrees)

	END
	ENDCASE

ENDIF ELSE BEGIN								; Rectangular coordinates

	CASE add2centroid OF
	0: BEGIN
		xy_centroid = 0.5*total(box,2)

		CASE shift_origin OF
		0: result = xy_centroid
		1: result = box-xy_centroid#[1,1]
		ENDCASE
	END
	1: result = centroid#[1,1]+box
	ENDCASE

ENDELSE

RETURN, result  &  END