;+ ; NAME: ; CO_REFRACT() ; ; PURPOSE: ; Calculate correction to altitude due to atmospheric refraction. ; ; DESCRIPTION: ; CO_REFRACT can calculate both apparent altitude from observed altitude and ; vice-versa. ; ; CALLING SEQUENCE: ; new_alt = CO_REFRACT(old_alt, [ ALTITUDE= , PRESSURE= , $ ; TEMPERATURE= , /TO_OBSERVED , EPSILON= ]) ; ; INPUT: ; old_alt - Observed (apparent) altitude, in DEGREES. (apparent if keyword ; /TO_OBSERVED set). May be scalar or vector. ; ; OUTPUT: ; Function returns apparent (observed) altitude, in DEGREES. (observed if ; keyword /TO_OBSERVED set). Will be of same type as input ; altitude(s). ; ; OPTIONAL KEYWORD INPUTS: ; ALTITUDE : The height of the observing location, in meters. This is ; only used to determine an approximate temperature and pressure, ; if these are not specified separately. [default=0, i.e. sea level] ; PRESSURE : The pressure at the observing location, in millibars. ; TEMPERATURE: The temperature at the observing location, in Kelvin. ; EPSILON: When keyword /TO_OBSERVED has been set, this is the accuracy ; to obtain via the iteration, in arcseconds [default = 0.25 ; arcseconds]. ; /TO_OBSERVED: Set this keyword to go from Apparent->Observed altitude, ; using the iterative technique. ; ; Note, if altitude is set, but temperature or pressure are not, the ; program will make an intelligent guess for the temperature and pressure. ; ; DESCRIPTION: ; ; Because the index of refraction of air is not precisely 1.0, the atmosphere ; bends all incoming light, making a star or other celestial object appear at ; a slightly different altitude (or elevation) than it really is. It is ; important to understand the following definitions: ; ; Observed Altitude: The altitude that a star is SEEN to BE, with a telescope. ; This is where it appears in the sky. This is always ; GREATER than the apparent altitude. ; ; Apparent Altitude: The altitude that a star would be at, if *there were no ; atmosphere* (sometimes called "true" altitude). This is ; usually calculated from an object's celestial coordinates. ; Apparent altitude is always LOWER than the observed ; altitude. ; ; Thus, for example, the Sun's apparent altitude when you see it right on the ; horizon is actually -34 arcminutes. ; ; This program uses couple simple formulae to estimate the effect for most ; optical and radio wavelengths. Typically, you know your observed altitude ; (from an observation), and want the apparent altitude. To go the other way, ; this program uses an iterative approach. ; ; EXAMPLE: ; The lower limb of the Sun is observed to have altitude of 0d 30'. ; Calculate the the true (=apparent) altitude of the Sun's lower limb using ; mean conditions of air pressure and temperature ; ; IDL> print, co_refract(0.5) ===> 0.025degrees (1.55') ; WAVELENGTH DEPENDENCE: ; This correction is 0 at zenith, about 1 arcminute at 45 degrees, and 34 ; arcminutes at the horizon FOR OPTICAL WAVELENGTHS. The correction is ; NON-NEGLIGIBLE at all wavelengths, but is not very easily calculable. ; These formulae assume a wavelength of 550 nm, and will be accurate to ; about 4 arcseconds for all visible wavelengths, for elevations of 10 ; degrees and higher. Amazingly, they are also ACCURATE FOR RADIO ; FREQUENCIES LESS THAN ~ 100 GHz. ; ; It is important to understand that these formulae really can't do better ; than about 30 arcseconds of accuracy very close to the horizon, as ; variable atmospheric effects become very important. ; ; REFERENCES: ; 1. Meeus, Astronomical Algorithms, Chapter 15. ; 2. Explanatory Supplement to the Astronomical Almanac, 1992. ; 3. Methods of Experimental Physics, Vol 12 Part B, Astrophysics, ; Radio Telescopes, Chapter 2.5, "Refraction Effects in the Neutral ; Atmosphere", by R.K. Crane. ; ; ; DEPENDENCIES: ; CO_REFRACT_FORWARD (contained in this file and automatically compiled). ; ; AUTHOR: ; Chris O'Dell ; Univ. of Wisconsin-Madison ; Observational Cosmology Laboratory ; Email: odell@cmb.physics.wisc.edu ; ; REVISION HISTORY: ; version 1 (May 31, 2002) ; Update iteration formula, W. Landsman June 2002 ; Corrected slight bug associated with scalar vs. vector temperature and ; pressure inputs. 6/10/2002 ;- function co_refract_forward, a, P=P, T=T ; INPUTS ; a = The observed (apparent) altitude, in DEGREES. ; May be scalar or vector. ; ; INPUT KEYWORDS ; P: Pressure [in millibars]. Default is 1010 millibars. [scalar or vector] ; T: Ground Temp [in Celsius]. Default is 0 Celsius. [scalar or vector] d2r = !dpi/180. if n_elements(P) eq 0 then P = 1010. if n_elements(T) eq 0 then T = 283. ; you have observed the altitude a, and would like to know what the "apparent" ; altitude is (the one behind the atmosphere). w = where(a LT 15.) R = 0.0166667/tan((a + 7.31/(a+4.4))*d2r) ;R = 1.02/tan((a + 10.3/(a+5.11))*d2r)/60. ; this formula goes the other direction! if w[0] ne -1 then R[w] = 3.569*(0.1594 + .0196*a[w] + $ .00002*a[w]^2)/(1.+.505*a[w]+.0845*a[w]^2) tpcor = P/1010. * 283/T R = tpcor * R return, R END function co_refract, a, altitude=altitude, pressure=pressure, $ temperature=temperature, To_observed=To_observed, epsilon=epsilon ; This is the main window. Calls co_refract_forward either iteratively or a ; single time depending on the direction we are going for refraction. o = keyword_set(To_observed) alpha = 0.0065 ; temp lapse rate [deg C per meter] if n_elements(altitude) eq 0 then altitude = 0. if n_elements(temperature) eq 0 then begin if altitude GT 11000 then temperature = 211.5 $ else temperature = 283.0 - alpha*altitude endif ; estimate Pressure based on altitude, using U.S. Standard Atmosphere formula. if n_elements(pressure) eq 0 then $ pressure = 1010.*(1-6.5/288000*altitude)^5.255 if n_elements(epsilon) eq 0 then $ epsilon = 0.25 ; accuracy of iteration for observed=1 case, in arcseconds if not o then begin aout = a - co_refract_forward(a,P=pressure,T=temperature) endif else begin aout = a*0. na = n_elements(a) ; if there are multiple elevations but only one temp and pressure entered, ; expand those to be arrays of the same size. P = pressure + a*0. & T = temperature + a*0. for i=0,na-1 do begin ;calculate initial refraction guess dr = co_refract_forward(a[i],P=pressure[i],T=temperature[i]) cur = a[i] + dr ; guess of observed location repeat begin last = cur dr = co_refract_forward(cur,P=pressure[i],T=temperature[i]) cur= a[i] + dr endrep until abs(last-cur)*3600. LT epsilon aout[i] = cur endfor endelse if N_elements(aout) GT 1 then return, reform(aout) else return, aout END