1 #include "astro.h"
2
3 void
merc(void)4 merc(void)
5 {
6 double pturbl, pturbr;
7 double lograd;
8 double dele, enom, vnom, nd, sl;
9 double q0, v0, t0, j0 , s0;
10 double lsun, elong, ci, dlong;
11
12 ecc = .20561421 + .00002046*capt - 0.03e-6*capt2;
13 incl = 7.0028806 + .0018608*capt - 18.3e-6*capt2;
14 node = 47.145944 + 1.185208*capt + .0001739*capt2;
15 argp = 75.899697 + 1.555490*capt + .0002947*capt2;
16 mrad = .3870986;
17 anom = 102.279381 + 4.0923344364*eday + 6.7e-6*capt2;
18 motion = 4.0923770233;
19
20 q0 = 102.28 + 4.092334429*eday;
21 v0 = 212.536 + 1.602126105*eday;
22 t0 = -1.45 + .985604737*eday;
23 j0 = 225.36 + .083086735*eday;
24 s0 = 175.68 + .033455441*eday;
25
26 q0 *= radian;
27 v0 *= radian;
28 t0 *= radian;
29 j0 *= radian;
30 s0 *= radian;
31
32 incl *= radian;
33 node *= radian;
34 argp *= radian;
35 anom = fmod(anom, 360.)*radian;
36
37
38 enom = anom + ecc*sin(anom);
39 do {
40 dele = (anom - enom + ecc * sin(enom)) /
41 (1. - ecc*cos(enom));
42 enom += dele;
43 } while(fabs(dele) > converge);
44 vnom = 2.*atan2(sqrt((1.+ecc)/(1.-ecc))*sin(enom/2.),
45 cos(enom/2.));
46 rad = mrad*(1. - ecc*cos(enom));
47
48 icosadd(mercfp, merccp);
49 pturbl = cosadd(2, q0, -v0);
50 pturbl += cosadd(2, q0, -t0);
51 pturbl += cosadd(2, q0, -j0);
52 pturbl += cosadd(2, q0, -s0);
53
54 pturbr = cosadd(2, q0, -v0);
55 pturbr += cosadd(2, q0, -t0);
56 pturbr += cosadd(2, q0, -j0);
57
58 /*
59 * reduce to the ecliptic
60 */
61
62 lambda = vnom + argp + pturbl*radsec;
63 nd = lambda - node;
64 lambda = node + atan2(sin(nd)*cos(incl), cos(nd));
65
66 sl = sin(incl)*sin(nd);
67 beta = atan2(sl, pyth(sl));
68
69 lograd = pturbr*2.30258509;
70 rad *= 1. + lograd;
71
72 motion *= radian*mrad*mrad/(rad*rad);
73 semi = 3.34;
74
75 lsun = 99.696678 + 0.9856473354*eday;
76 lsun *= radian;
77 elong = lambda - lsun;
78 ci = (rad - cos(elong))/sqrt(1. + rad*rad - 2.*rad*cos(elong));
79 dlong = atan2(pyth(ci), ci)/radian;
80 mag = -.003 + .01815*dlong + .0001023*dlong*dlong;
81
82 helio();
83 geo();
84 }
85