;+ ; NAME: ; HOR2EQ ; ; PURPOSE: ; Converts local horizon coords (alt-az) of something to equatorial (ra-dec). ; ; EXPLANATION: ; This is a nice code to calculate equatorial (ra,dec) coordinates from ; horizon (alt,az) coords. It is typically accurate to about 1 arcsecond ; or better (I have checked the output against the publicly available XEPHEM ; software). It preforms precession, nutation, abberation, and refraction ; corrections. The perhaps best thing about it is that it can take arrays ; as inputs, in all variables and keywords EXCEPT Lat, lon, and Altitude ; (the code assumes these aren't changing), and uses vector arithmetic in ; every calculation except when calculating the precession matrices. ; ; CALLING SEQUENCE: ; ; HOR2EQ, alt, az, jd, ra, dec, [ha, LAT= , LON= , /WS, OBSNAME= , $ ; /B1950 , PRECESS_= 0, NUTATE_= 0, REFRACT_= 0, $ ; ABERRATION_= 0, ALTITUDE= , /VERBOSE, _EXTRA= ] ; ; ; INPUT VARIABLES ; alt : altitude (in degrees) [scalar or vector] ; az : azimuth angle (in degrees, measured EAST from NORTH, but see ; keyword WS below.) [scalar or vector] ; JD : Julian Date [scalar or vector] ; Note: if RA and DEC are arrays, then alt and az will also be arrays. ; If RA and DEC are arrays, JD may be a scalar OR an array of ; the same dimensionality. ; ; OPTIONAL INPUT KEYWORDS: ; lat : north geodetic latitude of location in degrees ; lon : EAST longitude of location in degrees ; (Specify west longitude with a negative sign.) ; /WS : Set this to get the azimuth measured westward from south ; (not East of North). ; obsname : Set this to a valid observatory name to be used by the ; astrolib OBSERVATORY procedure, which will return the latitude ; and longitude to be used by this program. ; /B1950 : Set this if your ra and dec are specified in B1950, ; FK4 coordinates (instead of J2000, FK5) ; precess_ : Set this to 1 to force precession [default], 0 for no ; precession. ; nutate_ : Set this to 1 to force nutation [default], 0 for no nutation. ; aberration_ : Set this to 1 to force aberration correction [default], ; 0 for no correction. ; refract_ : Set to 1 to force refraction correction [default], 0 for ; no correction. ; altitude: The altitude of the observing location, in meters. [default=0]. ; /verbose: Set this for verbose output. The default is verbose=0. ; _extra: This is for setting TEMPERATURE or PRESSURE explicity, which are ; used by CO_REFRACT to calculate the refraction effect of the ; atmosphere. If you don't set these, the program will make an ; intelligent guess as to what they are (taking into account your ; altitude). See CO_REFRACT for more details. ; ; OUTPUT VARIABLES ; ra : Right Ascension of object (J2000) in degrees (FK5); scalar or ; vector. ; dec : Declination of object (J2000) in degrees (FK5), scalar or vector. ; ha : hour angle (in degrees) (optional) ; ; DEPENDENCIES: ; NUTATE, PRECESS, ADSTRING(), SUNPOS, OBSERVATORY (from the astrolib) ; CO_NUTATE, CO_ABERRATION, CO_REFRACT, HADEC2ALTAZ, VECTOR2RADEC ; ; BASIC STEPS ; Precess Ra-Dec to current equinox. ; Nutation Correction to Ra-Dec ; Aberration correction to Ra-Dec ; Calculate Local Mean Sidereal Time ; Calculate Local Apparent Sidereal Time ; Calculate Hour Angle ; Do Spherical Trig to find Apparent Alt-Az ; Apply refraction correction to find observed Alt. ; ;CORRECTIONS I DO NOT MAKE: ; * Deflection of Light by the sun due to GR. (typically milliarcseconds, ; can be arseconds within one degree of the sun) ; * The Effect of Annual Parallax (typically < 1 arcsecond) ; * and more (see below) ; ; TO DO ; * Better Refraction Correction. Need to put in wavelength dependence, ; and integrate through the atmosphere. ; * Topocentric Parallax Correction (will take into account elevation of ; the observatory) ; * Proper Motion (but this will require crazy lookup tables or something). ; * Difference between UTC and UT1 in determining LAST -- is this important? ; * Effect of Annual Parallax (is this the same as topocentric Parallax?) ; * Polar Motion ; * Better connection to Julian Date Calculator. ; ; EXAMPLE: ; ; You are at Kitt Peak National Observatory, looking at a star at azimuth ; angle 264d 55m 06s and elevation 37d 54m 41s (in the visible). Today is ; Dec 25, 2041 and the local time is 10 PM precisely. What is the ra and dec ; (J2000) of the star you're looking at? The temperature here is about 0 ; Celsius, and the pressure is 781 millibars. The Julian date for this ; time is 2466879.79083333 ; ; IDL> hor2eq, ten(37,54,41), ten(264,55,06), 2466879.7083333d, ra, dec, $ ; /verb, obs='kpno', pres=781.0, temp=273.0 ; ; The program produces this output (because the VERBOSE keyword was set): ; ; Latitude = +31 57 48.0 Longitude = *** 36 0.0 ; longitude prints weirdly b/c of negative input to ADSTRING!! ; Julian Date = 2466879.708333 ; Az, El = 17 39 40.4 +37 54 41.0 (Observer Coords) ; Az, El = 17 39 40.4 +37 53 39.6 (Apparent Coords) ; LMST = +03 53 54.1 ; LAST = +03 53 53.6 ; Hour Angle = +03 38 30.1 (hh:mm:ss) ; Ra, Dec: 00 15 23.5 +15 25 1.9 (Apparent Coords) ; Ra, Dec: 00 15 24.2 +15 25 0.1 (J2041.9841) ; Ra, Dec: 00 13 14.1 +15 11 0.3 (J2000) ; ; The star is therefore Algenib! Compare the derived Ra, Dec with what XEPHEM ; got: ; Ra, Dec: 00 13 14.2 +15 11 1.0 (J2000) ; ; AUTHOR: ; Chris O'Dell ; Univ. of Wisconsin-Madison ; Observational Cosmology Laboratory ; Email: odell@cmb.physics.wisc.edu ;- pro hor2eq, alt, az, jd, ra, dec, ha, lat=lat, lon=lon, WS=WS, obsname=obsname,$ B1950 = B1950, verbose=verbose, precess_=precess_, nutate_=nutate_, $ refract_ = refract_, aberration_ = aberration_, altitude=altitude, $ _extra = _extra if N_params() LT 4 then begin print,'Syntax - HOR2EQ, alt, az, jd, ra, dec, [ha, LAT= , LON= , /WS, ' print,' OBSNAME= ,/B1950 , PRECESS_= 0, NUTATE_= 0, REFRACT_= 0, ' print,' ABERRATION_= 0, ALTITUDE= , /VERBOSE, TEMPERATURE=, PRESSURE=' return endif ;******************************************************************************* ; INITIALIZE STUFF ; If no lat or lng entered, use Pine Bluff Observatory values if n_elements(lat) eq 0 then lat = 43.0783d ; (btw, this is the declination of the zenith) if n_elements(lon) eq 0 then lon = -89.865d if keyword_set(obsname) then begin ;override lat,lon if observatory name has been specified observatory, obsname, obs lat = obs.latitude lon = -1*obs.longitude ; minus sign is becase OBSERVATORY uses west ; ;longitude as positive. altitude = obs.altitude endif if n_elements(precess_) eq 0 then precess_ = 1 if n_elements(nutate_) eq 0 then nutate_ = 1 if n_elements(aberration_) eq 0 then aberration_ = 1 if n_elements(refract_) eq 0 then refract_ = 1 v = keyword_set(verbose) ; converstion factors d2r = !dpi/180. h2d = 15. alt_ = alt ;do this so we don't change ra, dec arrays. az_ = az if v then print, 'Latitude = ', adstring(lat), ' Longitude = ', adstring(lon) if v then print, 'Julian Date = ', jd, format='(A,f15.6)' if v then print,'Az, El = ', adstring(az_, alt_), ' (Observer Coords)' ;******************************************************************************* ; Make Correction for ATMOSPHERIC REFRACTION (visible only, NOT radio!) if refract_ then alt_ = co_refract(alt_, altitude=altitude, _extra=_extra) if v then print,'Az, El = ', adstring(az_, alt_), ' (Apparent Coords)' if keyword_set(WS) then az_ = az_ - 180. co_nutate, jd, 45.,45., dra1, ddec1, eps=eps, d_psi=d_psi ;****************************************************************************** ;Calculate LOCAL APPARENT SIDEREAL TIME ; first get local mean sidereal time (lmst) ; get LST (in hours) - note:this is indep of tzone since giving jd ct2lst, lmst, lon, 0, jd lmst = lmst*h2d ; convert LMST to degrees (btw, this is the RA of the zenith) ; calculate local APPARENT sidereal time (last) last = lmst + d_psi *cos(eps)/3600. ; add correction in degrees if v then print, 'LMST = ', adstring(lmst/15.) if v then print, 'LAST = ', adstring(last/15.) ;**************************************************************************** ; Now do the spheical trig to get APPARENT Hour Angle [degrees], and ; declination [degrees]. altaz2hadec, alt_, az_, lat, ha, dec ; Find Right Ascension (in degrees, from 0 to 360.) ra = (last - ha + 360.) mod 360. if v then print, 'Hour Angle = ', adstring(ha/15.), ' (hh:mm:ss)' if v then print, 'Ra, Dec: ', adstring(ra,dec), ' (Apparent Coords)' ;***************************************************************************** ; calculate NUTATION and ABERRATION Corrections to Ra-Dec co_nutate, jd, ra, dec, dra1, ddec1, eps=eps, d_psi=d_psi co_aberration, jd, ra, dec, dra2, ddec2, eps=eps ;****************************************************************************** ; Make Nutation and Aberration Corrections (if wanted) ra = ra - (dra1*nutate_ + dra2*aberration_)/3600. dec = dec - (ddec1*nutate_ + ddec2*aberration_)/3600. J_now = (JD - 2451545.)/365.25 + 2000.0 ; compute current equinox if v then print, 'Ra, Dec: ', adstring(ra,dec), ' (J'+ $ strcompress(string(J_now),/rem)+')' ;***************************************************************************** ; PRECESS coordinates to current date ; (uses astro lib procedure PRECESS.pro) if precess_ then begin if keyword_set(B1950) then begin for i=0,n_elements(jd)-1 do begin ra_i = ra[i] & dec_i = dec[i] precess, ra_i, dec_i, J_now[i], 1950.0, /FK4 ra[i] = ra_i & dec[i] = dec_i endfor endif else begin for i=0,n_elements(jd)-1 do begin ra_i = ra[i] & dec_i = dec[i] precess, ra_i, dec_i, J_now[i], 2000.0 ra[i] = ra_i & dec[i] = dec_i endfor endelse endif if keyword_set(B1950) then s_now=' (J1950)' else s_now=' (J2000)' if v then print, 'Ra, Dec: ', adstring(ra,dec), s_now Return END