;+
; 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

; Set up OS-dependent parameters.

@isitdos

; 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,pgmroot+dd+'pace.pol',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