;+
; NAME:
;	PACEPOLE
; PURPOSE:
;       Determine the geographic coordinates of the PACE geomagnetic
;	pole in the northern hemisphere.
; CATEGORY:
;	Coordinate transformations
; CALLING SEQUENCE:
;	PACEpole, HT, GLAT, GLONG
; INPUTS:
;       HT:	The height (km) for which the pole's coordinates
;		are desired.
; OPTIONAL INPUTS:
;	None.
; KEYWORD PARAMETERS:
;	None.
; OUTPUTS:
;       GLAT:   The geographic (geodetic) latitude of the PACE
;		pole, in degrees.
;       GLONG:  The geographic longitude of the PACE pole, in
;		degrees.
; OPTIONAL OUTPUTS:
;	None.
; COMMON BLOCKS:
;	None.
; SIDE EFFECTS:
;       A file (~cnsr3/pace/pace.pole) is read to get the pole
;	coordinates at 25 km altitude steps.
; RESTRICTIONS:
;       Northern hemisphere only.  PACE coordinate system only
;	defined below 600 km altitude.
; PROCEDURE:
;       If HT is tabulated in the lookup file, the tabulated GLAT
;       and GLONG are returned,, otherwise linear interpolation in
;       geocentric cartesian coordinates is used to estimate the
;	coordinates.
; EXAMPLE:
;	PACEpole,250.,glat,glong
; SEE ALSO:
;	CONVERT.PRO
; MODIFICATION HISTORY:
; 	Written by:	DPSteele, ISR, 950601
;-

PRO pacepole,ht,glat,glong

; Need help?

IF N_PARAMS() LT 3 THEN BEGIN
    doc_library,'pacepole'
    RETURN
    ENDIF

; Ensure a scalar altitude.

IF N_ELEMENTS(ht) GT 1 THEN BEGIN
    MESSAGE,'HT must be a scalar',/INFORMATIONAL
    glat=0
    glong=0
    RETURN
    ENDIF ELSE ht=ht(0)

; Ensure HT within range of validity.

IF (ht GT 600) OR (ht LT 0) THEN BEGIN
    MESSAGE,'PACE coordinates valid within 0..600 km only' $
	   ,/INFORMATIONAL
    glat=0
    glong=0
    RETURN
    ENDIF

; Read the lookup table.

rdcols,'~cnsr3/pace/pace.pole',n,alt,glatp,glongp	; for 950101

; Is there a table entry for the specified ALT?

wht=WHERE(alt EQ ht,nwht)

IF nwht EQ 1 THEN BEGIN	; The coordinates are already tabulated!

    glat=glatp(wht)
    glong=glongp(wht)

    ENDIF ELSE BEGIN	; Get them by interpolation

    polrec3d,relc(90.)+alt,90.-glatp,glongp,xp,yp,zp,/degrees
    xpole=interpol(xp,alt,ht)
    ypole=interpol(yp,alt,ht)
    zpole=interpol(zp,alt,ht)
    recpol3d,xpole,ypole,zpole,rpole,colatpole,glong,/degrees
    glat=90.-colatpole		; PACE pole at (glat,glong)

    ENDELSE

; Ensure -180 < glong < 180.

glong=glong-(glong GT 180.)*360.

glat=glat(0)
glong=glong(0)

RETURN
END