1*3e12c5d1SDavid du Colombier #include "astro.h"
2*3e12c5d1SDavid du Colombier
3*3e12c5d1SDavid du Colombier void
nutate(void)4*3e12c5d1SDavid du Colombier nutate(void)
5*3e12c5d1SDavid du Colombier {
6*3e12c5d1SDavid du Colombier
7*3e12c5d1SDavid du Colombier /*
8*3e12c5d1SDavid du Colombier * uses radian, radsec
9*3e12c5d1SDavid du Colombier * sets phi, eps, dphi, deps, obliq, gst, tobliq
10*3e12c5d1SDavid du Colombier */
11*3e12c5d1SDavid du Colombier
12*3e12c5d1SDavid du Colombier double msun, mnom, noded, dmoon;
13*3e12c5d1SDavid du Colombier
14*3e12c5d1SDavid du Colombier /*
15*3e12c5d1SDavid du Colombier * nutation of the equinoxes is a wobble of the pole
16*3e12c5d1SDavid du Colombier * of the earths rotation whose magnitude is about
17*3e12c5d1SDavid du Colombier * 9 seconds of arc and whose period is about 18.6 years.
18*3e12c5d1SDavid du Colombier *
19*3e12c5d1SDavid du Colombier * it depends upon the pull of the sun and moon on the
20*3e12c5d1SDavid du Colombier * equatorial bulge of the earth.
21*3e12c5d1SDavid du Colombier *
22*3e12c5d1SDavid du Colombier * phi and eps are the two angles which specify the
23*3e12c5d1SDavid du Colombier * true pole with respect to the mean pole.
24*3e12c5d1SDavid du Colombier *
25*3e12c5d1SDavid du Colombier * all coeffieients are from Exp. Supp. pp.44-45
26*3e12c5d1SDavid du Colombier */
27*3e12c5d1SDavid du Colombier
28*3e12c5d1SDavid du Colombier mnom = 296.104608 + 13.0649924465*eday + 9.192e-3*capt2
29*3e12c5d1SDavid du Colombier + 14.38e-6*capt3;
30*3e12c5d1SDavid du Colombier mnom *= radian;
31*3e12c5d1SDavid du Colombier msun = 358.475833 + .9856002669*eday - .150e-3*capt2
32*3e12c5d1SDavid du Colombier - 3.33e-6*capt3;
33*3e12c5d1SDavid du Colombier msun *= radian;
34*3e12c5d1SDavid du Colombier noded = 11.250889 + 13.2293504490*eday - 3.211e-3*capt2
35*3e12c5d1SDavid du Colombier - 0.33e-6*capt3;
36*3e12c5d1SDavid du Colombier noded *= radian;
37*3e12c5d1SDavid du Colombier dmoon = 350.737486 + 12.1907491914*eday - 1.436e-3*capt2
38*3e12c5d1SDavid du Colombier + 1.89e-6*capt3;
39*3e12c5d1SDavid du Colombier dmoon *= radian;
40*3e12c5d1SDavid du Colombier node = 259.183275 - .0529539222*eday + 2.078e-3*capt2
41*3e12c5d1SDavid du Colombier + 2.22e-6*capt3;
42*3e12c5d1SDavid du Colombier node *= radian;
43*3e12c5d1SDavid du Colombier
44*3e12c5d1SDavid du Colombier /*
45*3e12c5d1SDavid du Colombier * long period terms
46*3e12c5d1SDavid du Colombier */
47*3e12c5d1SDavid du Colombier
48*3e12c5d1SDavid du Colombier phi = 0.;
49*3e12c5d1SDavid du Colombier eps = 0.;
50*3e12c5d1SDavid du Colombier dphi = 0.;
51*3e12c5d1SDavid du Colombier deps = 0.;
52*3e12c5d1SDavid du Colombier
53*3e12c5d1SDavid du Colombier
54*3e12c5d1SDavid du Colombier icosadd(nutfp, nutcp);
55*3e12c5d1SDavid du Colombier phi = -(17.2327+.01737*capt)*sin(node);
56*3e12c5d1SDavid du Colombier phi += sinadd(4, node, noded, dmoon, msun);
57*3e12c5d1SDavid du Colombier
58*3e12c5d1SDavid du Colombier eps = cosadd(4, node, noded, dmoon, msun);
59*3e12c5d1SDavid du Colombier
60*3e12c5d1SDavid du Colombier /*
61*3e12c5d1SDavid du Colombier * short period terms
62*3e12c5d1SDavid du Colombier */
63*3e12c5d1SDavid du Colombier
64*3e12c5d1SDavid du Colombier
65*3e12c5d1SDavid du Colombier dphi = sinadd(4, node, noded, mnom, dmoon);
66*3e12c5d1SDavid du Colombier
67*3e12c5d1SDavid du Colombier deps = cosadd(3, node, noded, mnom);
68*3e12c5d1SDavid du Colombier
69*3e12c5d1SDavid du Colombier phi = (phi+dphi)*radsec;
70*3e12c5d1SDavid du Colombier eps = (eps+deps)*radsec;
71*3e12c5d1SDavid du Colombier dphi *= radsec;
72*3e12c5d1SDavid du Colombier deps *= radsec;
73*3e12c5d1SDavid du Colombier
74*3e12c5d1SDavid du Colombier obliq = 23.452294 - .0130125*capt - 1.64e-6*capt2
75*3e12c5d1SDavid du Colombier + 0.503e-6*capt3;
76*3e12c5d1SDavid du Colombier obliq *= radian;
77*3e12c5d1SDavid du Colombier tobliq = obliq + eps;
78*3e12c5d1SDavid du Colombier
79*3e12c5d1SDavid du Colombier gst = 99.690983 + 360.9856473354*eday + .000387*capt2;
80*3e12c5d1SDavid du Colombier gst -= 180.;
81*3e12c5d1SDavid du Colombier gst = fmod(gst, 360.);
82*3e12c5d1SDavid du Colombier if(gst < 0.)
83*3e12c5d1SDavid du Colombier gst += 360.;
84*3e12c5d1SDavid du Colombier gst *= radian;
85*3e12c5d1SDavid du Colombier gst += phi*cos(obliq);
86*3e12c5d1SDavid du Colombier }
87