source: BAORadio/libindi/v1.0.1/BAOControl/astro.cpp

#include "astro.h"

////////////////////
// Routines Astro //
////////////////////

Astro::Astro()
{
    Longitude=0.0;
    Latitude=0.0;
    Lst=0.0;
    JJ=0.0;
    CorP=0.0;
    CorEP=0.0;
    Annee=0.0;
    Mois=0.0;
    Jour=0.0;
    Heu=0.0;
    Min=0.0;
    Sec=0.0;
    ep=0.0;
    hs=0.0;
    UTCP=0.0;
    tsl=0.0;
}

Astro::~Astro()
{
}


//0<=Angle<=2*PI
double Astro::VerifAngle(double Angle)
{
    Angle=fmod(Angle, Pi2);
    if (Angle<0.0) Angle+=Pi2;

    return(Angle);
}


//////////////////////////////////////////////////////////////////////
//Nutation

void Astro::Nutation()
{
    double T,L,LP,M,MP,O,L2,LP2,O2;

    T=(JJ-2415020.0)/36525.0;

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

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

    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);

    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;
}


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

void Astro::CalculJJ(double Heure)
{

    long y, a, b, c, e, m;

    long year = (long)Annee;
    int month = (int)Mois;
    double day = Jour + 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;

    JJ = b + c + e + day - 32167.5;
}


///////////////////////////////////////////////////////////////////////
//Obliquité

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;
}

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

double Astro::TSL(double JJ, double HeureSiderale, double Longitude)
{
    double rd,T;

    T = (JJ - 2415020.0 ) / 36525.0;
    rd = 0.276919398 + ( 100.0021359 + 0.000001075 * T ) * T;
    rd += HeureSiderale;
    rd *= Pi2;
    rd += CorP * cos(ep);
    rd -= Longitude;

    return VerifAngle(rd);
}


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

    hs-=UTCP/24.0;

    CalculJJ(hs);

    Obliquite(JJ);

    Nutation();

    tsl=TSL(JJ, hs, Longitude);
}


///////////////////////////////////////////////////////////////////////
//Azimut et hauteur à partir de AR et Dec

void Astro::Azimut(double ts, double Latitude, double Ar, double De, double *azi, double *hau)
{
    double ah, ht, a1, az, zc, ac;

    ah=ts-Ar;

    zc=sin(Latitude)*sin(De)+cos(Latitude)*cos(De)*cos(ah);
    ht=atan(zc/sqrt((-zc*zc)+1.0));

    a1=cos(De)*sin(ah)/cos(ht);
    ac=(-cos(Latitude)*sin(De)+sin(Latitude)*cos(De)*cos(ah))/cos(ht);
    az=atan(a1/sqrt((-a1*a1)+1.0));
    if (ac<0.0) az=Pi-az;

    az=VerifAngle(Pi+az);

    *azi=az;
    *hau=ht;
}