/usr/share/gnudatalanguage/astrolib/hadec2altaz.pro is in gdl-astrolib 2018.02.16+dfsg-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 | PRO hadec2altaz, ha, dec, lat, alt, az, WS=WS
;+
; NAME:
; HADEC2ALTAZ
; PURPOSE:
; Converts Hour Angle and Declination to Horizon (alt-az) coordinates.
; EXPLANATION:
; Can deal with NCP/SCP singularity. Intended mainly to be used by
; program EQ2HOR
;
; CALLING SEQUENCE:
; HADEC2ALTAZ, ha, dec, lat ,alt ,az [ /WS ]
;
; INPUTS
; ha - the local apparent hour angle, in DEGREES, scalar or vector
; dec - the local apparent declination, in DEGREES, scalar or vector
; lat - the local latitude, in DEGREES, scalar or vector
;
; OUTPUTS
; alt - the local apparent altitude, in DEGREES.
; az - the local apparent azimuth, in DEGREES, all results in double
; precision
; OPTIONAL KEYWORD INPUT:
; /WS - Set this keyword for the output azimuth to be measured West from
; South. The default is to measure azimuth East from North.
;
; EXAMPLE:
; What were the apparent altitude and azimuth of the sun when it transited
; the local meridian at Pine Bluff Observatory (Lat=+43.07833 degrees) on
; April 21, 2002? An object transits the local meridian at 0 hour angle.
; Assume this will happen at roughly 1 PM local time (18:00 UTC).
;
; IDL> jdcnv, 2002, 4, 21, 18., jd ; get rough Julian date to determine
; ;Sun ra, dec.
; IDL> sunpos, jd, ra, dec
; IDL> hadec2altaz, 0., dec, 43.078333, alt, az
;
; ===> Altitude alt = 58.90
; Azimuth az = 180.0
; REVISION HISTORY:
; Written Chris O'Dell Univ. of Wisconsin-Madison May 2002
;-
if N_params() LT 4 then begin
print,'Syntax - HADEC2ALTAZ, ha, dec, lat ,alt ,az [ /WS ]'
return
endif
d2r = !dpi/180.
sh = sin(ha*d2r) & ch = cos(ha*d2r)
sd = sin(dec*d2r) & cd = cos(dec*d2r)
sl = sin(lat*d2r) & cl = cos(lat*d2r)
x = - ch * cd * sl + sd * cl
y = - sh * cd
z = ch * cd * cl + sd * sl
r = sqrt(x^2 + y^2)
; now get Alt, Az
az = atan(y,x) /d2r
alt = atan(z,r) / d2r
; correct for negative AZ
w = where(az LT 0)
if w[0] ne -1 then az[w] = az[w] + 360.
; convert AZ to West from South, if desired
if keyword_set(WS) then az = (az + 180.) mod 360.
END
|