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