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