source: BAORadio/libindi/libindi/communs/astro.cpp @ 682

////////////////////////////////
// Classe Astro               //
// Franck RICHARD             //
// franckrichard033@gmail.com //
// BAORadio                   //
// juin 2011                  //
////////////////////////////////


#include "astro.h"


/**************************************************************************************
** Constructeur de la classe Astro
**
***************************************************************************************/

Astro::Astro()
{
    //Initialisation des variables

    Longitude = 0.0;                     // Cordonnées du lieu d'observation
    Latitude  = 0.0;

    Annee     = 0.0;                     // Date et heure de l'observation
    Mois      = 0.0;
    Jour      = 0.0;
    Heure     = 0.0;
    Min       = 0.0;
    Sec       = 0.0;
    UTCP      = 0.0;                     // permet de gérer le fuseau horaire
    // mais semble inutile sous Linux

    hs        = 0.0;                     // fraction du jour en cours avec hs=(heure+min/60+secs/3600)/24.
    ep        = 0.0;                     // obliquité de l'écliptique
    tsl       = 0.0;                     // temps sidéral local
    JJ        = 0.0;                     // Jour Julien

    CorP      = 0.0;                     // Nutation en longitude
    CorEP     = 0.0;                     // Nutation en inclinaison

    Pression  = 1013.0;
    Temp      = 0.0;

    LongitudeSoleil = 0.0;
}


/**************************************************************************************
** Constructeur de la classe Astro
**
***************************************************************************************/

Astro::~Astro()
{
}


/**************************************************************************************
** Initialisation des variables globales de la classe Astro
**
***************************************************************************************/

void Astro::DefinirDateHeure(double A, double M, double J, double H, double Mi, double S)
{
    Annee = A;
    Mois  = M;
    Jour  = J;
    Heure = H;
    Min   = Mi;
    Sec   = S;
}

void Astro::DefinirPressionTemp(double P, double T)
{
    Pression = P;
    Temp     = T;
}

void Astro::DefinirLongitudeLatitude(double log, double lat)
{
    Longitude = log;
    Latitude  = lat;
}



/**************************************************************************************
** ParamÚtre : angle en radians
** retourne angle dans un intervalle 0 <= Angle <= 2*PI
***************************************************************************************/

double Astro::VerifAngle(double Angle)
{
    Angle=fmod(Angle, Pi2);
    if (Angle<0.0) Angle+=Pi2;

    return(Angle);
}

/**************************************************************************************
** ParamÚtre : angle en radians
** retourne angle dans un intervalle 0 <= Angle <= PI
***************************************************************************************/

double Astro::VerifDistance(double Angle)
{
    Angle=VerifAngle(Angle);
    if (Angle>Pi) Angle=Pi2-Angle;

    return(Angle);
}


/**************************************************************************************
** Distance entre deux points situés sur une sphÚre
**
***************************************************************************************/

double Astro::DistanceAngulaireEntre2Points(double az1, double ha1, double az2, double ha2)
{
    return acos(sin(ha1) * sin(ha2) + cos(ha1) * cos(ha2) * cos(az2 - az1));
}


/**************************************************************************************
** retourne l'arrondi du nombre donné en paramÚtre
** Exemples: Arrondi(5.7) =  6
**           Arrondi(-5.7)= -6
***************************************************************************************/

double Astro::Arrondi(double a)
{
    double b;

    if (a>=0.0)
    {
        b=a-floor(a);
        if (b>=0.5) return ceil(a);
        else return floor(a);
    }
    else
    {
        b=ceil(a)-a;
        if (b>=0.5) return floor(a);
        else return ceil(a);
    }
}



/**************************************************************************************
** converti un angle (exprimé en degrés) dans le format hh:mm:ss ou deg:mm:ss
** selon la valeur de HMS
***************************************************************************************/

string Astro::DHMS(double mema, bool HMS)
{
    double d, m, s, a;

    a=mema;

    if (HMS) a = a / 15.0;
    a=fabs(a * 3600.0);
    s=Arrondi(fmod(a, 60.0) * 10.0) / 10.0;
    a=floor(a / 60.0);
    m=fmod(a, 60.0);
    d=floor(a / 60.0);

    if (s == 60.0)
    {
        s=0.0;
        m+=1.0;
        if (m==60.0)
        {
            m=0.0;
            d+=1.0;

            if (HMS && d>23.0) d=0.0;
        }
    }

    stringstream strsTemp;

    if (mema < 0) strsTemp << '-';

    strsTemp << setfill('0') << setw(2) << d << ":" << setfill('0') << setw(2) << m << ":" << setfill('0') << setw(2) << s;

    return (strsTemp.str());
}

/**************************************************************************************
** Décomposition et vérification des dates, heures, AD et déclinaisons
**
** type = 0 -> chaîne contient une déclinaison ou une latitude
** type = 1 -> "        "      une longitude
** type = 2 -> "        "      une AD ou une heure
** type = 3 -> "        "      une date
**
** Exemple de Décomposition("15:23:12", 2, a, b, c) retourne a=15, b=23, c=12
** si la chaîne a un format incorrect, la fct retourne false
**
***************************************************************************************/

bool Astro::Decomposition(string chaine, char type, float *a1, float *a2, float *a3)
{
    string car, s;

    float a, b, c;

    // pour les heures et les coordonnées, on attend ":" comme caractÚre séparateur, sinon
    // on attend d'avoir des caractÚres "/" entre chaque éléments pour une date

    (type == 3) ? car="/" : car=":";

    // Y a-t-il 2 caractÚres ':' ou '/' dans la chaîne ?
    // C'est indispensable dans tous les cas

    int test = 0;

    for (size_t i=0; i<chaine.length(); i++) if (chaine[i] == car[0]) test++;

    if (test<2) return false;

    // Extraction des trois nombres

    s = chaine.substr(0, chaine.find(car));

    a = atoi(s.c_str());

    s = chaine.substr(chaine.find(car) + 1, chaine.rfind(car) - chaine.find(car) - 1);

    b = atoi(s.c_str());

    s = chaine.substr(chaine.rfind(car)+1);

    // 12.125454414 sec  -> 12.125 sec

    c = Arrondi(atof(s.c_str()) * 100.0) / 100.0;

    //vérification de la cohérence des infos contenues dans la chaîne

    if (type < 3 )
    {
        // pour une déclinaison
        if (( type == 0 ) && ( a>90.0  || a<-90.0 )) return false;
        // pour une AD
        if (( type == 1 ) && ( a>360.0 || a<0.0   )) return false;
        // pour une heure
        if (( type == 2 ) && ( a>23.0  || a<0.0   )) return false;
        // pour les minutes
        if ( b < 0.0 || b> 59.0 ) return false;
        //pour les secondes
        if ( c < 0.0 || c>=60.0 ) return false;
    }
    else
    {
        //pour les jours
        if ( a < 0.0 || a > 31.0 ) return false;
        //pour les mois
        if ( b < 0.0 || b > 12.0 ) return false;
    }

    if ( a1 != NULL ) *a1 = a;
    if ( a2 != NULL ) *a2 = b;
    if ( a3 != NULL ) *a3 = c;

    return true;
}

/**************************************************************************************
** Calcul de la longitude du Soleil
** La quantité PS n'est pas indispensable. Elle permet d'améliorer la précision
***************************************************************************************/

double Astro::CalculLongitudeSoleil()
{
    double T = (JJ - 2415020.0) / 36525.0;

    double A = 153.23 + 22518.7541 * T;
    double B = 216.57 + 44037.5082 * T;
    double C = 312.69 + 32964.3577 * T;
    double D = 350.74 + ( 445267.1142 - 0.00144 * T ) * T;
    double E = 231.19 + 20.2 * T;


    double L = 279.69668 + ( 36000.76892 + 0.0003025 * T ) * T;
    double M = 358.47583 + (((( -0.0000033 * T) -0.00015 )  * T ) + 35999.04975 ) * T;

    double CC =   (1.919460 + ((-0.000014 * T) - 0.004789 ) * T ) * sin(M * Pidiv180)
                  + (0.020094 - 0.0001 * T) * sin(2.0 * M * Pidiv180)
                  + 0.000293 * sin(3.0 * M * Pidiv180);

    double PS =  0.00134 * cos(A * Pidiv180)
                 +0.00154 * cos(B * Pidiv180)
                 +0.00200 * cos(C * Pidiv180)
                 +0.00179 * sin(D * Pidiv180)
                 +0.00178 * sin(E * Pidiv180);


    double LS = VerifAngle( (L + CC + PS) * Pidiv180);

    return LS;
}



/**************************************************************************************
** Calcul de l'AD et la déclinaison du Soleil
** indispensable pour exécuter la commande "goto sun"
***************************************************************************************/

void Astro::CalculARDecSoleil(CoordonneesHorairesDouble *Soleil)
{
    double ar;
    double de;

    ar = atan( cos(ep) * sin(LongitudeSoleil) / cos(LongitudeSoleil));

    if ( cos(LongitudeSoleil) < 0.0) ar += Pi;

    de = asin( sin(ep) * sin(LongitudeSoleil));

    Soleil->ar = VerifAngle(ar);

    Soleil->dec = de;
}




/**************************************************************************************
** Calcul de la nutation en longitude et en inclinaison
** Nécessaire pour calculer la position apparente d'un objet (étoile, galaxie etc...)
***************************************************************************************/

void Astro::Nutation()
{
    double T = (JJ-2415020.0) / 36525.0;

    double L  = (279.6967 + (36000.7689 + 0.000303 * T ) * T ) * Pidiv180;
    double LP = (270.4342 + (481267.8831 - 0.001133 * T ) * T ) * Pidiv180;
    double M  = (358.4758 + (35999.0498 - 0.000150 * T ) * T ) * Pidiv180;
    double MP = (296.1046 + (477198.8491 + 0.009192 * T ) * T ) * Pidiv180;
    double O  = (259.1833 + (-1934.1420 + 0.002078 * T ) * T ) * Pidiv180;

    double L2  = 2.0 * L;
    double LP2 = 2.0 * LP;
    double O2  = 2.0 * O;

    // Nutation en longitude

    CorP = -(17.2327 + 0.01737 * T) * sin(O)
           -(1.2729+0.00013 * T) * sin(L2)
           +0.2088 * sin(O2)
           -0.2037 * sin(LP2)
           +(0.1261 - 0.00031 * T) * sin(M)
           +0.0675 * sin(MP)
           -(0.0497 - 0.00012 * T) * sin(L2 + M)
           -0.0342 * sin(LP2 - O)
           -0.0261 * sin(LP2 + MP)
           +0.0214 * sin(L2 - M)
           -0.0149 * sin(L2 - LP2 + MP)
           +0.0124 * sin(L2 - O)
           +0.0114 * sin(LP2 - MP);

    //Nutation en inclinaison

    CorEP = (9.2100 + 0.00091 * T) * cos(O)
            +(0.5522 - 0.00029 * T) * cos(L2)
            -0.0904 * cos(O2)
            +0.0884 * cos(LP2)
            +0.0216 * cos(L2+M)
            +0.0183 * cos(LP2-O)
            +0.0113 * cos(LP2+MP)
            -0.0093 * cos(L2-M)
            -0.0066 * cos(L2-O);

    CorP  = CorP / 3600.0 * Pidiv180;
    CorEP = CorEP / 3600.0 * Pidiv180;
}



/**************************************************************************************
** Nutation des étoiles - on utilise les quantités calculées précédemment
** ar et de sont exprimés en radians
***************************************************************************************/

void Astro::NutationEtoile(double *ar, double *de)
{
    double a = (cos(ep) + sin(ep) * sin(*ar) * tan(*de)) * CorP - (cos(*ar) * tan(*de)) * CorEP;
    double b = (sin(ep) * cos(*ar)) * CorP + sin(*ar) * CorEP;

    *ar = VerifAngle(*ar+a);
    *de += b;
}


/**************************************************************************************
** Précession des équinoxes
** ar et de sont exprimés en radians
***************************************************************************************/

void Astro::Precession(double *ar, double *de)
{
    double t = (JJ - 2451544.5) / 36524.2199;

    double ze = ((((0.018 * t) + 0.302) * t + 2304.948) * t) / 3600.0 * Pidiv180;
    double z = ((((0.019 * t) + 1.093) * t + 2304.948) * t) / 3600.0 * Pidiv180;
    double delta = ((((-0.042 * t) - 0.426) * t + 2004.255) * t) / 3600.0 * Pidiv180;


    double A = cos(*de) * sin(*ar + ze);
    double B = cos(delta) * cos(*de) * cos (*ar + ze) - sin(delta) * sin(*de);
    double C = sin(delta) * cos(*de) * cos (*ar + ze) + cos(delta) * sin(*de);

    double AMZ = atan(A / B);

    if (B<0.0) AMZ += Pi;

    *ar=VerifAngle(AMZ + z);
    *de=asin(C);
}



/**************************************************************************************
** Obliquité
** calcul de l'inclinaison de l'axe terrestre par rapport à au plan de l'écliptique
** Quantité indispensable pour calculer le temps sidéral local et les coordonnées
** horaires du soleil
***************************************************************************************/

void Astro::Obliquite(double JJ)
{
    double T;

    T = (JJ - 2415020.0) / 36525.0;

    ep = (23.452294+ (((0.000000503 * T ) - 0.00000164 ) * T - 0.0130125 ) * T ) * Pidiv180;
}




/**************************************************************************************
** Aberration annuelle des étoiles
** ar et de sont exprimés en radians
***************************************************************************************/

void Astro::AberrationAnnuelle(double *ar, double *de)
{
    double c = -20.49 / 3600.0 * Pidiv180;
    double a = c * ( cos(*ar) * cos(LongitudeSoleil) * cos(ep)
                     + sin(*ar) * sin(LongitudeSoleil)
                   ) / cos(*de);

    double b = c * ( cos(LongitudeSoleil) * cos(ep) * (tan(ep) * cos(*de) - sin(*ar) * sin(*de))
                     + cos(*ar) * sin(*de) * sin(LongitudeSoleil));

    *ar = VerifAngle(*ar + a);
    *de += b;
}


/**************************************************************************************
** calcul Jour julien
***************************************************************************************/

void Astro::CalculJJ(double Heure)
{
    JJ = CalculJJ(Annee, Mois, Jour, Heure);
}

double Astro::CalculJJ(double A, double M, double J, double Heure)
{
    long y, a, b, c, e, m;

    long year = (long)A;
    int month = (int)M;
    double day = J + Heure;

    y = year + 4800;

    m = month;
    if ( m <= 2 )
    {
        m += 12;
        y -= 1;
    }
    e = (306 * (m+1))/10;

    a = y/100;
    if ( year <= 1582L )
    {
        if ( year == 1582L )
        {
            if ( month < 10 )
                goto julius;
            if ( month > 10)
                goto gregor;
            if ( day >= 15 )
                goto gregor;
        }
julius:

        b = -38;
    }
    else
    {

gregor:
        b = (a / 4) - a;
    }

    c = (36525L * y) / 100;

    return b + c + e + day - 32167.5;
}



/**************************************************************************************
** Temps sidéral local
***************************************************************************************/

double Astro::TSL(double JJ, double HeureSiderale, double Longitude)
{
    double T = (JJ /*- 66.185/3600.0/24.0*/ - 2415020.0 ) / 36525.0;
    double rd = 0.276919398 + ( 100.0021359 + 0.000001075 * T ) * T;
    rd += HeureSiderale;
    rd *= Pi2;
    // temps sidéral apparent
    rd += CorP * cos(ep);
    rd -= Longitude;

    return VerifAngle(rd);
}


/**************************************************************************************
** routine principale des calculs
***************************************************************************************/

void Astro::CalculTSL()
{
    hs = (Heure + Min / 60.0 + Sec / 3600.0) / 24.0;

    hs -= UTCP / 24.0;

    CalculJJ(hs);

    Obliquite(JJ);

    Nutation();

    tsl = TSL(JJ, hs, Longitude);

    LongitudeSoleil = CalculLongitudeSoleil();
}



/**************************************************************************************
** Calcule la hauteur et l'azimut des étoiles en fct du lieu d'observation
***************************************************************************************/

 void Astro::Azimut(double Ar, double De, double *azi, double *hau)
{
    double ah = tsl - Ar;
    
    double cosLatitude = cos(Latitude);
    double sinLatitude = sin(Latitude);
    double cosDe_cosAh = cos(De) * cos(ah);
    double sinDe       = sin(De);
   
    
    double zc = sinLatitude * sinDe + cosLatitude * cosDe_cosAh;
    double ht = atan(zc / sqrt((-zc * zc) + 1.0));
    
    double cosHt       = cos(ht);

    double a1 = cos(De) * sin(ah) / cosHt;
    double ac = (-cosLatitude * sinDe + sinLatitude * cosDe_cosAh) / cosHt;
    double az = atan(a1 / sqrt((-a1*a1)+1.0));

    if (ac<0.0) az = Pi - az;

    *azi = VerifAngle( Pi + az );
    *hau = ht;
}


/**************************************************************************************
** Calcule l'ascension droite et la déclinaison d'un objet situé à l'azimut az
** et à la déclinaison de
***************************************************************************************/

void Astro::AzHa2ADDe(double az, double ha, double *AD, double *Dec)
{
    double sindel, del, habar, cosha, ha10;

    sindel = sin(ha) * sin(Latitude) + cos(ha) * cos(Latitude) * cos(az);

    del = asin(sindel);

    if ( (sindel == 1.0) || ( sindel == -1.0) )
    {
        habar = 0.0;
    }
    else
    {
        if ( (sin(Latitude) == 1.0) || (sin(Latitude) == -1.0))
        {
            habar = az;
        }
        else
        {
            cosha = ( sin(ha) - sin(Latitude) * sindel ) / ( cos(Latitude) * cos(del) );
            if (cosha > 1.0  )  cosha =  1.0;
            if (cosha < -1.0 )  cosha = -1.0;

            habar = acos(cosha);
        }
    }

    if (sin(az) < 0.0)
    {
        ha10 = habar;
    }
    else
    {
        ha10 = Pi2 - habar ;
    }

    *AD = VerifAngle( tsl - ha10 );

    // if (*AD < 0.0) *AD += Pi2;

    *Dec = del;
}


/**************************************************************************************
** Recherche l'azimut de l'objet suivi au moment où celui-ci passera sous les 30°
** et deviendra donc hors d'atteinte de l'instrument.
***************************************************************************************/

void Astro::RechercheAzimutFinSuivi(double AD, double Dec, long int *azmin, long int *azmax)
{
    double az, ha;
    double min =  10000.0;
    double max = -10000.0;

    CalculTSL();

    for (int i=0; i < 120; i++)
    {
        tsl += i * 5.0 / 60.0 / 24.0 ;

        Azimut( AD, Dec, &az, &ha);

        if ( az > max ) max = az;

        if ( az < min ) min = az;

        if ( ha < HAUTMIN * Pidiv180) break;
    }

    *azmin = (long int) Arrondi( min * (double)NBREPASCODEURSAZ / 360.0) - NBREPASCODEURSAZ;

    *azmax = (long int) Arrondi( max * (double)NBREPASCODEURSAZ / 360.0) + NBREPASCODEURSAZ;

    CalculTSL();
}



/**************************************************************************************
** réfraction atmosphérique : formule simple de Jean Meeus
** la hauteur ht est exprimée en radians
***************************************************************************************/

double Astro::RefractionAtmospherique(double ht)
{
    double gamma, R0, a, b, PressionTempCalc, tanHT;

    if (Pression != 0.0)
    {
        PressionTempCalc = (Pression / 1013.25) * 273.0 / (273.0 + Temp);

        if (ht <= 15.0 * Pidiv180)
        {
            gamma = 2.6;
            a = 7.5262;
            b = -2.2204;

            if (ht >= 4.0 * Pidiv180)
            {
                gamma=-1.1;
                a=4.4010;
                b=-0.9603;
            }

            R0 = (pow((a+b*log(ht*N180divPi+gamma)), 2.0)) / 60.0 * Pidiv180;
            ht += R0 * PressionTempCalc;
        }
        else
        {
            tanHT = fabs(tan(Pidiv2 - ht));
            R0 = (0.0161877777777777777777-0.000022888888888888888 * tanHT * tanHT)
                 * tanHT * Pidiv180;
            ht += R0 * PressionTempCalc;
        }
    }

    return ht;
}




double Astro::slaDrange(double angle )
{
    return fmod(angle, Pi);
}

float  Astro::slaRange(float angle )
{
    return fmodf(angle, Pi);
}

void Astro::slaRefro ( double zobs, double hm, double tdk, double pmb,
                       double rh, double wl, double phi, double tlr,
                       double eps, double *ref )
/*
**  - - - - - - - - -
**   s l a R e f r o
**  - - - - - - - - -
**
**  Atmospheric refraction for radio and optical wavelengths.
**
**  Given:
**    zobs    double  observed zenith distance of the source (radian)
**    hm      double  height of the observer above sea level (metre)
**    tdk     double  ambient temperature at the observer (deg K)
**    pmb     double  pressure at the observer (millibar)
**    rh      double  relative humidity at the observer (range 0-1)
**    wl      double  effective wavelength of the source (micrometre)
**    phi     double  latitude of the observer (radian, astronomical)
**    tlr     double  temperature lapse rate in the troposphere (degK/met
**    eps     double  precision required to terminate iteration (radian)
**
**  Returned:
**    ref     double  refraction: in vacuo ZD minus observed ZD (radian)
**
**  Notes:
**
**  1  A suggested value for the tlr argument is 0.0065D0.  The
**     refraction is significantly affected by tlr, and if studies
**     of the local atmosphere have been carried out a better tlr
**     value may be available.
**
**  2  A suggested value for the eps argument is 1e-8.  The result is
**     usually at least two orders of magnitude more computationally
**     precise than the supplied eps value.
**
**  3  The routine computes the refraction for zenith distances up
**     to and a little beyond 90 deg using the method of Hohenkerk
**     and Sinclair (NAO Technical Notes 59 and 63, subsequently adopted
**     in the Explanatory Supplement, 1992 edition - see section 3.281).
**
**  4  The C code is an adaptation of the Fortran optical refraction
**     subroutine AREF of C.Hohenkerk (HMNAO, September 1984), with
**     extensions to support the radio case.  The following modifications
**     to the original HMNAO optical refraction algorithm have been made:
**
**     .  The angle arguments have been changed to radians.
**
**     .  Any value of zobs is allowed (see note 6, below).
**
**     .  Other argument values have been limited to safe values.
**
**     .  Murray's values for the gas constants have been used
**        (Vectorial Astrometry, Adam Hilger, 1983).
**
**     .  The numerical integration phase has been rearranged for
**        extra clarity.
**
**     .  A better model for Ps(T) has been adopted (taken from
**        Gill, Atmosphere-Ocean Dynamics, Academic Press, 1982).
**
**     .  More accurate expressions for Pwo have been adopted
**        (again from Gill 1982).
**
**     .  Provision for radio wavelengths has been added using
**        expressions devised by A.T.Sinclair, RGO (private
**        communication 1989), based on the Essen & Froome
**        refractivity formula adopted in Resolution 1 of the
**        13th International Geodesy Association General Assembly
**        (Bulletin Geodesique 1963 p390).
**
**     .  Various small changes have been made to gain speed.
**
**     None of the changes significantly affects the optical results
**     with respect to the algorithm given in the 1992 Explanatory
**     Supplement.  For example, at 70 deg zenith distance the present
**     routine agrees with the ES algorithm to better than 0.05 arcsec
**     for any reasonable combination of parameters.  However, the
**     improved water-vapour expressions do make a significant difference
**     in the radio band, at 70 deg zenith distance reaching almost
**     4 arcsec for a hot, humid, low-altitude site during a period of
**     low pressure.
**
**  5  The radio refraction is chosen by specifying wl > 100 micrometres.
**     Because the algorithm takes no account of the ionosphere, the
**     accuracy deteriorates at low frequencies, below about 30 MHz.
**
**  6  Before use, the value of zobs is expressed in the range +/- pi.
**     If this ranged zobs is -ve, the result ref is computed from its
**     absolute value before being made -ve to match.  In addition, if
**     it has an absolute value greater than 93 deg, a fixed ref value
**     equal to the result for zobs = 93 deg is returned, appropriately
**     signed.
**
**  7  As in the original Hohenkerk and Sinclair algorithm, fixed values
**     of the water vapour polytrope exponent, the height of the
**     tropopause, and the height at which refraction is negligible are
**     used.
**
**  8  The radio refraction has been tested against work done by
**     Iain Coulson, JACH, (private communication 1995) for the
**     James Clerk Maxwell Telescope, Mauna Kea.  For typical conditions,
**     agreement at the 0.1 arcsec level is achieved for moderate ZD,
**     worsening to perhaps 0.5-1.0 arcsec at ZD 80 deg.  At hot and
**     humid sea-level sites the accuracy will not be as good.
**
**  9  It should be noted that the relative humidity rh is formally
**     defined in terms of "mixing ratio" rather than pressures or
**     densities as is often stated.  It is the mass of water per unit
**     mass of dry air divided by that for saturated air at the same
**     temperature and pressure (see Gill 1982).
**
**  Called:  slaDrange, atmt, atms
**
**  Defined in slamac.h:  TRUE, FALSE
**
**  Last revision:   30 January 1997
**
**  Copyright P.T.Wallace.  All rights reserved.
*/
{
    /* Fixed parameters */

    static double d93 = 1.623156204; /* 93 degrees in radians        */
    static double gcr = 8314.32;     /* Universal gas constant       */
    static double dmd = 28.9644;     /* Molecular weight of dry air  */
    static double dmw = 18.0152;     /* Molecular weight of water
                                                             vapour */
    static double s = 6378120.0;     /* Mean Earth radius (metre)    */
    static double delta = 18.36;     /* Exponent of temperature
                                         dependence of water vapour
                                                           pressure */
    static double ht = 11000.0;      /* Height of tropopause (metre) */
    static double hs = 80000.0;      /* Upper limit for refractive
                                                    effects (metre) */

    /* Variables used when calling the internal routine atmt */
    double robs;   /* height of observer from centre of Earth (metre) */
    double tdkok;  /* temperature at the observer (deg K) */
    double alpha;  /* alpha          |        */
    double gamm2;  /* gamma minus 2  | see ES */
    double delm2;  /* delta minus 2  |        */
    double c1,c2,c3,c4,c5,c6;  /* various */

    /* Variables used when calling the internal routine atms */
    double rt;     /* height of tropopause from centre of Earth (metre) */
    double tt;     /* temperature at the tropopause (deg k) */
    double dnt;    /* refractive index at the tropopause */
    double gamal;  /* constant of the atmospheric model = g*md/r */

    int is, k, n, i, j, optic;
    double zobs1, zobs2, hmok, pmbok, rhok, wlok, tol, wlsq, gb,
           a, gamma, tdc, psat, pwo, w, tempo, dn0, rdndr0, sk0,
           f0, rdndrt, zt, ft, dnts, rdndrp, zts, fts, rs,
           dns, rdndrs, zs, fs, refold, z0, zrange, fb, ff, fo,
           fe, h, r, sz, rg, dr, tg, dn, rdndr, t, f, refp, reft;

    /* The refraction integrand */
#define refi(R,DN,RDNDR) ((RDNDR)/(DN+RDNDR));



    /* Transform zobs into the normal range */
    zobs1 = slaDrange ( zobs );
    zobs2 = fabs ( zobs1 );
    zobs2 = gmin ( zobs2, d93 );

    /* Keep other arguments within safe bounds */
    hmok = gmax ( hm, -1000.0 );
    hmok = gmin ( hmok, 10000.0 );
    tdkok = gmax ( tdk, 100.0 );
    tdkok = gmin ( tdkok, 500.0 );
    pmbok = gmax ( pmb, 0.0 );
    pmbok = gmin ( pmbok, 10000.0 );
    rhok  = gmax ( rh, 0.0 );
    rhok  = gmin ( rhok, 1.0 );
    wlok  = gmax ( wl, 0.1 );
    alpha = fabs ( tlr );
    alpha = gmax ( alpha, 0.001 );
    alpha = gmin ( alpha, 0.01 );

    /* Tolerance for iteration */
    w = fabs ( eps );
    tol = gmin ( w, 0.1 ) / 2.0;

    /* Decide whether optical or radio case - switch at 100 micron */
    optic = ( wlok <= 100.0 );

    /* Set up model atmosphere parameters defined at the observer */
    wlsq = wlok * wlok;
    gb = 9.784 * ( 1.0 - 0.0026 * cos ( 2.0 * phi ) - 2.8e-7 * hmok );
    a = ( optic ) ?
        ( ( 287.604 + 1.6288 / wlsq + 0.0136 / ( wlsq * wlsq ) )
          * 273.15 / 1013.25 ) * 1e-6
        :
        77.624e-6;
    gamal = gb * dmd / gcr;
    gamma = gamal / alpha;
    gamm2 = gamma - 2.0;
    delm2 = delta - 2.0;
    tdc = tdkok - 273.15;
    psat = pow ( 10.0, ( 0.7859 + 0.03477 * tdc ) /
                 ( 1.0 + 0.00412 * tdc ) ) *
           ( 1.0 + pmbok * ( 4.5e-6 + 6e-10 * tdc * tdc ) );
    pwo = ( pmbok > 0.0 ) ?
          rhok * psat / ( 1.0 - ( 1.0 - rhok ) * psat / pmbok ) :
          0.0;
    w = pwo * ( 1.0 - dmw / dmd ) * gamma / ( delta - gamma );
    c1 = a * ( pmbok + w ) / tdkok;
    c2 = ( a * w + ( optic ? 11.2684e-6 : 12.92e-6 ) * pwo ) / tdkok;
    c3 = ( gamma - 1.0 ) * alpha * c1 / tdkok;
    c4 = ( delta - 1.0 ) * alpha * c2 / tdkok;
    c5 = optic ? 0.0 : 371897e-6 * pwo / tdkok;
    c6 = c5 * delm2 * alpha / ( tdkok * tdkok );

    /* Conditions at the observer */
    robs = s + hmok;
    atmt ( robs, tdkok, alpha, gamm2, delm2, c1, c2, c3, c4, c5, c6, robs,
           &tempo, &dn0, &rdndr0 );
    sk0 = dn0 * robs * sin ( zobs2 );
    f0 = refi ( robs, dn0, rdndr0 );

    /* Conditions at the tropopause in the troposphere */
    rt = s + ht;
    atmt ( robs, tdkok, alpha, gamm2, delm2, c1, c2, c3, c4, c5, c6, rt,
           &tt, &dnt, &rdndrt );
    zt = asin ( sk0 / ( rt * dnt ) );
    ft = refi ( rt, dnt, rdndrt );

    /* Conditions at the tropopause in the stratosphere */
    atms ( rt, tt, dnt, gamal, rt, &dnts, &rdndrp );
    zts = asin ( sk0 / ( rt * dnts ) );
    fts = refi ( rt, dnts, rdndrp );

    /* At the stratosphere limit */
    rs = s + hs;
    atms ( rt, tt, dnt, gamal, rs, &dns, &rdndrs );
    zs = asin ( sk0 / ( rs * dns ) );
    fs = refi ( rs, dns, rdndrs );

    /*
    ** Integrate the refraction integral in two parts;  first in the
    ** troposphere (k=1), then in the stratosphere (k=2).
    */

    /* Initialize previous refraction to ensure at least two iterations */
    refold = 1e6;

    /*
    ** Start off with 8 strips for the troposphere integration, and then
    ** use the final troposphere value for the stratosphere integration,
    ** which tends to need more strips.
    */
    is = 8;

    /* Troposphere then stratosphere */
    for ( k = 1; k <= 2; k++ ) {

        /* Start z, z range, and start and end values */
        if ( k == 1 ) {
            z0 = zobs2;
            zrange = zt - z0;
            fb = f0;
            ff = ft;
        } else {
            z0 = zts;
            zrange = zs - z0;
            fb = fts;
            ff = fs;
        }

        /* Sums of odd and even values */
        fo = 0.0;
        fe = 0.0;

        /* First time through loop we have to do every point */
        n = 1;

        /* Start of iteration loop (terminates at specified precision) */
        for ( ; ; ) {

            /* Strip width */
            h = zrange / (double) is;

            /* Initialize distance from Earth centre for quadrature pass */
            r = ( k == 1 ) ? robs : rt;

            /* One pass (no need to compute evens after first time) */
            for ( i = 1; i < is; i += n ) {

                /* Sine of observed zenith distance */
                sz = sin ( z0 + h * (double) i );

                /* Find r (to nearest metre, maximum four iterations) */
                if ( sz > 1e-20 ) {
                    w = sk0 / sz;
                    rg = r;
                    j = 0;
                    do {
                        if ( k == 1 ) {
                            atmt ( robs, tdkok, alpha, gamm2, delm2,
                                   c1, c2, c3, c4, c5, c6, rg,
                                   &tg, &dn, &rdndr );
                        } else {
                            atms ( rt, tt, dnt, gamal, rg, &dn, &rdndr );
                        }
                        dr = ( rg * dn - w ) / ( dn + rdndr );
                        rg -= dr;
                    } while ( fabs ( dr ) > 1.0 && j++ <= 4 );
                    r = rg;
                }

                /* Find refractive index and integrand at r */
                if ( k == 1 ) {
                    atmt ( robs, tdkok, alpha, gamm2, delm2,
                           c1, c2, c3, c4, c5, c6, r,
                           &t, &dn, &rdndr );
                } else {
                    atms ( rt, tt, dnt, gamal, r, &dn, &rdndr );
                }
                f = refi ( r, dn, rdndr );

                /* Accumulate odd and (first time only) even values */
                if ( n == 1 && i % 2 == 0 ) {
                    fe += f;
                } else {
                    fo += f;
                }
            }

            /* Evaluate the integrand using Simpson's Rule */
            refp = h * ( fb + 4.0 * fo + 2.0 * fe + ff ) / 3.0;

            /* Has the required precision been reached? */
            if ( fabs ( refp - refold ) > tol ) {

                /* No: prepare for next iteration */
                refold = refp;   /* Save current value for convergence test */
                is += is;        /* Double the number of strips */
                fe += fo;        /* Sum of all = sum of evens next time */
                fo = 0.0;        /* Reset odds accumulator */
                n = 2;           /* Skip even values next time */

            } else {

                /* Yes: save troposphere component and terminate loop */
                if ( k == 1 ) reft = refp;
                break;
            }
        }
    }

    /* Result */
    *ref = reft + refp;
    if ( zobs1 < 0.0 ) *ref = - ( *ref );

    *ref+=Pidiv2-zobs1;
}

/*--------------------------------------------------------------------------*/

void Astro::atmt ( double robs, double tdkok, double alpha, double gamm2,
                   double delm2, double c1, double c2, double c3,
                   double c4, double c5, double c6, double r,
                   double *t, double *dn, double *rdndr )
/*
**  - - - - -
**   a t m t
**  - - - - -
**
**  Internal routine used by slaRefro:
**
**    refractive index and derivative with respect to height for the
**    troposphere.
**
**  Given:
**    robs    double   height of observer from centre of the Earth (metre)
**    tdkok   double   temperature at the observer (deg K)
**    alpha   double   alpha          )
**    gamm2   double   gamma minus 2  ) see ES
**    delm2   double   delta minus 2  )
**    c1      double   useful term  )
**    c2      double   useful term  )
**    c3      double   useful term  ) see source of
**    c4      double   useful term  ) slaRefro main routine
**    c5      double   useful term  )
**    c6      double   useful term  )
**    r       double   current distance from the centre of the Earth (metre)
**
**  Returned:
**    *t      double   temperature at r (deg K)
**    *dn     double   refractive index at r
**    *rdndr  double   r * rate the refractive index is changing at r
**
**  This routine is derived from the ATMOSTRO routine (C.Hohenkerk,
**  HMNAO), with enhancements specified by A.T.Sinclair (RGO) to
**  handle the radio case.
**
**  Note that in the optical case c5 and c6 are zero.
*/
{
    double w, tt0, tt0gm2, tt0dm2;

    w = tdkok - alpha * ( r - robs );
    w = gmin ( w, 320.0 );
    w = gmax ( w, 200.0 );
    tt0 = w / tdkok;
    tt0gm2 = pow ( tt0, gamm2 );
    tt0dm2 = pow ( tt0, delm2 );
    *t = w;
    *dn = 1.0 + ( c1 * tt0gm2 - ( c2 - c5 / w ) * tt0dm2 ) * tt0;
    *rdndr = r * ( - c3 * tt0gm2 + ( c4 - c6 / tt0 ) * tt0dm2 );
}

/*--------------------------------------------------------------------------*/

void Astro::atms ( double rt, double tt, double dnt, double gamal, double r,
                   double *dn, double *rdndr )
/*
**  - - - - -
**   a t m s
**  - - - - -
**
**  Internal routine used by slaRefro:
**
**   refractive index and derivative with respect to height for the
**   stratosphere.
**
**  Given:
**    rt      double   height of tropopause from centre of the Earth (metre)
**    tt      double   temperature at the tropopause (deg k)
**    dnt     double   refractive index at the tropopause
**    gamal   double   constant of the atmospheric model = g*md/r
**    r       double   current distance from the centre of the Earth (metre)
**
**  Returned:
**    *dn     double   refractive index at r
**    *rdndr  double   r * rate the refractive index is changing at r
**
**  This routine is derived from the ATMOSSTR routine (C.Hohenkerk, HMNAO).
*/
{
    double b, w;

    b = gamal / tt;
    w = ( dnt - 1.0 ) * exp ( - b * ( r - rt ) );
    *dn = 1.0 + w;
    *rdndr = - r * b * w;
}