![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to merc.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 merc.c CVS log | Up to [Research Unix] / researchv10no / cmd / sky |
1.1 root 1: #include "sky.h"
2:
3: extern struct merctab
4: {
5: float f[2];
6: char c[3];
7: } merctab[];
8:
9: merc()
10: {
11: double pturbl, pturbb, pturbr;
12: double lograd;
13: double dele, enom, vnom, nd, sl;
14: double q0, v0, t0, m0, j0 , s0;
15: double lsun, elong, ci, dlong;
16: double planp[7];
17: struct merctab *pp = &merctab[0];
18: double olong;
19: double temp;
20:
21: /*
22: * The arguments nnd coefficients are taken from
23: * Simon Newcomb, Tables of the Heliocentric Motion
24: * of Mercury
25: * A.P.A.E. VI, part 2 (1895).
26: *
27: * Here are the mean orbital elements.
28: */
29:
30: object = "Mercury ";
31: ecc = .20561421 + .00002046*capt - 0.03e-6*capt2;
32: ecc -= .00000250; /*empirical*/
33: ecc -= .00000003; /*empirical*/
34: incl = 7.0028806 + .0018608*capt - 18.3e-6*capt2;
35: incl += .00025;/* empirical */
36: node = 47.145944 + 1.185208*capt + .0001739*capt2;
37: node -= 2.2/3600.; /* empirical */
38: argp = 75.899697 + 1.555490*capt + .0002947*capt2;
39: argp += 0.06/3600. + 1.89/3600.; /* empirical */
40: mrad = .3870986;
41: anom = 102.279381 + 4.0923344364*eday + 6.7e-6*capt2;
42: anom += 2.07/3600. - 1.89/3600.; /* empirical */
43: motion = 4.0923770233;
44:
45: incl *= radian;
46: node *= radian;
47: argp *= radian;
48: anom = fmod(anom, 360.)*radian;
49: motion *= radian;
50:
51: /*
52: * Conventional mean anomalies of perturbing planets.
53: */
54:
55: q0 = 102.28 + 4.092334429*eday;
56: v0 = 212.536 + 1.602126105*eday;
57: t0 = -1.45 + .985604737*eday;
58: m0 = 319.66 + .524028480*eday;
59: j0 = 225.36 + .083086735*eday;
60: s0 = 175.68 + .033455441*eday;
61:
62: q0 *= radian;
63: v0 *= radian;
64: t0 *= radian;
65: m0 *= radian;
66: j0 *= radian;
67: s0 *= radian;
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:
76: /*
77: * Computation of long period terms affecting the mean anomaly.
78: */
79:
80: anom = anom;
81:
82: /*
83: * Computation of elliptic orbit.
84: */
85:
86: enom = anom + ecc*sin(anom);
87: do {
88: dele = (anom - enom + ecc * sin(enom)) /
89: (1. - ecc*cos(enom));
90: enom += dele;
91: } while(fabs(dele) > 1.e-8);
92: vnom = 2.*atan2(sqrt((1.+ecc)/(1.-ecc))*sin(enom/2.),
93: cos(enom/2.));
94: rad = mrad*(1. - ecc*cos(enom));
95:
96: /*
97: * Perturbations in longitude.
98: */
99:
100: pturbl = 0.;
101: for(;;){
102: if(pp->f[0]==0.){
103: pp++;
104: break;
105: }
106: pturbl += pp->f[0]*cos(pp->f[1] + pp->c[0]*q0 + pp->c[1]*planp[pp->c[2]]);
107: pp++;
108: }
109:
110: /*
111: * Perturbations in latitude.
112: */
113:
114: pturbb = 0.;
115: for(;;){
116: if(pp->f[0]==0.){
117: pp++;
118: break;
119: }
120: pturbb += pp->f[0]*cos(pp->f[1] + pp->c[0]*q0 + pp->c[1]*planp[pp->c[2]]);
121: pp++;
122: }
123:
124: /*
125: * Perturbations in log radius vector.
126: */
127:
128: pturbr = 0.;
129: for(;;){
130: if(pp->f[0]==0.){
131: pp++;
132: break;
133: }
134: pturbr += pp->f[0]*cos(pp->f[1] + pp->c[0]*q0 + pp->c[1]*planp[pp->c[2]]);
135: pp++;
136: }
137: pturbr *= 1.e-6;
138:
139: /*
140: * reduce to the ecliptic
141: */
142:
143: olong = vnom + argp + pturbl*radsec;
144: nd = olong - node;
145: lambda = node + atan2(sin(nd)*cos(incl), cos(nd));
146:
147: sl = sin(incl)*sin(nd);
148: beta = atan2(sl, sqrt(1.-sl*sl)) + pturbb*radsec;
149:
150: lograd = pturbr*2.30258509;
151: rad *= 1. + lograd;
152:
153: /*
154: * Compute motion for planetary aberration.
155: */
156:
157: temp = motion*mrad*mrad*sqrt(1.-ecc*ecc)/(rad*rad);
158: ldot = temp*sin(2.*(lambda-node))/sin(2.*(olong-node));
159: bdot = temp*sin(incl)*cos(lambda-node);
160: rdot = motion*mrad*ecc*sin(olong-argp)/sqrt(1.-ecc*ecc);
161:
162: /*
163: * Compute magnitude.
164: */
165:
166: lsun = 99.696678 + 0.9856473354*eday;
167: lsun *= radian;
168: elong = lambda - lsun;
169: ci = (rad - cos(elong))/sqrt(1. + rad*rad - 2.*rad*cos(elong));
170: dlong = atan2(sqrt(1.-ci*ci), ci)/radian;
171: mag = -.003 + .01815*dlong + .0001023*dlong*dlong;
172:
173: semi = 3.34;
174:
175: helio();
176: geo();
177:
178: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.