[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

# astro coords pro: equatorial to horizon (rd2azel.pro)

• Subject: astro coords pro: equatorial to horizon (rd2azel.pro)
• From: John Keck <jwk(at)phys.columbia.edu>
• Date: Thu, 02 Dec 1999 15:17:27 -0500
• Newsgroups: comp.lang.idl-pvwave
• Organization: Columbia Astrophysics Laboratory
• Xref: news.doit.wisc.edu comp.lang.idl-pvwave:17594

Hi,

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
```