source: Sophya/trunk/SophyaExt/XephemAstroLib/refract.c@ 3632

#include <stdio.h>
#include <math.h>

#include "astro.h"

static void unrefractLT15 (double pr, double tr, double aa, double *ta);
static void unrefractGE15 (double pr, double tr, double aa, double *ta);

void
unrefract (double pr, double tr, double aa, double *ta)
{
#define LTLIM   14.5
#define GELIM   15.5

        double aadeg = raddeg(aa);

        if (aadeg < LTLIM)
            unrefractLT15 (pr, tr, aa, ta);
        else if (aadeg >= GELIM)
            unrefractGE15 (pr, tr, aa, ta);
        else {
            /* smooth blend -- important for inverse */
            double taLT, taGE, p;

            unrefractLT15 (pr, tr, aa, &taLT);
            unrefractGE15 (pr, tr, aa, &taGE);
            p = (aadeg - LTLIM)/(GELIM - LTLIM);
            *ta = taLT + (taGE - taLT)*p;
        }
}

static void
unrefractGE15 (double pr, double tr, double aa, double *ta)
{
        double r;
        
        r = 7.888888e-5*pr/((273+tr)*tan(aa));
        *ta  =  aa - r;
}

static void
unrefractLT15 (double pr, double tr, double aa, double *ta)
{
        double aadeg = raddeg(aa);
        double r, a, b;

        a = ((2e-5*aadeg+1.96e-2)*aadeg+1.594e-1)*pr;
        b = (273+tr)*((8.45e-2*aadeg+5.05e-1)*aadeg+1);
        r = degrad(a/b);

        *ta  =  (aa < 0 && r < 0) ? aa : aa - r;        /* 0 below ~5 degs */
}

/* 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.
 */
void
refract (double pr, double tr, double ta, double *aa)
{
#define MAXRERR degrad(0.1/3600.)       /* desired accuracy, rads */

        double d, t, t0, a;

        /* first guess of error is to go backwards.
         * make use that we know delta-apparent is always < delta-true.
         */
        unrefract (pr, tr, ta, &t);
        d = 0.8*(ta - t);
        t0 = t;
        a = ta;

        /* use secant method to discover a value that unrefracts to ta.
         * max=7 ave=2.4 loops in hundreds of test cases.
         */
        while (1) {
            a += d;
            unrefract (pr, tr, a, &t);
            if (fabs(ta-t) <= MAXRERR)
                break;
            d *= -(ta - t)/(t0 - t);
            t0 = t;
        }

        *aa = a;

#undef  MAXRERR
}

/* For RCS Only -- Do Not Edit */
static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: refract.c,v $ $Date: 2008-03-25 17:45:18 $ $Revision: 1.8 $ $Name: not supported by cvs2svn $"};