;+
; procedure to calculate the azimuths and zenith angles
; of stars from a list of stars and their right ascensions and
; declinations.
;
; Calling sequence:
;       PRO STEPHEM,yr,mo,day,hr,min,sec,lat,lon,nstars,desig,az,za,mag
;
; Inputs:
;       yr,mo,day,hr,min,sec,lat,lon :
;               Universal Time of exposure of image, and N. latitude  
;               and E. longitude of camera.
;               Note that YR must have 4 digits.
; Outputs:
;       nstars: number of stars in the list which are visible above the
;               horizon at the given time
;       desig : array of designations of these stars
;       az :    array of azimuths of the visible stars
;       za :    array of zenith angles of visible stars
;       mag :   array of magnitudes of visible stars
;
; Written by David van der Schee, U. Calgary, May 1993.
;-

PRO STEPHEM,yr,mo,day,hr,min,sec,lat,lon,nstars,desig,az,za,mag
ut = (sec/3600. + min/60. + hr)/24.     ; fraction of a day
jday = 0.
jday = JULDAY(mo, day, yr)
lsidt = LMST(jday, ut, lon)             ; local mean sidereal time, as
                                        ; fraction of a day
                                        ; load star catalog
rdybsc4, n,desig,rtasc,dec,mag;,path='[steele.idl]ybsc4s.dat'
rtasc = rtasc*!DTOR
dec = dec * !DTOR
lat = lat * !DTOR
lsidt = lsidt*2*!PI
hangle = lsidt - rtasc                  ; hour angle, radians
sinalt = SIN(dec)*SIN(lat) + COS(hangle)*COS(dec)*COS(lat)
za = ACOS(sinalt)
sinaz = -COS(dec)*SIN(hangle)/SIN(za)
cosaz = (SIN(dec) - sinalt*sin(lat))/(SIN(za)*COS(lat))
az = ATAN(sinaz, cosaz)                 ; azimuth, radians
vis = WHERE(sinalt GT 0.)               ; visible stars
za = za(vis)
mag=mag(vis)
az = az(vis)
desig = desig(vis)
nstars = N_ELEMENTS(vis)
rtasc = rtasc(vis)
dec=dec(vis)
az = az*!RADEG
za = za*!RADEG
RETURN
END