![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to mars.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 mars.c CVS log | Up to [Research Unix] / researchv10no / cmd / sky |
1.1 root 1: #include "sky.h"
2:
3: extern struct marst
4: {
5: float f[2];
6: char c[3];
7: } marst[];
8:
9: mars()
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;
15: double lsun, elong, ci, dlong;
16: double planp[8];
17: struct marst *pp = &marst[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 Mars
25: * A.P.A.E. VI, part 4 (1895).
26: *
27: * Here are the mean orbital elements.
28: */
29:
30: object = "Mars ";
31: ecc = .09331290 + .000092064*capt - 0.077e-6*capt2;
32: incl = 1.850333 - 6.75e-4*capt + 12.61e-6*capt2;
33: node = 48.786442 + .770992*capt - 1.39e-6*capt2
34: - 5.33e-6*capt3;
35: argp = 334.218203 + 1.840758*capt + 1.299e-4*capt2
36: - 1.19e-6*capt3;
37: mrad = 1.52368840;
38: anom = 319.529425 + .5240207666*eday + 1.808e-4*capt2
39: + 1.19e-6*capt3;
40: motion = 0.5240711638;
41:
42: incl *= radian;
43: node *= radian;
44: argp *= radian;
45: anom = fmod(anom, 360.)*radian;
46: motion *= radian;
47:
48: /*
49: * Conventional mean anomalies of perturbing planets.
50: */
51:
52: q0 = 102.28 + 4.092334429*eday;
53: v0 = 212.388 + 1.60211831*eday;
54: t0 = 358.415 + .98559696*eday;
55: m0 = 319.530 + .52402078*eday;
56: j0 = 225.209 + .08307904*eday + 0.332*sin((134.4+38.5*capt)*radian);
57: s0 = 175.533 + .03344747*eday - 0.808*sin((134.4+38.5*capt)*radian);
58: u0 = 74.188 + 0.0117193*eday;
59:
60: q0 *= radian;
61: v0 *= radian;
62: t0 *= radian;
63: m0 *= radian;
64: j0 *= radian;
65: s0 *= radian;
66: u0 *= radian;
67:
68: planp[1] = q0;
69: planp[2] = v0;
70: planp[3] = t0;
71: planp[4] = m0;
72: planp[5] = j0;
73: planp[6] = s0;
74: planp[7] = u0;
75:
76: /*
77: * Computation of long period terms affecting the mean anomaly.
78: * 4*mars - 7*earth + 3*venus
79: * 3*jupiter - 8*mars + 4*earth
80: * 2*jupiter - - 6*mars + 3*earth
81: * 2*saturn - 2*mars + earth
82: * jupiter - 2*mars + earth
83: * 5*saturn - 2*jupiter
84: */
85:
86: anom = anom
87: - (37.05 + 13.50*capt)*radsec
88: + 0.606*radsec*sin((212.87+119.051*capt)*radian)
89: + 52.490*radsec*sin((47.48+19.771*capt)*radian)
90: + 0.319*radsec*sin((116.88+773.444*capt)*radian)
91: + 0.130*radsec*sin((74.00+163.00*capt)*radian)
92: + 0.009*radsec*sin((325.00+753.67*capt)*radian)
93: + 0.280*radsec*sin((300.00+40.8*capt)*radian);
94:
95: /*
96: * Computation of elliptic orbit.
97: */
98:
99: enom = anom + ecc*sin(anom);
100: do {
101: dele = (anom - enom + ecc * sin(enom)) /
102: (1. - ecc*cos(enom));
103: enom += dele;
104: } while(fabs(dele) > 1.e-8);
105: vnom = 2.*atan2(sqrt((1.+ecc)/(1.-ecc))*sin(enom/2.),
106: cos(enom/2.));
107: rad = mrad*(1. - ecc*cos(enom));
108:
109: /*
110: * Perturbations in longitude.
111: */
112:
113: pturbl = 0.043*sin(2.*anom);
114: for(;;){
115: if(pp->f[0]==0.){
116: pp++;
117: break;
118: }
119: pturbl += pp->f[0]*cos(pp->f[1] + pp->c[0]*m0 + pp->c[1]*planp[pp->c[2]]);
120: pp++;
121: }
122:
123: /*
124: * Perturbations in latitude.
125: */
126:
127: pturbb = 0.;
128: for(;;){
129: if(pp->f[0]==0.){
130: pp++;
131: break;
132: }
133: pturbb += pp->f[0]*cos(pp->f[1] + pp->c[0]*m0 + pp->c[1]*planp[pp->c[2]]);
134: pp++;
135: }
136:
137: /*
138: * Perturbations in log radius vector.
139: */
140:
141: pturbr = 0.;
142: for(;;){
143: if(pp->f[0]==0.){
144: pp++;
145: break;
146: }
147: pturbr += pp->f[0]*cos(pp->f[1] + pp->c[0]*m0 + pp->c[1]*planp[pp->c[2]]);
148: pp++;
149: }
150: pturbr *= 1.e-6;
151:
152: /*
153: * reduce to the ecliptic
154: */
155:
156: olong = vnom + argp + pturbl*radsec;
157: nd = olong - node;
158: lambda = node + atan2(sin(nd)*cos(incl), cos(nd));
159:
160: sl = sin(incl)*sin(nd);
161: beta = atan2(sl, sqrt(1.-sl*sl)) + pturbb*radsec;
162:
163: lograd = pturbr*2.30258509;
164: rad *= 1. + lograd;
165:
166: /*
167: * Compute motion for planetary aberration.
168: */
169:
170: temp = motion*mrad*mrad*sqrt(1.-ecc*ecc)/(rad*rad);
171: ldot = temp*sin(2.*(lambda-node))/sin(2.*(olong-node));
172: bdot = temp*sin(incl)*cos(lambda-node);
173: rdot = motion*mrad*ecc*sin(olong-argp)/sqrt(1.-ecc*ecc);
174:
175: /*
176: * Compute magnitude.
177: */
178:
179: lsun = 99.696678 + 0.9856473354*eday;
180: lsun *= radian;
181: elong = lambda - lsun;
182: ci = (rad - cos(elong))/sqrt(1. + rad*rad - 2.*rad*cos(elong));
183: dlong = atan2(sqrt(1.-ci*ci), ci)/radian;
184: mag = -1.30 + .01486*dlong;
185:
186: semi = 4.68;
187:
188: helio();
189: geo();
190:
191: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.