;+
; NAME:
;       ANG_BAR
;
; PURPOSE:
;	Computes the coordinates of the baricenter, or geometrical center, on the sphere,
;   of a set of points P.
;
; CATEGORY:
;
; CALLING SEQUENCE:
;
;	Result = ANG_BAR(P, [/DEGREES], [/CENTER])
;
; INPUTS:
;       P:  array (2,N), where N is the number of points of the form
;		[lon, lat] in RADIANS.
;
; KEYWORDS:
;
;       DEGREES: If set, accept P in degrees and return the baricenter position in degrees.
;
;       CENTER:  If set, return the geometrical center instead of the baricenter.
;
; OUTPUTS:
;   Result: vector of coordinates of baricenter [lon, lat] in RADIANS.
;	NOTE: if the points are distributed point-symmetrical way respect to the
;	center of the sphere, then no baricentrical points exist on the sphere
;	surface and a default value of [lon,lat]=[0,0] is returned.
;
; MODIFICATION HISTORY:
;          G. Li Causi, O.A.R.: 25-08-2000
;-

FUNCTION ang_bar, P, DEGREES=DEGREES, CENTER=CENTER

npoints = n_elements(P(0, *))
IF npoints EQ 1 THEN RETURN, P(*,0)			;se c'e' un solo punto

IF KEYWORD_SET(DEGREES) THEN P = P * !DTOR

;Conversione da coord. angolari a coord. x,y,z:
x = cos(P(0,*)) * cos(P(1,*))
y = sin(P(0,*)) * cos(P(1,*))
z = sin(P(1,*))

xbar = TOTAL(x)/npoints
ybar = TOTAL(y)/npoints
zbar = TOTAL(z)/npoints

IF KEYWORD_SET(CENTER) THEN BEGIN
	xbar = (max(x) + min(x)) / 2.
	ybar = (max(y) + min(y)) / 2.
	zbar = (max(z) + min(z)) / 2.
ENDIF

modulus = SQRT(xbar^2 + ybar^2 + zbar^2)
xbar = xbar / modulus
ybar = ybar / modulus
zbar = zbar / modulus

IF xbar EQ 0 AND ybar EQ 0 AND zbar EQ 0 THEN RETURN, [0,0]

;Conversione da coord. x,y,z, a coord. angolari:
lat = asin(zbar)
IF cos(lat) EQ 0 THEN lon = 0 ELSE lon = atan(ybar/cos(lat),xbar/cos(lat))

bar = [lon, lat]

IF KEYWORD_SET(DEGREES) THEN bar = bar * !RADEG

RETURN, bar

END