;+
; NAME:
;	SUNAZEL
; PURPOSE:
;	Calculation of solar direction at a given place and time.
; CATEGORY:
;
; CALLING SEQUENCE:
;	SUNAZEL, YEAR, MONTH, DAY, HOUR, MINUTE, SECOND $
;	       , AZIMUTH, ELEVATION, SRA, SDEC [,/REFRACTION] $
;             [, LOC = LOC]
; INPUTS:
;	YEAR:	year, in either 4-digit or 2-digit representation
;	MONTH:	month, a number from 1 to 12
;	DAY:	day of month
;	HOUR:	hour, in Universal Time
;	MINUTE:	minute, in UT
;	SECOND:	second, in UT
; OPTIONAL INPUTS:
;       LOC:    If supplied, this parameter must be a 2-element
;               array giving the [N latitude, E longitude] of the
;               place for which the Sun's direction is to be
;		computed.
; KEYWORD PARAMETERS:
;       /REFRACTION: Set this keyword to apply a refraction correction
;		     near the horizon.
; OUTPUTS:
;       AZIMUTH: The angle East of geographic North of the solar
;		 meridian
;       ELEVATION: The angle up from the horizon to the Sun in the
;                  meridian specified by AZIMUTH (ELEVATION < 90 deg)
; OPTIONAL OUTPUTS:
;	SRA:	The solar right ascension
;	SDEC:	The solar declination
; COMMON BLOCKS:
;	None.
; SIDE EFFECTS:
;	None.
; RESTRICTIONS:
;       The algorithm giving the solar position is valid from 1901
;       through 2099 AD.  The place for which the solar direction
;       is calculated is assumed to be the location of the Polar
;	Camera, and this is inferred from the specified date.
; PROCEDURE:
;       The location of the Polar Camera is established from the
;       date specified.  The direction to the Sun in inertial
;       coordinates is calculated using an algorithm published by
;       C. T. Russell.  The rest of the calculation is based on
;	equations in the Observer's Handbook.
; EXAMPLE:
;	sunazel,94,10,12,0,0,0,az,el,/ref	; calculate solar
;						; az and el for 
;						; October 12, 1994,
;						; 00:00 UT (Eureka).
; SEE ALSO:
;	sunpos.pro
; MODIFICATION HISTORY:
; 	Written by:	D P Steele, ISR, U of Calgary, March 1994.
;	Documented:	September 1994.
;-
PRO sunazel,year,month,day,hour,minute,second,azimuth,elevation $
	   ,sra,sdec                                            $
	   ,refraction=refraction,loc=loc

; This routine computes the APPROXIMATE solar azimuth and elevation
; for a given observing location.  If enough arguments are specified,
; the solar right ascension and declination returned by sunpos.pro
; are also returned.

IF N_PARAMS() LT 7 THEN BEGIN
    doc_library,'sunazel'
    RETURN
    ENDIF

utsecs=(hour*60.+minute)*60.+second

PCsite=PoCasite(year,month,day)
ns=N_ELEMENTS(PCsite)

IF N_ELEMENTS(loc) EQ 0 THEN BEGIN

    ; No coordinates supplied, assume Polar Camera location

    lat=REPLICATE(51.08D0,ns)       ; default to Calgary
    lon=REPLICATE(245.8681D0,ns)

    rabb=WHERE(PCsite EQ 2,nrabb)   ; Rabbit Lake coordinates
    IF nrabb GT 0 THEN BEGIN
        lat(rabb)=58.23D0
        lon(rabb)=256.333D0
        ENDIF
    eur=WHERE(PCsite EQ 3,neur)     ; Eureka coordinates
    IF neur GT 0 THEN BEGIN
        lat(eur)=80.0533D0
        lon(eur)=273.5839D0
        ENDIF

    other=WHERE(ABS(PCsite-2) GT 1,nother)
    IF nother GT 0 THEN BEGIN
        PRINT,'An unrecognized Polar Camera site is indicated!'
        PRINT,'Enter the N lat and E long (deg)'
        READ,olat,olon
        lat(other)=olat
        lon(other)=olon
        END

    ENDIF ELSE BEGIN

    ; Use the supplied coordinates

    lat=REPLICATE(loc(0),ns)
    lon=REPLICATE(loc(1),ns)
    ENDELSE

IF ns EQ 1 THEN BEGIN           ; scalar inputs, scalar outputs
    lat=lat(0)
    lon=lon(0)
    ENDIF

sunpos,year,month,day,utsecs,gmst,eclong,sra,sdec	; sra, sdec in degrees
sra=sra*!DTOR						; convert to radians
sdec=sdec*!DTOR

jdf=ymdh2jdf(year,month,day,utsecs/3600.)
t2000=(DOUBLE(jdf)-2451545.0D0)/36525.0D0

locmeanst=lmst(ROUND(jdf),utsecs/86400.,lon)*360.D0	; lmst in degrees

hourangle=(locmeanst-!RADEG*sra+360.D0) MOD 360.0D0	; hourangle in degrees

elevation=!RADEG*ASIN(SIN(sdec)*SIN(!DTOR*lat) $
	             +COS(!DTOR*hourangle)*COS(sdec)*COS(!DTOR*lat))

azimuth=!RADEG*ASIN(-COS(sdec)*SIN(!DTOR*hourangle)/COS(!DTOR*elevation))

ca=WHERE((((0.0D0 LE azimuth) AND (azimuth LT 90.0D0)) AND $
          ((270.0D0 LE hourangle) AND (hourangle LT 360.0D0))) OR $
         (((-90.0D0 LE azimuth) AND (azimuth LT 0.0D0)) AND $
          ((0.0D0 LE hourangle) AND (hourangle LT 90.0D0))),nca)
IF nca GT 0 THEN azimuth(ca)=180.0D0-azimuth(ca)

lz=WHERE(azimuth LT 0.0D0,nlz)
IF nlz GT 0 THEN azimuth(lz)=azimuth(lz)+360.0D0

IF KEYWORD_SET(refraction) THEN elevation=refract(elevation)

RETURN
END