;------------------------------------------------------------------------------
;+
; NAME:
;       AZEL2RADEC
;
; PURPOSE:
;       To convert from horizon coordinates (azimuth-elevation)
;	to celestial coordinates (right ascension-declination)
;
; CALLING SEQUENCE:
;	AZEL2RADEC, az, el, ra, dec, year, month, day, time, lat, lng
;	
; INPUTS:
;	az:	azimuth (degrees)
;	el:	elevation (degrees)
;	year:	year (e.g. 1998)
;	month:  month, in numerical format (1-12, e.g. 8 for August)
;	day:	day of the month (1-31)
;	time:	Universal Time in decimal hours from 00:00:00 
;		(e.g. 4.5 = 04h 30m 00s AM)
;	lat:	geographical latitude (degrees)
;	lng:	geographical longitude (degrees) measured towards East
;
; OUTPUTS:
;       ra:	Right Ascension (hours)
;	dec:	declination (degrees)
;
; NOTES: 
;	The azimuth of the North Celestial Pole is 180 degrees
; 
; EXAMPLE:
;	Find the RA-dec coordinates corresponding to az=140 degs, el=40 degs
;	on 1998 August 2, 04h 00s 00s UT, in Palestine (Texas), lat=31.8 degs,
;	lng=264.3 degs E (or -95.7 degs E)
;
;	IDL> azel2radec,140.,40.,ra,dec,1998,8,2,4.,31.8,264.3
;	IDL> print,ra,dec
;               14.037406       56.872673
;
; PROCEDURES USED:
;       CT2LST - Convert from Civil Time to Local Sidereal Time
;
; MODIFICATION HISTORY:
;	Created, Amedeo Balbi, August 1998 (based partly on material
;	by Pedro Gil Ferreira)	
;       Modified, Amedeo Balbi, October 1998, to accept vectors as input
;-
;------------------------------------------------------------------------------
pro AZEL2RADEC,az,el,ra,dec,year,month,day,time,lat,lng

; converting from Universal Time to Greenwich Sidereal Time - this is done
; by calling CT2LST with longitude=0 and time zone=0 
; (the reason why I don't use CT2LST to obtain Local Sidereal Time directly
; is because I want to work with UT as input in the main procedure and I don't
; want to use the time zone as an additional input parameter)
CT2LST, gst, 0., 0., time, day, month, year

; converting from Greenwich Sidereal Time to Local Sidereal Time
; the formula is LST = GST + long
lst = float(gst)+lng/15. ; the long. has a "+" because is measured towards East

; converting from degrees to radians
rlat = lat*!dtor
raz = az*!dtor
rel = el*!dtor

; working out declination
rdec = asin( sin(rel)*sin(rlat)-cos(rel)*cos(raz)*cos(rlat) )
dec = rdec*!radeg

; working out right ascension
ha = atan(cos(rel)*sin(raz), sin(rel)*cos(rlat)+cos(rel)*cos(raz)*sin(rlat))
ha = 24.*ha/2./!pi
ra = lst-ha

index = where(ra lt 0.0, count)
if count gt 0 then ra(index) = 24.0 + (ra(index) mod 24)
ra = ra mod 24.0

return
end