;+
; NAME:
;	sphere_great_arc
; PURPOSE:
;	Returns array of points along great circle
;	connecting two points on sphere
; CATEGORY:
;	sat/widget/qimage_cw
; CALLING SEQUENCE:
	FUNCTION sphere_great_arc, los1, los2, $
		npoints = npoints	, $
		degrees = degrees
; INPUTS:
;	los1
;	los2
; OPTIONAL INPUT PARAMETERS:
;	npoints = npoints
;	/degrees
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
;	InitVar, ToRadians, SyncArgs, sphere_distance, gridgen
; PROCEDURE:
; MODIFICATION HISTORY:
;	JUL-2008, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;-

InitVar, npoints, 101
rpm = ToRadians(degrees=degrees)

v1 = los1
v2 = los2

SyncArgs, v1, v2

sz = size(v1)
nlos = sz[sz[0]+2]/sz[1]			; sz[1] = 2 or 3

v1 = reform(v1, sz[1], nlos, /overwrite)
v2 = reform(v2, sz[1], nlos, /overwrite)

v1 = v1[0:1,*]*rpm					; [2,nlos]
v2 = v2[0:1,*]*rpm					; [2,nlos]

eta = sphere_distance(v1,v2,eps=eps); [nlos]
one = replicate(1.0d0,npoints)		; [npoints]

eta = gridgen(npoints,/one)#eta 	; [npoints,nlos]
eps = one#eps						; [npoints,nlos]
phi = one#reform(v1[0,*])			; [npoints,nlos]
lat = one#reform(v1[1,*])			; [npoints,nlos]

rr = dblarr(2,npoints*nlos,/nozero)
rr[0,*] = phi+atan(sin(eps)*sin(eta),cos(eta)*cos(lat)-sin(eta)*sin(lat)*cos(eps)) 
rr[1,*] = !dpi/2.0d0-acos( cos(eta)*sin(lat)+sin(eta)*cos(lat)*cos(eps) )

CASE sz[0] EQ 1 OF
0: rr = reform(rr,[2,npoints,sz[2:sz[0]]],/overwrite)
1: rr = reform(rr,/overwrite)
ENDCASE

RETURN, rr/rpm  &  END