![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to venus.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 venus.c CVS log | Up to [Research Unix] / researchv10no / cmd / sky |
1.1 root 1: #include "sky.h"
2:
3: extern struct venust
4: {
5: float f[2];
6: char c[3];
7: } venust[];
8:
9: venus()
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 venust *pp = &venust[0];
18: double olong;
19: double temp;
20: double temp1;
21:
22: /*
23: * The arguments nnd coefficients are taken from
24: * Simon Newcomb, Tables of the Heliocentric Motion
25: * of Venus
26: * A.P.A.E. VI, part 3 (1895).
27: *
28: * Here are the mean orbital elements.
29: */
30:
31: object = "Venus ";
32: ecc = .00682069 - .00004774*capt + 0.091e-6*capt2;
33: incl = 3.393631 + .0010058*capt - 0.97e-6*capt2;
34: node = 75.779647 + .89985*capt + .00041*capt2;
35: argp = 130.163833 + 1.408036*capt - .0009763*capt2;
36: mrad = .7233316;
37: anom = 212.603219 + 1.6021301540*eday + .00128605*capt2;
38: motion = 1.6021687039;
39:
40: incl *= radian;
41: node *= radian;
42: argp *= radian;
43: anom = fmod(anom, 360.)*radian;
44: motion *= radian;
45:
46: /*
47: * Conventional mean anomalies of perturbing planets.
48: */
49:
50: q0 = 102.35 + 4.092338439*eday;
51: v0 = 212.60 + 1.602130154*eday;
52: t0 = 358.63 + .985608747*eday;
53: m0 = 319.74 + 0.524032490*eday;
54: j0 = 225.43 + .083090842*eday;
55: s0 = 175.8 + .033459258*eday;
56:
57: q0 *= radian;
58: v0 *= radian;
59: t0 *= radian;
60: m0 *= radian;
61: j0 *= radian;
62: s0 *= radian;
63:
64: planp[1] = q0;
65: planp[2] = v0;
66: planp[3] = t0;
67: planp[4] = m0;
68: planp[5] = j0;
69: planp[6] = s0;
70:
71: /*
72: * Computation of long period terms affecting the mean anomaly.
73: * 13*earth - 8.*venus
74: * 4*mars - 7.*earth + 3.*venus
75: * saturn
76: */
77:
78: anom +=
79: (2.761-0.022*capt)*radsec*sin((237.24+150.27*capt)*radian)
80: + 0.269*radsec*sin((212.2+119.05*capt)*radian)
81: - 0.208*radsec*sin((175.8+1223.5*capt)*radian);
82:
83: /*
84: * Computation of elliptic orbit.
85: */
86:
87: enom = anom + ecc*sin(anom);
88: do {
89: dele = (anom - enom + ecc * sin(enom)) /
90: (1. - ecc*cos(enom));
91: enom += dele;
92: } while(fabs(dele) > 1.e-8);
93: vnom = 2.*atan2(sqrt((1.+ecc)/(1.-ecc))*sin(enom/2.),
94: cos(enom/2.));
95: rad = mrad*(1. - ecc*cos(enom));
96:
97: /*
98: * Perturbations in longitude.
99: */
100:
101: pturbl = 0.;
102: for(;;){
103: if(pp->f[0]==0.){
104: pp++;
105: break;
106: }
107: pturbl += pp->f[0]*cos(pp->f[1] + pp->c[0]*v0 + pp->c[1]*planp[pp->c[2]]);
108: pp++;
109: }
110:
111: /*
112: * Perturbations in latitude.
113: */
114:
115: pturbb = 0.;
116: for(;;){
117: if(pp->f[0]==0.){
118: pp++;
119: break;
120: }
121: pturbb += pp->f[0]*cos(pp->f[1] + pp->c[0]*v0 + pp->c[1]*planp[pp->c[2]]);
122: pp++;
123: }
124:
125: /*
126: * Perturbations in log radius vector.
127: */
128:
129: pturbr = 0.;
130: for(;;){
131: if(pp->f[0]==0.){
132: pp++;
133: break;
134: }
135: pturbr += pp->f[0]*cos(pp->f[1] + pp->c[0]*v0 + pp->c[1]*planp[pp->c[2]]);
136: pp++;
137: }
138: pturbr *= 1.e-6;
139:
140: /*
141: * reduce to the ecliptic
142: */
143:
144: olong = vnom + argp + pturbl*radsec;
145: nd = olong - node;
146: lambda = node + atan2(sin(nd)*cos(incl), cos(nd));
147:
148: sl = sin(incl)*sin(nd);
149: beta = atan2(sl, sqrt(1.-sl*sl)) + pturbb*radsec;
150:
151: lograd = pturbr*2.30258509;
152: rad *= 1. + lograd;
153:
154: /*
155: * Compute motion for planetary aberration.
156: */
157:
158: temp = motion*mrad*mrad*sqrt(1.-ecc*ecc)/(rad*rad);
159: ldot = temp*sin(2.*(lambda-node))/sin(2.*(olong-node));
160: bdot = temp*sin(incl)*cos(lambda-node);
161: rdot = motion*mrad*ecc*sin(olong-argp)/sqrt(1.-ecc*ecc);
162:
163: /*
164: * Compute magnitude.
165: */
166:
167: lsun = 99.696678 + 0.9856473354*eday;
168: lsun *= radian;
169: elong = lambda - lsun;
170: ci = (rad - cos(elong))/sqrt(1. + rad*rad - 2.*rad*cos(elong));
171: dlong = atan2(sqrt(1.-ci*ci), ci)/radian;
172: mag = -4.00 + .01322*dlong + .0000004247*dlong*dlong*dlong;
173:
174: semi = 8.41;
175:
176: helio();
177: geo();
178:
179: /*
180: * transit computation
181: */
182:
183: if(!((flags&GEO)||(flags&HELIO))){
184: temp1 = sin(sundec)*sin(decl2) + cos(sundec)*cos(decl2)
185: *cos(sunra-ra);
186: temp1 = atan2(sqrt(1.-temp1*temp1),temp1)/radsec;
187: temp1 = temp1/(semi2+sunsd);
188: if(temp1 <= 1.0){
189: printf("Transit of Venus: Dist. = %.4f\n", temp1);
190: }else if(temp1 < 1.1){
191: printf("Near Transit of Venus: Dist. = %.4f\n", temp1);
192: }
193: }
194:
195: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.