FUNCTION refract,true_altitude $
		,pressure=pr,temperature=tr,radians=rad,debug=dbg

;/* correct the true altitude, ta, for refraction to the apparent altitude, aa,
; * each in radians, given the local atmospheric pressure, pr, in mbars, and
; * the temperature, tr, in degrees C. 
; */
; Angles are assumed to be in degrees unless the radians keyword is set.

;void refract ( double pr, double tr, double ta, double *aa )
;{
;	double r;	/* refraction correction*/

sz=SIZE(true_altitude)
scalar=(sz(0) EQ 0)

IF scalar THEN true_altitude=FLTARR(1)+true_altitude

IF KEYWORD_SET(rad) THEN ta=true_altitude ELSE ta=true_altitude*!DTOR

a=FLTARR(N_ELEMENTS(ta))
b=a
r=a

IF NOT KEYWORD_SET(pr) THEN pr=1000.	; default pressure
IF NOT KEYWORD_SET(tr) THEN tr=-25.	; default temperature

; There are three cases to handle: ta < -5 deg, -5 deg < ta < 15 deg, 15 deg < ta

; if (ta >= degrad(15.)) {
g15=WHERE(15.*!DTOR LE ta,ng15)
; /* model for altitudes at least 15 degrees above horizon */
IF ng15 GT 0 THEN r(g15)=7.888888E-5*pr/((273+tr)*TAN(ta(g15)))

; altitudes near the horizon
mid=WHERE((-5.*!DTOR LE ta) AND (ta LT 15.*!DTOR),nmid)
; /* hairier model for altitudes at least -5 and below 15 degrees */
IF nmid GT 0 THEN BEGIN
    tad=ta(mid)*!RADEG
    a(mid)=((2E-5*tad+1.96E-2)*tad+1.594E-1)*pr
    b(mid)=(273+tr)*((8.45E-2*tad+5.05E-1)*tad+1)
    r(mid)=!DTOR*a(mid)/b(mid)
    ENDIF

; We don't have to do anything for altitudes more than 5 deg below
; the horizon

IF KEYWORD_SET(rad) $
    THEN apparent_altitude=ta+r $
    ELSE apparent_altitude=(ta+r)*!RADEG

IF scalar THEN apparent_altitude=apparent_altitude(0)

IF KEYWORD_SET(dbg) THEN BEGIN
    HELP
    STOP
    ENDIF

RETURN,apparent_altitude
END