1*59cc4ca5SDavid du Colombier #include "astro.h"
2*59cc4ca5SDavid du Colombier
3*59cc4ca5SDavid du Colombier static double elem[] =
4*59cc4ca5SDavid du Colombier {
5*59cc4ca5SDavid du Colombier 36525.0, // [0] eday of epoc
6*59cc4ca5SDavid du Colombier
7*59cc4ca5SDavid du Colombier 39.48168677, // [1] semi major axis (au)
8*59cc4ca5SDavid du Colombier 0.24880766, // [2] eccentricity
9*59cc4ca5SDavid du Colombier 17.14175, // [3] inclination (deg)
10*59cc4ca5SDavid du Colombier 110.30347, // [4] longitude of the ascending node (deg)
11*59cc4ca5SDavid du Colombier 224.06676, // [5] longitude of perihelion (deg)
12*59cc4ca5SDavid du Colombier 238.92881, // [6] mean longitude (deg)
13*59cc4ca5SDavid du Colombier
14*59cc4ca5SDavid du Colombier -0.00076912, // [1+6] (au/julian century)
15*59cc4ca5SDavid du Colombier 0.00006465, // [2+6] (e/julian century)
16*59cc4ca5SDavid du Colombier 11.07, // [3+6] (arcsec/julian century)
17*59cc4ca5SDavid du Colombier -37.33, // [4+6] (arcsec/julian century)
18*59cc4ca5SDavid du Colombier -132.25, // [5+6] (arcsec/julian century)
19*59cc4ca5SDavid du Colombier 522747.90, // [6+6] (arcsec/julian century)
20*59cc4ca5SDavid du Colombier };
21*59cc4ca5SDavid du Colombier
22*59cc4ca5SDavid du Colombier void
plut(void)23*59cc4ca5SDavid du Colombier plut(void)
24*59cc4ca5SDavid du Colombier {
25*59cc4ca5SDavid du Colombier double pturbl, pturbb, pturbr;
26*59cc4ca5SDavid du Colombier double lograd;
27*59cc4ca5SDavid du Colombier double dele, enom, vnom, nd, sl;
28*59cc4ca5SDavid du Colombier
29*59cc4ca5SDavid du Colombier double capj, capn, eye, comg, omg;
30*59cc4ca5SDavid du Colombier double sb, su, cu, u, b, up;
31*59cc4ca5SDavid du Colombier double sd, ca, sa;
32*59cc4ca5SDavid du Colombier
33*59cc4ca5SDavid du Colombier double cy;
34*59cc4ca5SDavid du Colombier
35*59cc4ca5SDavid du Colombier cy = (eday - elem[0]) / 36525.; // per julian century
36*59cc4ca5SDavid du Colombier
37*59cc4ca5SDavid du Colombier mrad = elem[1] + elem[1+6]*cy;
38*59cc4ca5SDavid du Colombier ecc = elem[2] + elem[2+6]*cy;
39*59cc4ca5SDavid du Colombier
40*59cc4ca5SDavid du Colombier cy = cy / 3600; // arcsec/deg per julian century
41*59cc4ca5SDavid du Colombier incl = elem[3] + elem[3+6]*cy;
42*59cc4ca5SDavid du Colombier node = elem[4] + elem[4+6]*cy;
43*59cc4ca5SDavid du Colombier argp = elem[5] + elem[5+6]*cy;
44*59cc4ca5SDavid du Colombier
45*59cc4ca5SDavid du Colombier anom = elem[6] + elem[6+6]*cy - argp;
46*59cc4ca5SDavid du Colombier motion = elem[6+6] / 36525. / 3600;
47*59cc4ca5SDavid du Colombier
48*59cc4ca5SDavid du Colombier incl *= radian;
49*59cc4ca5SDavid du Colombier node *= radian;
50*59cc4ca5SDavid du Colombier argp *= radian;
51*59cc4ca5SDavid du Colombier anom = fmod(anom,360.)*radian;
52*59cc4ca5SDavid du Colombier
53*59cc4ca5SDavid du Colombier enom = anom + ecc*sin(anom);
54*59cc4ca5SDavid du Colombier do {
55*59cc4ca5SDavid du Colombier dele = (anom - enom + ecc * sin(enom)) /
56*59cc4ca5SDavid du Colombier (1. - ecc*cos(enom));
57*59cc4ca5SDavid du Colombier enom += dele;
58*59cc4ca5SDavid du Colombier } while(fabs(dele) > converge);
59*59cc4ca5SDavid du Colombier vnom = 2.*atan2(sqrt((1.+ecc)/(1.-ecc))*sin(enom/2.),
60*59cc4ca5SDavid du Colombier cos(enom/2.));
61*59cc4ca5SDavid du Colombier rad = mrad*(1. - ecc*cos(enom));
62*59cc4ca5SDavid du Colombier
63*59cc4ca5SDavid du Colombier lambda = vnom + argp;
64*59cc4ca5SDavid du Colombier pturbl = 0.;
65*59cc4ca5SDavid du Colombier lambda += pturbl*radsec;
66*59cc4ca5SDavid du Colombier pturbb = 0.;
67*59cc4ca5SDavid du Colombier pturbr = 0.;
68*59cc4ca5SDavid du Colombier
69*59cc4ca5SDavid du Colombier /*
70*59cc4ca5SDavid du Colombier * reduce to the ecliptic
71*59cc4ca5SDavid du Colombier */
72*59cc4ca5SDavid du Colombier
73*59cc4ca5SDavid du Colombier nd = lambda - node;
74*59cc4ca5SDavid du Colombier lambda = node + atan2(sin(nd)*cos(incl),cos(nd));
75*59cc4ca5SDavid du Colombier
76*59cc4ca5SDavid du Colombier sl = sin(incl)*sin(nd) + pturbb*radsec;
77*59cc4ca5SDavid du Colombier beta = atan2(sl, pyth(sl));
78*59cc4ca5SDavid du Colombier
79*59cc4ca5SDavid du Colombier lograd = pturbr*2.30258509;
80*59cc4ca5SDavid du Colombier rad *= 1. + lograd;
81*59cc4ca5SDavid du Colombier
82*59cc4ca5SDavid du Colombier
83*59cc4ca5SDavid du Colombier lambda -= 1185.*radsec;
84*59cc4ca5SDavid du Colombier beta -= 51.*radsec;
85*59cc4ca5SDavid du Colombier
86*59cc4ca5SDavid du Colombier motion *= radian*mrad*mrad/(rad*rad);
87*59cc4ca5SDavid du Colombier semi = 83.33;
88*59cc4ca5SDavid du Colombier
89*59cc4ca5SDavid du Colombier /*
90*59cc4ca5SDavid du Colombier * here begins the computation of magnitude
91*59cc4ca5SDavid du Colombier * first find the geocentric equatorial coordinates of Saturn
92*59cc4ca5SDavid du Colombier */
93*59cc4ca5SDavid du Colombier
94*59cc4ca5SDavid du Colombier sd = rad*(cos(beta)*sin(lambda)*sin(obliq) +
95*59cc4ca5SDavid du Colombier sin(beta)*cos(obliq));
96*59cc4ca5SDavid du Colombier sa = rad*(cos(beta)*sin(lambda)*cos(obliq) -
97*59cc4ca5SDavid du Colombier sin(beta)*sin(obliq));
98*59cc4ca5SDavid du Colombier ca = rad*cos(beta)*cos(lambda);
99*59cc4ca5SDavid du Colombier sd += zms;
100*59cc4ca5SDavid du Colombier sa += yms;
101*59cc4ca5SDavid du Colombier ca += xms;
102*59cc4ca5SDavid du Colombier alpha = atan2(sa,ca);
103*59cc4ca5SDavid du Colombier delta = atan2(sd,sqrt(sa*sa+ca*ca));
104*59cc4ca5SDavid du Colombier
105*59cc4ca5SDavid du Colombier /*
106*59cc4ca5SDavid du Colombier * here are the necessary elements of Saturn's rings
107*59cc4ca5SDavid du Colombier * cf. Exp. Supp. p. 363ff.
108*59cc4ca5SDavid du Colombier */
109*59cc4ca5SDavid du Colombier
110*59cc4ca5SDavid du Colombier capj = 6.9056 - 0.4322*capt;
111*59cc4ca5SDavid du Colombier capn = 126.3615 + 3.9894*capt + 0.2403*capt2;
112*59cc4ca5SDavid du Colombier eye = 28.0743 - 0.0128*capt;
113*59cc4ca5SDavid du Colombier comg = 168.1179 + 1.3936*capt;
114*59cc4ca5SDavid du Colombier omg = 42.9236 - 2.7390*capt - 0.2344*capt2;
115*59cc4ca5SDavid du Colombier
116*59cc4ca5SDavid du Colombier capj *= radian;
117*59cc4ca5SDavid du Colombier capn *= radian;
118*59cc4ca5SDavid du Colombier eye *= radian;
119*59cc4ca5SDavid du Colombier comg *= radian;
120*59cc4ca5SDavid du Colombier omg *= radian;
121*59cc4ca5SDavid du Colombier
122*59cc4ca5SDavid du Colombier /*
123*59cc4ca5SDavid du Colombier * now find saturnicentric ring-plane coords of the earth
124*59cc4ca5SDavid du Colombier */
125*59cc4ca5SDavid du Colombier
126*59cc4ca5SDavid du Colombier sb = sin(capj)*cos(delta)*sin(alpha-capn) -
127*59cc4ca5SDavid du Colombier cos(capj)*sin(delta);
128*59cc4ca5SDavid du Colombier su = cos(capj)*cos(delta)*sin(alpha-capn) +
129*59cc4ca5SDavid du Colombier sin(capj)*sin(delta);
130*59cc4ca5SDavid du Colombier cu = cos(delta)*cos(alpha-capn);
131*59cc4ca5SDavid du Colombier u = atan2(su,cu);
132*59cc4ca5SDavid du Colombier b = atan2(sb,sqrt(su*su+cu*cu));
133*59cc4ca5SDavid du Colombier
134*59cc4ca5SDavid du Colombier /*
135*59cc4ca5SDavid du Colombier * and then the saturnicentric ring-plane coords of the sun
136*59cc4ca5SDavid du Colombier */
137*59cc4ca5SDavid du Colombier
138*59cc4ca5SDavid du Colombier su = sin(eye)*sin(beta) +
139*59cc4ca5SDavid du Colombier cos(eye)*cos(beta)*sin(lambda-comg);
140*59cc4ca5SDavid du Colombier cu = cos(beta)*cos(lambda-comg);
141*59cc4ca5SDavid du Colombier up = atan2(su,cu);
142*59cc4ca5SDavid du Colombier
143*59cc4ca5SDavid du Colombier /*
144*59cc4ca5SDavid du Colombier * at last, the magnitude
145*59cc4ca5SDavid du Colombier */
146*59cc4ca5SDavid du Colombier
147*59cc4ca5SDavid du Colombier
148*59cc4ca5SDavid du Colombier sb = sin(b);
149*59cc4ca5SDavid du Colombier mag = -8.68 +2.52*fabs(up+omg-u)-
150*59cc4ca5SDavid du Colombier 2.60*fabs(sb) + 1.25*(sb*sb);
151*59cc4ca5SDavid du Colombier
152*59cc4ca5SDavid du Colombier helio();
153*59cc4ca5SDavid du Colombier geo();
154*59cc4ca5SDavid du Colombier }
155