A while back (more than a year ago, maybe two), I posted a request for a program to convert (RA,Dec) to (azimuth, elevation). Some people replied with helpful suggestions and at least one asked for a copy of anything I came up with (unfortunately I lost his address). No one had what I needed, so I ended up writing my own program. The only external modules it requires are in the standard IDL astro library.
I have finally worked out the bugs and written up documentation and I present it here because I suspect that it may be useful to more than a few people out there.
Cheers,
John Keck
PRO de2h, ha, dec, phi, az, el
;+
; PURPOSE
; Equatorial to horizon coordinates: HA,Dec to Az,El
; INPUT
; ha hour angle, radians
; dec declination, radians
; phi latitude of observer, radians
;
; OUTPUT
; az azimuth, radians
; el elevation, radians
;
; *** ALL ARGUMENTS IN RADIANS ***
;
; EXAMPLE
; DE2H, ha, dec, phi, az, el
;
; adapted from a C version of slaDe2h, John Keck, May 1998
;-
; Useful trig functions
sh = SIN ( ha );
ch = COS ( ha );
sd = SIN ( dec );
cd = COS ( dec );
sp = SIN ( phi );
cp = COS ( phi );
; Az,El as x,y,z
x = - ch * cd * sp + sd * cp;
y = - sh * cd;
z = ch * cd * cp + sd * sp;
; To spherical */
r = SQRT ( x * x + y * y );
; a = ( r == 0.0 ) ? 0.0 : atan2 ( y, x ) ;
IF r EQ 0 THEN az = 0. ELSE az = ATAN2(y,x)
; *az = ( a < 0.0 ) ? a + D2PI : a;
IF az LT 0 THEN az = az + 2*!pi
; *el = atan2 ( z, r );
el = ATAN2 ( z, r )
RETURN
END
;----------------------------------------------------------------------
PRO rd2azel, ra,dec,u_time,lon,lat,az,el $
, DEBUG= db
;+
; RD2AZEL
; convery RA and Dec to Azimuth and Elevation
; (equatorial to horizon coordinates)
; INPUT
; ra decimal hours
; dec decimal degrees of arc
; u_time time structure with tags:
; (year, month, day, hour, minute, sec)
; lon degrees E
; lat degrees N
;
; OUTPUT
; az azimuth, degrees
; el elevation, degrees
;
; ***** ALL ARGUMENTS IN DEGREES, *****
; ** EXCEPT R.A., WHICH IS IN HOURS ***
; ***** (u_time has its own units) ****
;
; all arguments are scalars in this version
;
; EXAMPLE
; lat = -23.8 & lon = 133.40 ; Alice Springs, NT
; ra = 18.020083 & dec = -25.743333 ; GRS 1758-258 for epoch 2000
; u_time = {UT, year: 1995 $
; , Month: 10 $
; , day: 17 $
; , hour: 5, minute: 43, sec: 0.0}
;
; RD2AZEL,ra,dec,u_time,lon,lat,az,el
; PRINT, az, el
; 100.147 66.4507
;
; REQUISITES
; DE2H (included in this module)
; standard IDL astronomy library
; (viz. JULDATE, JDCNV, CT2LST)
;
; TBD
; make it accept vector arguments
;
; HISTORY
; adapted from a C program:
; John Keck (jwk@phys.columbia.edu),
; Columbia Astrophysics Lab, May 1998
;
;-
IF N_ELEMENTS(db) LT 1 THEN db = 0
; c60 = 1./[1,60,3600]
ut = TOTAL([u_time.hour,u_time.minute,u_time.sec]/[1.,60.,3600.])
; ut = REFORM([[u_time.hour],[u_time.minute],[u_time.sec]]#c60)
JULDATE,[u_time.year,u_time.month,u_time.day],rjd2
JDCNV, u_time.year,u_time.month,u_time.day,0,jd
CT2LST, gmst $ ; output
, 0 $ ; longitude for Greenwich
, 0 $ ; time zone
, jd+ut/24. ; Julian date
; gmst2 = LMST(jd, ut, 0)
gmst = gmst *(15*!dtor)
IF db THEN HELP,jd,gmst
eqeqx_corr = 0 ; approx.
hour_angle=(gmst + lon*!dtor + eqeqx_corr - ra*15*!dtor);
IF db THEN HELP,hour_angle
DE2H,hour_angle,dec*!dtor,lat*!dtor,az,el
az = az*!radeg
el = el*!radeg
RETURN
END