;
;+
; NAME:
;
;       POLE_SPHERE_CIRCLE
;
; PURPOSE:
;
;		Computes the sherical coordinates of the north pole of a circle
;		for three points on a sphere or of the maximum circle for two points.
;
; CATEGORY:
;
;       Mathematics
;
; CALLING SEQUENCE:
;
;		POLE_SPHERE_CIRCLE, p1,p2,p3, lon_pol,lat_pol [, /north] [, /south]
;
; INPUT:
;
;		p1, p2, p3: 2-elements vectors of the form [lon, lat] of the coordinates of the
;				three points on the sphere (degrees).
;				If p3 is a zero scalar i.e. p3=0, then the maximum circle passing for p1 and p2
;				is considered.
;
; OUTPUT:
;
;		lon_pol, lat_pol: computed longitude and latitude of the north pole of the circle (degrees).
;				"North pole" means the pole from which the three points are viewed in clockwise order.
;
;
; KEYWORDS:
;
;		NORTH: if set, the pole in northern emisphere is computed
;
;		SOUTH: if set, the pole in southern emisphere is computed
;
;
; PROCEDURE:
;
;		The plane passing for the three points, or for p1, p2 and the sphere center, is computed.
;		The positive normal of the plane is the direction of the pole, whose [lon, lat] coordinates
;		are so computed.
;
;
; MODIFICATION HISTORY:
;
;       March 2004 - 	Gianluca Li Causi, INAF - Rome Astronomical Observatory
;					licausi@mporzio.astro.it
;					http://www.mporzio.astro.it/~licausi/
;-

PRO Pole_Sphere_Circle, p1,p2,p3, lon_pol,lat_pol, north=north, south=south

p1 = double(p1)
p2 = double(p2)
p3 = double(p3)

center = 0
IF n_elements(p3) EQ 1 THEN BEGIN
	center = 1
	p3 = [0,0]
ENDIF

p1r = p1 * !dtor	;To radians
p2r = p2 * !dtor
p3r = p3 * !dtor

lon = [p1r[0], p2r[0], p3r[0]]
lat = [p1r[1], p2r[1], p3r[1]]
x = cos(lat) * cos(lon)		;To xyz space
y = cos(lat)* sin(lon)
z = sin(lat)

IF center EQ 1 THEN BEGIN
	x[2] = 0.
	y[2] = 0.
	z[2] = 0.
ENDIF

;Plane passing for the three points:
aa = (y(1)-y(0))*(z(2)-z(0))-(y(2)-y(0))*(z(1)-z(0))
bb = -((x(1)-x(0))*(z(2)-z(0))-(x(2)-x(0))*(z(1)-z(0)))
cc = (x(1)-x(0))*(y(2)-y(0))-(x(2)-x(0))*(y(1)-y(0))
;;dd = aa*x(0)+bb*y(0)+cc*z(0)

;Positive normal to the plane
m = sqrt(aa^2+bb^2+cc^2)
nx = aa/m		;normale al piano
ny = bb/m
nz = cc/m

IF KEYWORD_SET(north) THEN BEGIN
	if nz lt 0 then begin
		nx = -nx
		ny = -ny
		nz = -nz
	endif
ENDIF
IF KEYWORD_SET(south) THEN BEGIN
	if nz gt 0 then begin
		nx = -nx
		ny = -ny
		nz = -nz
	endif
ENDIF


;coordinate del polo del cerchio massimo per i due punti
lat_pol = !radeg*asin(nz/(sqrt(nx^2+ny^2+nz^2)))
if nx eq 0 then lon_pol = 90+180*(ny lt 0) $
	else $
	lon_pol = !radeg*atan(ny, nx)

while lat_pol gt 90 do lat_pol = 180-lat_pol
while lat_pol lt (-90) do lat_pol = -180-lat_pol
lon_pol = lon_pol mod 360

END