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