FUNCTION PA_Pole, T, lngsun=LngSun, equatorial=Equatorial, southpole=SouthPole, degrees=degrees
;+
; NAME:
;	PA_Pole
; PURPOSE:
;	Calculate the position angle of the solar north pole
; CATEGORY:
;	Astronomy: celestial physics
; CALLING SEQUENCE:
;	PA = PA_POLE(T, lngsun=LngSun)
; INPUTS:
;	T		array; type: time structure
;				times (UT)
;	lngsun=LngSun
;			array; type: float
;				ecliptic longitude of Sun as viewed from observer
; OPTIONAL INPUTS:
;	/equatorial		if set, the return value is relative to equatorial
;					north (default is ecliptic north)
;	/southpole		if set, return values for solar south pole
;					(default: solar north pole)
;	/degrees		if set, all angles are in degrees (default is radians)
; OUTPUTS:
;	PA_POLE	array; type: float
;				position angle
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
;	InitVar, ToRadians, TimeGet, TimeUnit
; PROCEDURE:
;	For given time T and ecliptic longitude of the Sun (as seen
;	from an observer in the ecliptic plane) the position angle
;	of the solar pole is calculated, either relative to ecliptic or to
;	equatiorial North.
;	See Spherical Astronomy pp. 430-433 by Green, Robin M., Cambridge UP
;	(1985) (The position angle is measured counterclockwise, i.e. toward
;	the east, from North)
; MODIFICATION HISTORY:
;	JUL-1992, Tom Davidson (UCSD/CASS)
;	SEP-1992, Paul Hick (UCSD/CASS), converted from F77 to IDL
;	FEB-1998, Paul Hick (UCSD/CASSu)
;		added keywords /degrees and /southpole
;	NOV-1999; Paul Hick (UCSD/CASS, pphick@ucsd.edu); switched to using
;		time structures instead of yr,doy
;-

InitVar, Equatorial, /key
InitVar, SouthPole , /key

rpd = ToRadians(/degrees)	; Conversion degree to radians

I =  7.25*rpd				; Inclination of Solar Eq. to ecliptic
E = 23.44*rpd				; Obliquity of the ecliptic

Yr  = TimeGet(T, timeunit(/year), /scalar)
Doy = TimeGet(T, timeunit(/day ), /fraction, /scalar)

O = Yr-1850+ (Doy-1)/365.25	; # years since 1850
O = (73.+2./3. + .5025/36.*O)*rpd
							; Longitude ascending node of solar equator on ecliptic
rpd = ToRadians(degrees=degrees)

PA = -atan(tan(I)*cos(O-LngSun*rpd))
IF Equatorial THEN PA = PA-atan(cos(LngSun*rpd)*tan(E))
IF SouthPole  THEN PA = PA+!pi

RETURN, PA/rpd  &  END