;
;+
; NAME:
;	MOONAZEL
; PURPOSE:
;	Calculation of the direction to the Moon at a given time and
;	place.
; CATEGORY:
;
; CALLING SEQUENCE:
;	MOONAZEL, YEAR, MONTH, DAY, HOUR, MINUTE, SECOND $
;		, AZIMUTH, ELEVATION [, MRA, MDEC, PHASE ] $
;	       [, /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:
;	None.
; KEYWORD PARAMETERS:
;       /REFRACTION: Set this keyword to apply a refraction correction
;                    near the horizon.
;	LOC:    If supplied, this parameter must be a 2-element
;		array giving the [N latitude, E longitude] of the
;		place for which the Moon's direction is to be
;		computed.
; OUTPUTS:
;       AZIMUTH: The angle East of geographic North of the lunar
;		meridian
;       ELEVATION: The angle up from the horizon to the Moon in the
;		meridian specified by AZIMUTH (ELEVATION < 90 deg)
; OPTIONAL OUTPUTS:
;       MRA:    The lunar right ascension
;       MDEC:   The lunar declination
;	PHASE:	The geocentric angle between the Sun and the Moon.
; COMMON BLOCKS:
;	None.
; SIDE EFFECTS:
;	None.
; RESTRICTIONS:
;       The algorithm giving the lunar position is valid for an
;       uncertain period of time.  The accuracy of the position is
;       typically better than 0.3 deg in RA, and 0.2 deg in
;       declination.  The place for which the lunar 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 Moon in inertial
;       coordinates is calculated using a routine coded by Michael
;       R. Greason, Hughes STX Corporation, October 1988.  The
;	rest of the calculation is based on equations in the
;	Observer's Handbook.
; EXAMPLE:
;       moonazel,94,10,12,0,0,0,az,el,/ref      ; calculate lunar
;                                               ; az and el for
;                                               ; October 12, 1994,
;                                               ; 00:00 UT (Eureka).
; SEE ALSO:
;	moonpos.pro
; MODIFICATION HISTORY:
; 	Written by:	D P Steele, ISR, U of Calgary, March 1994.
;       Documented:     September 1994.
;-
PRO moonazel,year,month,day,hour,minute,second,azimuth,elevation $
	    ,mra,mdec,phase $
	    ,refraction=refraction,loc=loc

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

fhour=hour+(second/60.+minute)/60.

yr=year
l1900=WHERE(yr LT 1900.,nl1900)
IF nl1900 GT 0 THEN BEGIN
    l50=WHERE(yr(l1900) LT 50.,nl50)
    IF nl50 GT 0 THEN yr(l1900(l50))=yr(l1900(l50))+2000.
    IF nl50 LT nl1900 THEN BEGIN
	g50=setminus(INDGEN(N_ELEMENTS(l1900)),l50)
        yr(l1900(g50))=yr(l1900(g50))+1900.
	ENDIF
    ENDIF

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

jdf=ymdh2jdf(yr,month,day,fhour)
moonpos,jdf,mra,mdec,mx,my,mz	; mra, mdec in radians; mx, my, mz in Re

t2000=(DOUBLE(jdf)-2451545.0D0)/36525.0D0

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

xp=mx-COS(lat*!DTOR)*COS(locmeanst*!DTOR)		; topocentric
yp=my-COS(lat*!DTOR)*SIN(locmeanst*!DTOR)		; EI coordinates
zp=mz-SIN(lat*!DTOR)					; wrt lat, lon

rp=SQRT(xp*xp+yp*yp+zp*zp)
alphap=ATAN(yp,xp)					; apparent RA
wneg=WHERE(alphap LT 0,nwneg)
IF nwneg GT 0 THEN alphap(wneg)=alphap(wneg)+2.0*!DPI
deltap=ASIN(zp/rp)					; apparent dec

mra=alphap
mdec=deltap

hourangle=(locmeanst-!RADEG*mra) MOD 360.0D0		; ditto for hourangle
whan=WHERE(hourangle LT 0.0D0,nwhan)
IF nwhan GT 0 THEN hourangle(whan)=hourangle(whan)+360.0D0

elevation=!RADEG*ASIN(SIN(mdec)*SIN(!DTOR*lat) $
	             +COS(!DTOR*hourangle)*COS(mdec)*COS(!DTOR*lat))
IF KEYWORD_SET(refraction) THEN elevation=refract(elevation)

azimuth=!RADEG*ASIN(-COS(mdec)*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 N_PARAMS() EQ 11 THEN BEGIN	; moon phase is desired
    sunpos,year,month,day,fhour*3600.,gmst,eclong,sra,sdec
    sra=sra*!DTOR		; convert solar RA and Dec
    sdec=sdec*!DTOR		; to radians
    radiff=(mra+2*!DPI-sra) MOD (2*!DPI)
    ; phase below is geocentric angle between Moon and Sun
    phase=ACOS(SIN(sdec)*SIN(mdec)+COS(sdec)*COS(mdec)*COS(sra-mra))*!RADEG
    pastfull=(radiff GT !DPI)
    ; moon phase continues to increase past full moon
    phase=360.D0*pastfull+(1-2*pastfull)*phase
    ENDIF ELSE phase=-1.D0

RETURN
END