Algorithm series (18) using astronomical method to calculate Shiber (lower)

"Take the Last Piece"

The coordinate value is the theoretical value, or the geometrical position of the celestial body, but the FK5 system is a visual system, that is to say, the human eye observation effect (optical position), which needs to be further modified according to the Earth's physical environment, atmosphere environment and other information. So that it is consistent with the observations of human beings observing stars on Earth.

First, a nutation correction is required. Nutation refers to a slight wobble caused by the uneven distribution of material on earth and the perturbation of the Moon and other planets as the Earth rotates along its axis of rotation around the ecliptic's extremely slow rotation. British astronomer James Bradley (1693-1762) found the first nutation, nutation can be decomposed along the ecliptic into horizontal components and vertical components, the level of the ecliptic is recorded as Δψ, called Huang Movement, which affects the longitude of all celestial bodies on the celestial sphere. The vertical component of the ecliptic, known as δε, is called the angle nutation, which affects the obliquity of the intersection. At present, the nutation theory based on the compilation of Astronomical Almanac is founded by Woolard in 1953, and it is based on the rigid Earth model. In 1977, an International Astronomical Union panel of Experts recommended the introduction of a non-rigid Earth model-Molodenski II model to replace the rigid earth model calculation nutation, which was formally adopted by the 17th session of the International Astronomical Union in 1979 and decided to be formally implemented in 1984.

The Earth nutation is mainly caused by the lunar movement, also has a certain periodicity, can be described as a number of periodic items and the period of the main term is 6798.4 days (18.6 years), but the other items are some short period (less than 10 days). The calculation method used in this paper is taken from the IAU1980 nutation theory of the International Astronomical Union, and the periodic term factor data is derived from the table 21-a of the 21st chapter of the astronomical algorithm, which ignores the periodic term of the coefficient less than 0.0003 in the IAU1980 nutation theory, so there are only 63 items. Each cycle term includes the calculation of the sine coefficient (phase Huang) of the Motion (Δψ), calculates the (Δε) cosine coefficient (phase Wei Xiang coefficient) of the intersection nutation and the linear combination coefficients of the 5 basic angular distances (m, M, D, F, Ω) for calculating the angle of convergence. The formula for calculating 5 basic angular distances is:

Angle (angular distance between the sun and the Earth):
D = 297.85036 + 455267.111480 * T-0.0019142 * T2 + t3/189474 (3.10-type)
Sun (Earth) flat near point angle:
M = 357.52772 + 35999.050340 * T-0.0001603 * t2-t3/300000 (3.11-type)
Close point angle of lunar flat
M ' = 134.96298 + 477198.867398 * T + 0.0086972 * T2 + t3/56250 (3.12-type)

Lunar Latitude Parameters:
F = 93.27191 + 483202.017538 * T-0.0036825 * T2 + t3/327270 (3.13-type)
The ecliptic and the lunar orbit rise to the intersection of the Yellow meridian:
ω= 125.04452-1934.136261 * T + 0.0020708 * T2 + t3/450000 (3.14-type)

The above all kinds of T is the number of Confucianism century, calculated the 5 basic angular distance units are degrees, in the calculation of sine or cosine to be converted to radians units. The calculation of the Huang dynamic process for each cycle term is done by first combining the 3.10-3.14 formula with the corresponding 5 basic angular distance coefficients, and calculating the spoke angle. Taking the seventh item in the Nutation Cycle factor table used in this paper as an example, the corresponding coefficients of the 5 basic angular distances are 1, 0, 2, 2 and 2, and the value of the angular θ is: -2d + M + 2F + 2Ω. The value of the period item can be calculated after the spoke angle is calculated:

S = (s1+ S2 * T) * sin (theta) (3.15-type)

Still, for example, in the seventh item, the value of S is ( -517 + 1.2 * T) * sin (θ). The value of each item S is cumulative to get Huang move, Unit is 0.0001 ". The calculation method of the angle nutation is similar to that of the Huang motion, and the value of the angular θ is the same, except that the cosine factor is used in the calculation of nutation:

C = (C1 + C2 * T) * cos (theta) (3.16-type)

The Calcearthlongitudenutation () function is the implementation code that calculates the Huang motion:

Double calcearthlongitudenutation (double dt)   
       
{   
       
    double T = dt *;   
       
    Double D,m,mp,f,omega;   
       
       
       
    Getearthnutationparameter (DT, &d, &m, &MP, &f, &omega);   
       
       
       
    Double resulte = 0.0;   
       
    for (int i = 0; i < sizeof (nutation)/sizeof (nutation[0)); i++)   
       
    {   
       
        double sita = nutation[i]. D * d + nutation[i]. M * m + nutation[i]. MP * MP + 

Nutation[i]. F * f + nutation[i].omega * Omega;   
       
               
       
        Resulte + = (nutation[i].sine1 + nutation[i].sine2 * T) * sin (SITA);   
       
    }   
       
       
       
    /* First times the factor of the nutation table 0.0001, and then the unit of conversion into degrees
       
    /return resulte * 0.0001/3600.0;   
       
}

The Getearthnutationparameter () auxiliary function is used to calculate 5 basic angular distances:

void 

Getearthnutationparameter (double dt, double *d, double *m, double *mp, double *f, double 

*omega)   
       
{   
       
    Double T = dt * 10; /*t is the number of centuries from the J2000 of Confucianism
       
    /Double T2 = t * t;   
       
    Double T3 = T2 * T;   
       
       
       
    /* Flat pitch angle (such as the angular distance of the Earth's core) * * *
       
    *d = 297.85036 + 445267.111480 * T-0.0019142 * T2 + t3/189474.0;   
       
       
       
    * * SUN (Earth) flat near point angle * *
       
    *m = 357.52772 + 35999.050340 * T-0.0001603 * t2-t3/300000.0;   
       
       
       
    * * Moon Flat near Point angle * *
       
    *MP = 134.96298 + 477198.867398 * T + 0.0086972 * T2 + t3/56250.0;   
       
       
       
    /* Moon Latitude Parameter * *
       
    *f = 93.27191 + 483202.017538 * T-0.0036825 * T2 + t3/327270.0;   
       
       
       
    * * The ecliptic and the moon orbit rise the intersection of the Yellow meridian * * *
       
    *omega = 125.04452-1934.136261 * T + 0.0020708 * T2 + t3/450000.0;   
       
}