;+
; NAME:
;	SPHERE_CIRCLE
; PURPOSE:
;	Calculate spherical coordinates of points at fixed angular distance
;	from given center
; CATEGORY:
;	Spherical geometry
; CALLING SEQUENCE:
;	SPHERE_CIRCLE,latitude=LAT,elongation=ELONG,DLON,DLAT
; INPUTS:
;	LAT		latitude of center (degrees in [-90,90])
;	ELONG		angular radius (degrees)
;	AMIN		minimum value for angle A (see PROCEDURE)
;	AMAX		maximum value for angle A
; OUTPUTS:
;	DLON(180)	relative longitude of points on circle (degrees)
;	DLAT(180)	latitude of points on circle
; RESTRICTIONS:
;	Latitude must be in range [-90,90]
;	Latitude must be in range (0,180)
; PROCEDURE:
;    >	The array DLON returns the relative longitude, i.e. the difference
;	between between longitude of a point on the circumference of the circle
;	and the longitude of the center.
;    >	The calculation is performed in spherical triangle with corners at
;	the north pole P, the centre of the circle C and a point on the 
;	circumference X. ELONG is arc(CX), LAT is 90-arc(PC). Required
;	are DLAT=90-arc(PX) and DLON=angle(CPX).
;    >	AMAX-AMIN+1 points on the circle are calculated by varying the angle
;	A=angle(PCX) from AMIN to AMAX degrees, increasing in the clockwise 
;	direction from the CP meridian. Results are stored in DLON and DLAT.
; MODIFICATION HISTORY:
;	JUN-1992, Paul Hick (UCSD)
 pro SPHERE_CIRCLE,latitude=LAT,elongation=ELONG,DLON,DLAT, $
	angmin=AMIN,angmax=AMAX,plot=APLOT
;-

if n_elements(LAT)   eq 0 then LAT = 30.	; Check input
if abs(LAT) gt 90 then message, 'Invalid latitude'

if n_elements(ELONG) eq 0 then ELONG = 90.
if ELONG le 0 or ELONG ge 180 then message, 'Invalid elongation'

if not keyword_set(AMIN) then AMIN = 0
if not keyword_set(AMAX) then AMAX = 180
AMIN = AMIN mod 360  &  AMAX = AMAX mod 360
if AMIN gt AMAX then begin A = AMIN  &  AMIN = AMAX  &  AMAX = A  &  endif

if not keyword_set(APLOT) then APLOT = 0

RADEG = asin(1.d0)/90.d0			; Degree to radians

COLAT = 90-LAT					; Convert to colatitude
GAM = ELONG*RADEG  &  SGAM = sin(GAM)  &  CGAM = cos(GAM) ; Elongation
DEL = COLAT*RADEG  &  SDEL = sin(DEL)  &  CDEL = cos(DEL) ; Colatitude

A = (AMIN+indgen(AMAX-AMIN+1))*RADEG  &  CA = cos(A)  &  SA = sin(A)

DLAT = (CDEL*CGAM+SDEL*SGAM*CA) > (-1) < 1	; Cosine rule
DLAT = acos(DLAT)  &  DLAT = DLAT/RADEG		; Colatitude
DLAT = 90-DLAT					; Latitude

SA = SGAM*SA					; Sine rule
CA = SDEL*CGAM-CDEL*SGAM*CA			; Analogue ruleC
XA = where(SA ne 0 or CA ne 0)			; Exclude poles
DLON = fltarr(AMAX-AMIN+1)  &  DLON(XA) = atan(SA(XA),CA(XA))
DLON = DLON/RADEG				; Relative longitude

CA = 0  &  SA = 0  &  XA = 0			; Destroy arrays

if ELONG eq     COLAT and (AMIN le    0 or AMAX ge 360) then DLON(0)   = 90 ; Set north pole
if ELONG eq 180-COLAT and (AMIN le -180 or AMAX ge 180) then DLON(180) = 90 ; Set south pole

if APLOT then plot, DLON,DLAT

return  &  end