1 #include "astro.h"
2
3 void
sat(void)4 sat(void)
5 {
6 double pturbl, pturbb, pturbr;
7 double lograd;
8 double dele, enom, vnom, nd, sl;
9
10 double capj, capn, eye, comg, omg;
11 double sb, su, cu, u, b, up;
12 double sd, ca, sa;
13
14 ecc = .0558900 - .000347*capt;
15 incl = 2.49256 - .0044*capt;
16 node = 112.78364 + .87306*capt;
17 argp = 91.08897 + 1.95917*capt;
18 mrad = 9.538843;
19 anom = 175.47630 + .03345972*eday - .56527*capt;
20 motion = 120.4550/3600.;
21
22 incl *= radian;
23 node *= radian;
24 argp *= radian;
25 anom = fmod(anom, 360.)*radian;
26
27 enom = anom + ecc*sin(anom);
28 do {
29 dele = (anom - enom + ecc * sin(enom)) /
30 (1. - ecc*cos(enom));
31 enom += dele;
32 } while(fabs(dele) > converge);
33 vnom = 2.*atan2(sqrt((1.+ecc)/(1.-ecc))*sin(enom/2.),
34 cos(enom/2.));
35 rad = mrad*(1. - ecc*cos(enom));
36
37 lambda = vnom + argp;
38 pturbl = 0.;
39 lambda += pturbl*radsec;
40 pturbb = 0.;
41 pturbr = 0.;
42
43 /*
44 * reduce to the ecliptic
45 */
46
47 nd = lambda - node;
48 lambda = node + atan2(sin(nd)*cos(incl),cos(nd));
49
50 sl = sin(incl)*sin(nd) + pturbb*radsec;
51 beta = atan2(sl, pyth(sl));
52
53 lograd = pturbr*2.30258509;
54 rad *= 1. + lograd;
55
56
57 lambda -= 1185.*radsec;
58 beta -= 51.*radsec;
59
60 motion *= radian*mrad*mrad/(rad*rad);
61 semi = 83.33;
62
63 /*
64 * here begins the computation of magnitude
65 * first find the geocentric equatorial coordinates of Saturn
66 */
67
68 sd = rad*(cos(beta)*sin(lambda)*sin(obliq) +
69 sin(beta)*cos(obliq));
70 sa = rad*(cos(beta)*sin(lambda)*cos(obliq) -
71 sin(beta)*sin(obliq));
72 ca = rad*cos(beta)*cos(lambda);
73 sd += zms;
74 sa += yms;
75 ca += xms;
76 alpha = atan2(sa,ca);
77 delta = atan2(sd,sqrt(sa*sa+ca*ca));
78
79 /*
80 * here are the necessary elements of Saturn's rings
81 * cf. Exp. Supp. p. 363ff.
82 */
83
84 capj = 6.9056 - 0.4322*capt;
85 capn = 126.3615 + 3.9894*capt + 0.2403*capt2;
86 eye = 28.0743 - 0.0128*capt;
87 comg = 168.1179 + 1.3936*capt;
88 omg = 42.9236 - 2.7390*capt - 0.2344*capt2;
89
90 capj *= radian;
91 capn *= radian;
92 eye *= radian;
93 comg *= radian;
94 omg *= radian;
95
96 /*
97 * now find saturnicentric ring-plane coords of the earth
98 */
99
100 sb = sin(capj)*cos(delta)*sin(alpha-capn) -
101 cos(capj)*sin(delta);
102 su = cos(capj)*cos(delta)*sin(alpha-capn) +
103 sin(capj)*sin(delta);
104 cu = cos(delta)*cos(alpha-capn);
105 u = atan2(su,cu);
106 b = atan2(sb,sqrt(su*su+cu*cu));
107
108 /*
109 * and then the saturnicentric ring-plane coords of the sun
110 */
111
112 su = sin(eye)*sin(beta) +
113 cos(eye)*cos(beta)*sin(lambda-comg);
114 cu = cos(beta)*cos(lambda-comg);
115 up = atan2(su,cu);
116
117 /*
118 * at last, the magnitude
119 */
120
121
122 sb = sin(b);
123 mag = -8.68 +2.52*fabs(up+omg-u)-
124 2.60*fabs(sb) + 1.25*(sb*sb);
125
126 helio();
127 geo();
128 }
129