1 #include "astro.h"
2 
3 
4 void
5 venus(void)
6 {
7 	double pturbl, pturbb, pturbr;
8 	double lograd;
9 	double dele, enom, vnom, nd, sl;
10 	double v0, t0, m0, j0, s0;
11 	double lsun, elong, ci, dlong;
12 
13 /*
14  *	here are the mean orbital elements
15  */
16 
17 	ecc = .00682069 - .00004774*capt + 0.091e-6*capt2;
18 	incl = 3.393631 + .0010058*capt - 0.97e-6*capt2;
19 	node = 75.779647 + .89985*capt + .00041*capt2;
20 	argp = 130.163833 + 1.408036*capt - .0009763*capt2;
21 	mrad = .7233316;
22 	anom = 212.603219 + 1.6021301540*eday + .00128605*capt2;
23 	motion = 1.6021687039;
24 
25 /*
26  *	mean anomalies of perturbing planets
27  */
28 
29 	v0 = 212.60 + 1.602130154*eday;
30 	t0 = 358.63  + .985608747*eday;
31 	m0 = 319.74 + 0.524032490*eday;
32 	j0 = 225.43 + .083090842*eday;
33 	s0 = 175.8  + .033459258*eday;
34 
35 	v0 *= radian;
36 	t0 *= radian;
37 	m0 *= radian;
38 	j0 *= radian;
39 	s0 *= radian;
40 
41 	incl *= radian;
42 	node *= radian;
43 	argp *= radian;
44 	anom = fmod(anom, 360.)*radian;
45 
46 /*
47  *	computation of long period terms affecting the mean anomaly
48  */
49 
50 	anom +=
51 		   (2.761-0.022*capt)*radsec*sin(
52 		  13.*t0 - 8.*v0 + 43.83*radian + 4.52*radian*capt)
53 		 + 0.268*radsec*cos(4.*m0 - 7.*t0 + 3.*v0)
54 		 + 0.019*radsec*sin(4.*m0 - 7.*t0 + 3.*v0)
55 		 - 0.208*radsec*sin(s0 + 1.4*radian*capt);
56 
57 /*
58  *	computation of elliptic orbit
59  */
60 
61 	enom = anom + ecc*sin(anom);
62 	do {
63 		dele = (anom - enom + ecc * sin(enom)) /
64 			(1 - ecc*cos(enom));
65 		enom += dele;
66 	} while(fabs(dele) > converge);
67 	vnom = 2*atan2(sqrt((1+ecc)/(1-ecc))*sin(enom/2),
68 		cos(enom/2));
69 	rad = mrad*(1 - ecc*cos(enom));
70 
71 	lambda = vnom + argp;
72 
73 /*
74  *	perturbations in longitude
75  */
76 
77 	icosadd(venfp, vencp);
78 	pturbl = cosadd(4, v0, t0, m0, j0);
79 	pturbl *= radsec;
80 
81 /*
82  *	perturbations in latidude
83  */
84 
85 	pturbb = cosadd(3, v0, t0, j0);
86 	pturbb *= radsec;
87 
88 /*
89  *	perturbations in log radius vector
90  */
91 
92 	pturbr = cosadd(4, v0, t0, m0, j0);
93 
94 /*
95  *	reduction to the ecliptic
96  */
97 
98 	lambda += pturbl;
99 	nd = lambda - node;
100 	lambda = node + atan2(sin(nd)*cos(incl),cos(nd));
101 
102 	sl = sin(incl)*sin(nd);
103 	beta = atan2(sl, pyth(sl)) + pturbb;
104 
105 	lograd = pturbr*2.30258509;
106 	rad *= 1 + lograd;
107 
108 
109 	motion *= radian*mrad*mrad/(rad*rad);
110 
111 /*
112  *	computation of magnitude
113  */
114 
115 	lsun = 99.696678 + 0.9856473354*eday;
116 	lsun *= radian;
117 	elong = lambda - lsun;
118 	ci = (rad - cos(elong))/sqrt(1 + rad*rad - 2*rad*cos(elong));
119 	dlong = atan2(pyth(ci), ci)/radian;
120 	mag = -4 + .01322*dlong + .0000004247*dlong*dlong*dlong;
121 
122 	semi = 8.41;
123 
124 	helio();
125 	geo();
126 }