co_refract.pro

basic median filter simulation

;+
; 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: [e-mail]
;
; 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
;    Fixed problem with vector input when /TO_OBSERVED set W. Landsman Dec 2005
;-
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=P[i],T=T[i])
            cur = a[i] + dr ; guess of observed location

                repeat begin
                  last = cur
                  dr = co_refract_forward(cur,P=P[i],T=T[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