123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135 |
- /*
- * This file is part of the UCB release of Plan 9. It is subject to the license
- * terms in the LICENSE file found in the top-level directory of this
- * distribution and at http://akaros.cs.berkeley.edu/files/Plan9License. No
- * part of the UCB release of Plan 9, including this file, may be copied,
- * modified, propagated, or distributed except according to the terms contained
- * in the LICENSE file.
- */
- #include "astro.h"
- void
- venus(void)
- {
- double pturbl, pturbb, pturbr;
- double lograd;
- double dele, enom, vnom, nd, sl;
- double v0, t0, m0, j0, s0;
- double lsun, elong, ci, dlong;
- /*
- * here are the mean orbital elements
- */
- ecc = .00682069 - .00004774*capt + 0.091e-6*capt2;
- incl = 3.393631 + .0010058*capt - 0.97e-6*capt2;
- node = 75.779647 + .89985*capt + .00041*capt2;
- argp = 130.163833 + 1.408036*capt - .0009763*capt2;
- mrad = .7233316;
- anom = 212.603219 + 1.6021301540*eday + .00128605*capt2;
- motion = 1.6021687039;
- /*
- * mean anomalies of perturbing planets
- */
- v0 = 212.60 + 1.602130154*eday;
- t0 = 358.63 + .985608747*eday;
- m0 = 319.74 + 0.524032490*eday;
- j0 = 225.43 + .083090842*eday;
- s0 = 175.8 + .033459258*eday;
- v0 *= radian;
- t0 *= radian;
- m0 *= radian;
- j0 *= radian;
- s0 *= radian;
- incl *= radian;
- node *= radian;
- argp *= radian;
- anom = fmod(anom, 360.)*radian;
- /*
- * computation of long period terms affecting the mean anomaly
- */
- anom +=
- (2.761-0.022*capt)*radsec*sin(
- 13.*t0 - 8.*v0 + 43.83*radian + 4.52*radian*capt)
- + 0.268*radsec*cos(4.*m0 - 7.*t0 + 3.*v0)
- + 0.019*radsec*sin(4.*m0 - 7.*t0 + 3.*v0)
- - 0.208*radsec*sin(s0 + 1.4*radian*capt);
- /*
- * computation of elliptic orbit
- */
- enom = anom + ecc*sin(anom);
- do {
- dele = (anom - enom + ecc * sin(enom)) /
- (1 - ecc*cos(enom));
- enom += dele;
- } while(fabs(dele) > converge);
- vnom = 2*atan2(sqrt((1+ecc)/(1-ecc))*sin(enom/2),
- cos(enom/2));
- rad = mrad*(1 - ecc*cos(enom));
- lambda = vnom + argp;
- /*
- * perturbations in longitude
- */
- icosadd(venfp, vencp);
- pturbl = cosadd(4, v0, t0, m0, j0);
- pturbl *= radsec;
- /*
- * perturbations in latidude
- */
- pturbb = cosadd(3, v0, t0, j0);
- pturbb *= radsec;
- /*
- * perturbations in log radius vector
- */
- pturbr = cosadd(4, v0, t0, m0, j0);
- /*
- * reduction to the ecliptic
- */
- lambda += pturbl;
- nd = lambda - node;
- lambda = node + atan2(sin(nd)*cos(incl),cos(nd));
- sl = sin(incl)*sin(nd);
- beta = atan2(sl, pyth(sl)) + pturbb;
- lograd = pturbr*2.30258509;
- rad *= 1 + lograd;
- motion *= radian*mrad*mrad/(rad*rad);
- /*
- * computation of magnitude
- */
- lsun = 99.696678 + 0.9856473354*eday;
- lsun *= radian;
- elong = lambda - lsun;
- ci = (rad - cos(elong))/sqrt(1 + rad*rad - 2*rad*cos(elong));
- dlong = atan2(pyth(ci), ci)/radian;
- mag = -4 + .01322*dlong + .0000004247*dlong*dlong*dlong;
- semi = 8.41;
- helio();
- geo();
- }
|