![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to sun.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 sun.c CVS log | Up to [Research Unix] / researchv10no / cmd / sky |
1.1 root 1: #include "sky.h"
2:
3: extern struct suntab
4: {
5: float f[2];
6: char c[3];
7: } suntab[];
8: /*
9: * WARNING: Requires sign extension of characters.
10: */
11:
12: sun()
13: {
14: double mmerc, mven, merth, mmars, mjup, msat;
15: double dmoon, mmoon, gmoon;
16: double pturbb, pturbl, pturbr, lograd;
17: double planp[7];
18: struct suntab *pp;
19: double Eday, capT, capT2, capT3;
20: pp = &suntab[0];
21: Eday = eday - 36525.;
22: capT = Eday/36525.;
23: capT2 = capT*capT;
24: capT3 = capT2*capT;
25:
26: /*
27: * The arugments and coefficients are taken from
28: * Simon Newcomb, "Tables of the Motion of the Earth upon
29: * its Axis and About the Sun", A.P.A.E. VI (1895)
30: */
31:
32: object = "sun ";
33: incl = 0.;
34: ecc = .01670911 - 4.205e-5 * capT - 1.26e-7*capT2;
35: node = 0.;
36: argp = 282.938352 + .0000470774*Eday + .000462*capT2
37: + .000003*capT3;
38: mrad = 1.0;
39: anom = 357.527723 + .985600283*Eday - .000159*capT2
40: - .000003*capT3;
41: motion = .9856473354;
42:
43: argp -= 0.00085; /*empirical*/
44: anom += 0.00085; /*empirical*/
45:
46: dmoon = 350.737681+12.1907491914*eday-.001436*capt2;
47: gmoon = 11.250889 + 13.2293504490*eday - .003212*capt2;
48: mmoon = 296.104608 + 13.0649924465*eday + 9.192e-3*capt2;
49: mmerc = 102.19 + 4.092329959*eday + .19;
50: mven = 212.448 + 1.602121635*eday + .17;
51: merth = 358.476 + 0.985600267*eday + .19;
52: mmars = 319.590 + .524024095*eday + .11;
53: mjup = 225.269 + .083082362*eday + .34;
54: msat = 175.593 + .033450794*eday + .24;
55: /* muran = 72.248 + .011722663*eday + .19;*/
56:
57: dmoon = fmod(dmoon, 360.)*radian;
58: gmoon = fmod(gmoon, 360.)*radian;
59: mmoon = fmod(mmoon, 360.)*radian;
60: mmerc *= radian;
61: mven *= radian;
62: merth *= radian;
63: mmars *= radian;
64: mjup *= radian;
65: msat *= radian;
66:
67:
68: /*
69: long period terms in the mean anomaly - the arguments
70: of the terms are:
71: 4*mars - 7*earth + 3*venus
72: 3*jup - 8*mars + 4*earth
73: 13*earth - 8*venus
74: 8*earth - 15*mars
75: */
76:
77: anom = anom
78: + 0.266/3600.*sin((31.8 + 119.0*capt)*radian)
79: + 6.40 /3600.*sin((231.19 + 20.20*capt)*radian)
80: + 1.870/3600.*sin((57.24+150.27*capt)*radian);
81:
82: incl *= radian;
83: node *= radian;
84: argp *= radian;
85: anom = fmod(anom, 360.)*radian;
86:
87: /*
88: * computation of elliptic orbit
89: */
90:
91: lambda = anom + argp;
92:
93: pturbl = (2.*ecc - ecc*ecc*ecc/4.)*sin(anom)
94: + (ecc*ecc*5./4. - ecc*ecc*ecc*ecc*11./24.)*sin(2.*anom)
95: + (ecc*ecc*ecc*13./12.)*sin(3.*anom)
96: + (ecc*ecc*ecc*ecc*103./96.)*sin(4.*anom);
97:
98: pturbl -= 0.27*radsec*sin(anom);
99:
100: lambda += pturbl;
101:
102: beta = 0.;
103:
104: lograd = (ecc*ecc/4. + ecc*ecc*ecc*ecc/32.)
105: + (-ecc + ecc*ecc*ecc*3./8.)*cos(anom)
106: + (-ecc*ecc*3./4. + ecc*ecc*ecc*ecc*11./24.)*cos(2.*anom)
107: + (-ecc*ecc*ecc*17./24.)*cos(3.*anom)
108: + (-ecc*ecc*ecc*ecc*71./96.)*cos(4.*anom);
109:
110: lograd += 0.00000035; /*empirical*/
111: lograd += 0.00000070*cos(anom); /*empirical*/
112:
113: lograd = lograd/2.30258509;
114:
115: /*
116: * Computation of perturbations to elliptic orbit.
117: */
118:
119: planp[1] = mmerc;
120: planp[2] = mven;
121: planp[3] = merth;
122: planp[4] = mmars;
123: planp[5] = mjup;
124: planp[6] = msat;
125:
126: pturbl = 0.;
127: for(;;){
128: if(pp->f[0]==0.){
129: pp++;
130: break;
131: }
132: pturbl += pp->f[0]*cos(pp->f[1] + pp->c[0]*merth + pp->c[1]*planp[pp->c[2]]);
133: pp++;
134: }
135:
136: pturbl = pturbl
137: + 6.454*sin(dmoon)
138: + 0.015*sin(dmoon) /* empirical */
139: + 0.013*sin(3.*dmoon)
140: + 0.177*sin(dmoon + mmoon)
141: - 0.424*sin(dmoon - mmoon)
142: + 0.039*sin(3.*dmoon - mmoon)
143: - 0.064*sin(dmoon + merth)
144: + 0.172*sin(dmoon - merth)
145: /*
146: - 0.013*sin(dmoon-mmoon-merth)
147: - 0.013*sin(2.*gmoon-2.*dmoon)
148: */
149: ;
150:
151: pturbl *= radsec;
152:
153:
154: pturbb = 0.;
155: for(;;){
156: if(pp->f[0]==0.){
157: pp++;
158: break;
159: }
160: pturbb += pp->f[0]*cos(pp->f[1] + pp->c[0]*merth + pp->c[1]*planp[pp->c[2]]);
161: pp++;
162: }
163:
164:
165: pturbb = pturbb
166: + 0.576*sin(gmoon)
167: + 0.016*sin(gmoon + mmoon)
168: - 0.047*sin(gmoon - mmoon)
169: + 0.021*sin(2.*dmoon - gmoon)
170: /*
171: + 0.005*sin(2.*dmoon-gmoon-mmoon)
172: + 0.005*sin(gmoon+merth)
173: + 0.005*sin(gmoon-merth)
174: */
175: ;
176:
177: pturbb *= radsec;
178:
179:
180: pturbr = 0.;
181: for(;;){
182: if(pp->f[0]==0.){
183: pp++;
184: break;
185: }
186: pturbr += pp->f[0]*cos(pp->f[1] + pp->c[0]*merth + pp->c[1]*planp[pp->c[2]]);
187: pp++;
188: }
189:
190:
191: pturbr = pturbr
192: +13.360e-6*cos(dmoon)
193: + 0.030e-6*cos(3.*dmoon)
194: + 0.370e-6*cos(dmoon + mmoon)
195: - 1.330e-6*cos(dmoon - mmoon)
196: + 0.080e-6*cos(3.*dmoon - mmoon)
197: - 0.140e-6*cos(dmoon + merth)
198: + 0.360e-6*cos(dmoon - merth)
199: /*
200: - 0.030e-6*cos(dmoon-mmoon-merth)
201: + 0.030e-6*cos(2.*gmoon-2.*dmoon)
202: */
203: ;
204: lambda += pturbl;
205: if(lambda > 2.*pi)
206: lambda -= 2.*pi;
207:
208: beta += pturbb;
209:
210: lograd = (lograd+pturbr) * 2.30258509;
211: rad = 1. + lograd * (1. + lograd * (.5 + lograd/6.));
212:
213: motion *= radian*mrad*mrad/(rad*rad);
214:
215: semi = 961.182;
216: if(flags & OCCULT)
217: semi = 959.63;
218: mag = -26.5;
219: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.