![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to nutate.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [Research Unix] / researchv10no / cmd / sky|
Return to nutate.c CVS log | Up to [Research Unix] / researchv10no / cmd / sky |
1.1 root 1: #include "sky.h"
2:
3: extern struct nuttab
4: {
5: float f[2];
6: char c[5];
7: } nuttab[];
8:
9: nutate()
10: {
11:
12: /*
13: * uses radian, radsec, eday
14: * sets psi, eps, dpsi, deps, obliq, gst, tobliq
15: */
16:
17: double msun, mnom, noded, dmoon;
18: struct nuttab *pp;
19: double arg;
20: double Eday;
21: double capT, capT2, capT3;
22:
23: /*
24: * nutation of the equinoxes is a irregular wobble of the
25: * pole of the earths rotation whose magnitude is about
26: * 9 seconds of arc and whose period is about 18.6 years.
27: *
28: * it depends upon the pull of the sun and moon on the
29: * equatorial bulge of the earth.
30: *
31: * psi and eps are the two angles which specify the
32: * true pole with respect to the mean pole.
33: *
34: * all coeffieients are from Supplement to A.E. 1984
35: */
36:
37: pp = &nuttab[0];
38: Eday = eday - 36525;
39: capT = capt-1;
40: capT2 = capT*capT;
41: capT3 = capT2*capT;
42:
43: mnom = 134.96298139 + 13.06499294724*Eday + 8.6972e-3*capT2
44: + 17.8e-6*capT3;
45: msun = 357.52772333 + 0.98560028309*Eday - 0.16028e-3*capT2
46: - 3.33e-6*capT3;
47: noded = 93.27191028 + 13.22935024060*Eday - 3.6825e-3*capT2
48: + 3.05e-6*capT3;
49: dmoon = 297.85036306 + 12.19074911650*Eday - 1.91417e-3*capT2
50: + 5.28e-6*capT3;
51: node = 125.04452222 - .05295376484*Eday + 2.07083e-3*capT2
52: + 2.22e-6*capT3;
53: mnom *= radian;
54: msun *= radian;
55: noded *= radian;
56: dmoon *= radian;
57: node *= radian;
58:
59: /*
60: * long period terms
61: */
62:
63:
64: psi = -(17.1996+.01742*capT)*sin(node);
65: eps = (9.2025+0.00089*capT)*cos(node);
66: for(;;){
67: if(pp->f[0]==0.){
68: pp++;
69: break;
70: }
71: arg = pp->c[0]*mnom + pp->c[1]*msun + pp->c[2]*noded
72: + pp->c[3]*dmoon + pp->c[4]*node;
73: psi += pp->f[0]*sin(arg);
74: eps += pp->f[1]*cos(arg);
75: pp++;
76: }
77:
78: /*
79: * short period terms
80: */
81:
82: dpsi = 0.;
83: deps = 0.;
84: if((flags&APPARENT)==0){
85: for(;;){
86: if(pp->f[0]==0.){
87: pp++;
88: break;
89: }
90: arg = pp->c[0]*mnom + pp->c[1]*msun + pp->c[2]*noded
91: + pp->c[3]*dmoon + pp->c[4]*node;
92: dpsi += pp->f[0]*sin(arg);
93: deps += pp->f[1]*cos(arg);
94: pp++;
95: }
96: }
97:
98:
99: psi = (psi+dpsi)*radsec;
100: eps = (eps+deps)*radsec;
101: dpsi *= radsec;
102: deps *= radsec;
103:
104: obliq = 23.43929111 - .013004167*capT - 0.164e-6*capT2
105: + 0.5036e-6*capT3;
106: obliq *= radian;
107: tobliq = obliq + eps;
108:
109: gst = 100.460618375 + 360.98564736629*Eday + 0.3882667e-3*capT2
110: - 0.025e-6*capT3;
111: gst -= 180.;
112: gst = fmod(gst, 360.);
113: if(gst < 0.)
114: gst += 360.;
115: gst *= radian;
116: gst += psi*cos(tobliq);
117:
118: /*
119: * print true obliquity, mean obliquity, mean GST, and true GST.
120: */
121:
122: if(flags&GOOBIE){
123: printf("Mean G. Sid. Time: ");
124: prhms(gst - psi*cos(obliq), 4);
125: printf("\n");
126:
127: printf("App. G. Sid. Time: ");
128: prhms(gst,4);
129: printf("\n");
130:
131: printf("Goobie: %6.3f %6.3f %6.3f %6.3f\n", psi/radsec,
132: eps/radsec, dpsi/radsec, deps/radsec);
133:
134: printf("Obliquity: ");
135: prdms(tobliq,3);
136: printf("\n");
137: }
138:
139: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.