![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to jup.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 jup.c CVS log | Up to [Research Unix] / researchv10no / cmd / sky |
1.1 root 1: #include "sky.h"
2:
3: extern struct juptab
4: {
5: float f[6];
6: char c[3];
7: } juptab[];
8:
9: jup()
10: {
11: double pturbl, pturbb, pturbr;
12: double lograd;
13: double dele, enom, vnom, nd, sl;
14: double q0, v0, t0, m0, j0 , s0, u0, n0;
15: double lsun, elong, ci, dlong;
16: double planp[9];
17: struct juptab *pp = &juptab[0];
18: double olong;
19: double grin;
20: double temp;
21:
22: /*
23: * The arguments nnd coefficients are taken from
24: *
25: * Here are the mean orbital elements.
26: */
27:
28: object = "Jupiter ";
29: ecc = 0.04833748 + 1.63903e-4*capt - 0.4642e-6*capt2 - 1.593e-9*capt3;
30: incl = 1.3086429 - 0.005696*capt + 3.89e-6*capt2;
31: node = 99.4431901 + 1.0105300*capt + 3.522e-4*capt2 - 8.51e-6*capt3;
32: argp = 13.4119487 + 0.21344495*capt + 7.466e-4*capt2 - 3.7946e-6*capt3;
33: anom = 225.3350378 + 0.0830853474*eday - 8.332e-4*capt2 + 3.9876e-6*capt3;
34: motion = .083091212;
35: mrad = 5.202803;
36:
37: incl *= radian;
38: node *= radian;
39: argp *= radian;
40: anom = fmod(anom, 360.)*radian;
41: motion *= radian;
42:
43: /*
44: * Longitudes of perturbing planets,
45: * They are of epoch Jan 0.5, 1900, and are
46: * referred to the fixed qquinox of that date.
47: */
48:
49: q0 = 178.179 + 4.092338817*eday;
50: v0 = 342.767 + 1.602130491*eday;
51: t0 = 99.697 + 0.985609114*eday;
52: m0 = 293.748 + 0.524032950*eday;
53: j0 = 238.050 + 0.083091230*eday;
54: s0 = 266.280 + 0.033459743*eday;
55: u0 = 243.370 + 0.011731421*eday;
56: n0 = 85.183 + 0.005987356*eday;
57:
58: q0 *= radian;
59: v0 *= radian;
60: t0 *= radian;
61: m0 *= radian;
62: j0 *= radian;
63: s0 *= radian;
64: u0 *= radian;
65: n0 *= radian;
66:
67: grin = 5.*s0 - 2.*j0 - 8079.0*radsec*capt;
68:
69: planp[1] = q0;
70: planp[2] = v0;
71: planp[3] = t0;
72: planp[4] = m0;
73: planp[5] = j0;
74: planp[6] = s0;
75: planp[7] = u0;
76: planp[8] = n0;
77:
78: /*
79: * Computation of long period terms affecting the mean anomaly.
80: */
81:
82: anom += 0.
83: +(1192.85-6.076*capt-0.0400*capt2+0.00075*capt3)*radsec*sin(grin)
84: +(-23.80-0.192*capt+0.0226*capt2-0.00080*capt3)*radsec*cos(grin)
85: +(-11.04-0.060*capt-0.0072*capt2+0.00021*capt3)*radsec*sin(2.*grin)
86: +(1.44-0.086*capt+0.0004*capt2+0.00006*capt3)*radsec*cos(2.*grin)
87:
88: + (8.22-0.120*capt-0.002*capt2)*radsec*sin(2.*j0-6.*s0+3.*u0)
89: + (0.55 + 0.420*capt - 0.0079*capt2)*radsec*cos(2.*j0-6.*s0+3.*u0)
90: ;
91:
92: /*
93: * Computation of elliptic orbit.
94: */
95:
96: enom = anom + ecc*sin(anom);
97: do {
98: dele = (anom - enom + ecc * sin(enom)) /
99: (1. - ecc*cos(enom));
100: enom += dele;
101: } while(fabs(dele) > 1.e-8);
102: vnom = 2.*atan2(sqrt((1.+ecc)/(1.-ecc))*sin(enom/2.),
103: cos(enom/2.));
104: rad = mrad*(1. - ecc*cos(enom));
105:
106: /*
107: * Perturbations in longitude.
108: */
109:
110: pturbl = 0.
111: +(-83.79-1.222*capt+0.0097*capt2)*sin(grin-j0)
112: +(137.08-1.508*capt-0.0069*capt2)*cos(grin-j0)
113: ;
114:
115: for(;;){
116: if(pp->c[2]==0){
117: pp++;
118: break;
119: }
120: temp = planp[pp->c[2]]*pp->c[0] + j0*pp->c[1];
121: pturbl += (pp->f[0]+pp->f[1]*capt+pp->f[2]*capt2)*sin(temp)
122: + (pp->f[3]+pp->f[4]*capt+pp->f[5]*capt2)*cos(temp);
123: pp++;
124: }
125:
126: /*
127: * Perturbations in latitude.
128: */
129:
130: pturbb = 0.;
131: /*
132: for(;;){
133: if(pp->f[0]==0.){
134: pp++;
135: break;
136: }
137: pturbb += pp->f[0]*cos(pp->f[1] + pp->c[0]*j0 + pp->c[1]*planp[pp->c[2]]);
138: pp++;
139: }
140: */
141:
142: /*
143: * Perturbations in log radius vector.
144: */
145:
146: pturbr = 0.;
147: /*
148: for(;;){
149: if(pp->f[0]==0.){
150: pp++;
151: break;
152: }
153: pturbr += pp->f[0]*cos(pp->f[1] + pp->c[0]*j0 + pp->c[1]*planp[pp->c[2]]);
154: pp++;
155: }
156: */
157: pturbr *= 1.e-6;
158:
159: /*
160: * reduce to the ecliptic
161: */
162:
163: olong = vnom + argp + pturbl*radsec;
164: nd = olong - node;
165: lambda = node + atan2(sin(nd)*cos(incl), cos(nd));
166:
167: sl = sin(incl)*sin(nd);
168: beta = atan2(sl, sqrt(1.-sl*sl)) + pturbb*radsec;
169:
170: lograd = pturbr*2.30258509;
171: rad *= 1. + lograd;
172:
173: /*
174: * Compute motion for planetary aberration.
175: */
176:
177: temp = motion*mrad*mrad*sqrt(1.-ecc*ecc)/(rad*rad);
178: ldot = temp*sin(2.*(lambda-node))/sin(2.*(olong-node));
179: bdot = temp*sin(incl)*cos(lambda-node);
180: rdot = motion*mrad*ecc*sin(olong-argp)/sqrt(1.-ecc*ecc);
181:
182: /*
183: * Compute magnitude.
184: */
185:
186: mag = -8.93;
187:
188: semi = 98.57;
189:
190: helio();
191: geo();
192:
193: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.