;+
; NAME:
;	sphere_distance
; PURPOSE:
;	Get angle between two vectors
; CATEGORY:
;	sat/idl/toolbox/math
; CALLING SEQUENCE:
	FUNCTION sphere_distance, Los1, Los2, Los3, Los4, Los5, Los6, $
		degrees 	= degrees		, $
		rectangular = rectangular	, $
		eps 		= eps
; INPUTS:
;	Los1, Los2		any arrays with 1st dimension of 2 or 3 elements
;					containing vector coordinates in spherical
;					(longitude, latitude, distance) or rectangular (x,y,z)
;					coordinates.
; OPTIONAL INPUT PARAMETERS:
;	/degrees		if set, all angles are in degrees (default: radians)
;	/rectangular	if set, Los1 and Los2 are in rectangular (x,y,z) coordinates
; OUTPUTS:
;	E				distances in range [0,180] degrees.
;					The array size is the same as for the input array with
;					the 1st dimension removed.
; OPTIONAL OUTPUT PARAMETERS:
;	eps=eps 		position angle of Los2 relative to north at Los2
;					(only works if /rectangular is NOT set)
; INCLUDE:
	@compile_opt.pro			; On error, return to caller
; CALLS:
;	ToRadians, SyncArgs, SyncDims
; EXAMPLE:
;	Specify vectors in spherical coordinates (all in the equatorial plane).
;	Every vector in los2 is combined with the vector in los1. All angles are in degrees.
;
;	los1 = [0,0,1]
;	los2 = [[0,0,1], [45,0,1], [90,0,1], [135,0,1], [180,0,1]]
;	e = sphere_distance(los1, los2, /degrees)
;	print, e
;		0.000000      45.0000      90.0000      135.000      180.000
;	help, los1, los2, e
;		LOS1            INT       = Array[3]
;		LOS2            INT       = Array[3, 5]
;		E               FLOAT     = Array[5]
; PROCEDURE:
;	Apply the following to the longitude and latitude angles:
;	L = acos( sin(lat1)*sin(lat2)+cos(lat1)*cos(lat2)*cos(lng1-lng2) )
; MODIFICATION HISTORY:
;	AUG-2000, Paul Hick (UCSD/CASS)
;	DEC-2000, Paul Hick (UCSD/CASS)
;		Modified to accept vector arrays of any number of dimensions
;	JUL-2006, Paul Hick (UCSD/CASS)
;		Added option to provide components of vectors seperately
;		(using 4 or 6 arguments).
;	JUL-2008, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Added keyword 'eps'. Renamed from 'elongation' to
;		'sphere_distance'
;-

rpm = ToRadians(degrees=degrees)

CASE n_params() OF
2: BEGIN
	L = Los1
	V = Los2
	InitVar, rectangular, /key
END
4: BEGIN
	L = SuperArray(Los1,2,/lead)
	L = SubArray(L,element=1,replace=Los2)
	V = SuperArray(Los3,2,/lead)
	V = SubArray(V,element=1,replace=Los4)
	rectangular = 0
END
6: BEGIN
	L = SuperArray(Los1,3,/lead)
	L = SubArray(L,element=1,replace=Los2)
	L = SubArray(L,element=2,replace=Los3)
	V = SuperArray(Los4,3,/lead)
	V = SubArray(V,element=1,replace=Los5)
	V = SubArray(V,element=2,replace=Los6)
	rectangular = 1
END
ELSE: BEGIN
	message, /info, 'Synax: elo = sphere_distance(Los1,Los2 [,/degrees,/rectangular])'
	message, /info,	'Synax: elo = sphere_distance(Los1p,Los1d,Los2p,Los2d [,/degrees])'
	message,		'Synax: elo = sphere_distance(Los1x,Los1y,Los1z,Los2x,Los2y,Los2z)'
END 
ENDCASE

SyncArgs, L, V

sza = size(L)

n = sza[sza[0]+2]/sza[1]				; # vector pairs

L = reform(L, [sza[1], n])
V = reform(V, [sza[1], n])

CASE rectangular OF
0: BEGIN
	L = L[0:1,*]*rpm
	V = V[0:1,*]*rpm

	phi = V[0,*]-L[0,*]

	IF arg_present(eps) THEN	$
		eps = atan( sin(phi)*cos(V[1,*]) , sin(V[1,*])*cos(L[1,*])-cos(V[1,*])*sin(L[1,*])*cos(phi) )/rpm

	L = sin(L[1,*])*sin(V[1,*])+cos(L[1,*])*cos(V[1,*])*cos(phi)

END
1: BEGIN
	L = L/([1,1,1]#sqrt(total(L*L,1)))	; Unit vector along line of sight
	V = V/([1,1,1]#sqrt(total(V*V,1)))	; Unit vector along line of sight

	L = total(L*V,1)					; Inner product = cos(angular distance)
END
ENDCASE

L = (L < 1) > (-1)
L = acos(L)/rpm

newsz = [sza[0]-1,sza[2:sza[0]+1],n]
SyncDims, L, size=newsz
IF IsType(eps,/defined) THEN SyncDims, eps, size=newsz

RETURN, L  &  END