venus.c (2417B)
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 }