tnutate.c - plan9port - [fork] Plan 9 from user space
(HTM) git clone git://src.adamsgaard.dk/plan9port
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) README
(DIR) LICENSE
---
tnutate.c (1754B)
---
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 }