1 cd5bae78 2004-04-21 devnull #include "astro.h"
3 cd5bae78 2004-04-21 devnull static double elem[] =
5 cbeb0b26 2006-04-01 devnull 36525.0, /* [0] eday of epoc */
7 cbeb0b26 2006-04-01 devnull 30.06896348, /* [1] semi major axis (au) */
8 cbeb0b26 2006-04-01 devnull 0.00858587, /* [2] eccentricity */
9 cbeb0b26 2006-04-01 devnull 1.76917, /* [3] inclination (deg) */
10 cbeb0b26 2006-04-01 devnull 131.72169, /* [4] longitude of the ascending node (deg) */
11 cbeb0b26 2006-04-01 devnull 44.97135, /* [5] longitude of perihelion (deg) */
12 cbeb0b26 2006-04-01 devnull 304.88003, /* [6] mean longitude (deg) */
14 cbeb0b26 2006-04-01 devnull -0.00125196, /* [1+6] (au/julian century) */
15 cbeb0b26 2006-04-01 devnull 0.0000251, /* [2+6] (e/julian century) */
16 cbeb0b26 2006-04-01 devnull -3.64, /* [3+6] (arcsec/julian century) */
17 cbeb0b26 2006-04-01 devnull -151.25, /* [4+6] (arcsec/julian century) */
18 cbeb0b26 2006-04-01 devnull -844.43, /* [5+6] (arcsec/julian century) */
19 cbeb0b26 2006-04-01 devnull 786449.21, /* [6+6] (arcsec/julian century) */
23 cd5bae78 2004-04-21 devnull nept(void)
25 cd5bae78 2004-04-21 devnull double pturbl, pturbb, pturbr;
26 cd5bae78 2004-04-21 devnull double lograd;
27 cd5bae78 2004-04-21 devnull double dele, enom, vnom, nd, sl;
29 cd5bae78 2004-04-21 devnull double capj, capn, eye, comg, omg;
30 cd5bae78 2004-04-21 devnull double sb, su, cu, u, b, up;
31 cd5bae78 2004-04-21 devnull double sd, ca, sa;
33 cd5bae78 2004-04-21 devnull double cy;
35 cbeb0b26 2006-04-01 devnull cy = (eday - elem[0]) / 36525.; /* per julian century */
37 cd5bae78 2004-04-21 devnull mrad = elem[1] + elem[1+6]*cy;
38 cd5bae78 2004-04-21 devnull ecc = elem[2] + elem[2+6]*cy;
40 cbeb0b26 2006-04-01 devnull cy = cy / 3600; /* arcsec/deg per julian century */
41 cd5bae78 2004-04-21 devnull incl = elem[3] + elem[3+6]*cy;
42 cd5bae78 2004-04-21 devnull node = elem[4] + elem[4+6]*cy;
43 cd5bae78 2004-04-21 devnull argp = elem[5] + elem[5+6]*cy;
45 cd5bae78 2004-04-21 devnull anom = elem[6] + elem[6+6]*cy - argp;
46 cd5bae78 2004-04-21 devnull motion = elem[6+6] / 36525. / 3600;
48 cd5bae78 2004-04-21 devnull incl *= radian;
49 cd5bae78 2004-04-21 devnull node *= radian;
50 cd5bae78 2004-04-21 devnull argp *= radian;
51 cd5bae78 2004-04-21 devnull anom = fmod(anom,360.)*radian;
53 cd5bae78 2004-04-21 devnull enom = anom + ecc*sin(anom);
55 cd5bae78 2004-04-21 devnull dele = (anom - enom + ecc * sin(enom)) /
56 cd5bae78 2004-04-21 devnull (1. - ecc*cos(enom));
57 cd5bae78 2004-04-21 devnull enom += dele;
58 cd5bae78 2004-04-21 devnull } while(fabs(dele) > converge);
59 cd5bae78 2004-04-21 devnull vnom = 2.*atan2(sqrt((1.+ecc)/(1.-ecc))*sin(enom/2.),
60 cd5bae78 2004-04-21 devnull cos(enom/2.));
61 cd5bae78 2004-04-21 devnull rad = mrad*(1. - ecc*cos(enom));
63 cd5bae78 2004-04-21 devnull lambda = vnom + argp;
64 cd5bae78 2004-04-21 devnull pturbl = 0.;
65 cd5bae78 2004-04-21 devnull lambda += pturbl*radsec;
66 cd5bae78 2004-04-21 devnull pturbb = 0.;
67 cd5bae78 2004-04-21 devnull pturbr = 0.;
70 cd5bae78 2004-04-21 devnull * reduce to the ecliptic
73 cd5bae78 2004-04-21 devnull nd = lambda - node;
74 cd5bae78 2004-04-21 devnull lambda = node + atan2(sin(nd)*cos(incl),cos(nd));
76 cd5bae78 2004-04-21 devnull sl = sin(incl)*sin(nd) + pturbb*radsec;
77 cd5bae78 2004-04-21 devnull beta = atan2(sl, pyth(sl));
79 cd5bae78 2004-04-21 devnull lograd = pturbr*2.30258509;
80 cd5bae78 2004-04-21 devnull rad *= 1. + lograd;
83 cd5bae78 2004-04-21 devnull lambda -= 1185.*radsec;
84 cd5bae78 2004-04-21 devnull beta -= 51.*radsec;
86 cd5bae78 2004-04-21 devnull motion *= radian*mrad*mrad/(rad*rad);
87 cd5bae78 2004-04-21 devnull semi = 83.33;
90 cd5bae78 2004-04-21 devnull * here begins the computation of magnitude
91 cd5bae78 2004-04-21 devnull * first find the geocentric equatorial coordinates of Saturn
94 cd5bae78 2004-04-21 devnull sd = rad*(cos(beta)*sin(lambda)*sin(obliq) +
95 cd5bae78 2004-04-21 devnull sin(beta)*cos(obliq));
96 cd5bae78 2004-04-21 devnull sa = rad*(cos(beta)*sin(lambda)*cos(obliq) -
97 cd5bae78 2004-04-21 devnull sin(beta)*sin(obliq));
98 cd5bae78 2004-04-21 devnull ca = rad*cos(beta)*cos(lambda);
99 cd5bae78 2004-04-21 devnull sd += zms;
100 cd5bae78 2004-04-21 devnull sa += yms;
101 cd5bae78 2004-04-21 devnull ca += xms;
102 cd5bae78 2004-04-21 devnull alpha = atan2(sa,ca);
103 cd5bae78 2004-04-21 devnull delta = atan2(sd,sqrt(sa*sa+ca*ca));
106 cd5bae78 2004-04-21 devnull * here are the necessary elements of Saturn's rings
107 cd5bae78 2004-04-21 devnull * cf. Exp. Supp. p. 363ff.
110 cd5bae78 2004-04-21 devnull capj = 6.9056 - 0.4322*capt;
111 cd5bae78 2004-04-21 devnull capn = 126.3615 + 3.9894*capt + 0.2403*capt2;
112 cd5bae78 2004-04-21 devnull eye = 28.0743 - 0.0128*capt;
113 cd5bae78 2004-04-21 devnull comg = 168.1179 + 1.3936*capt;
114 cd5bae78 2004-04-21 devnull omg = 42.9236 - 2.7390*capt - 0.2344*capt2;
116 cd5bae78 2004-04-21 devnull capj *= radian;
117 cd5bae78 2004-04-21 devnull capn *= radian;
118 cd5bae78 2004-04-21 devnull eye *= radian;
119 cd5bae78 2004-04-21 devnull comg *= radian;
120 cd5bae78 2004-04-21 devnull omg *= radian;
123 cd5bae78 2004-04-21 devnull * now find saturnicentric ring-plane coords of the earth
126 cd5bae78 2004-04-21 devnull sb = sin(capj)*cos(delta)*sin(alpha-capn) -
127 cd5bae78 2004-04-21 devnull cos(capj)*sin(delta);
128 cd5bae78 2004-04-21 devnull su = cos(capj)*cos(delta)*sin(alpha-capn) +
129 cd5bae78 2004-04-21 devnull sin(capj)*sin(delta);
130 cd5bae78 2004-04-21 devnull cu = cos(delta)*cos(alpha-capn);
131 cd5bae78 2004-04-21 devnull u = atan2(su,cu);
132 cd5bae78 2004-04-21 devnull b = atan2(sb,sqrt(su*su+cu*cu));
135 cd5bae78 2004-04-21 devnull * and then the saturnicentric ring-plane coords of the sun
138 cd5bae78 2004-04-21 devnull su = sin(eye)*sin(beta) +
139 cd5bae78 2004-04-21 devnull cos(eye)*cos(beta)*sin(lambda-comg);
140 cd5bae78 2004-04-21 devnull cu = cos(beta)*cos(lambda-comg);
141 cd5bae78 2004-04-21 devnull up = atan2(su,cu);
144 cd5bae78 2004-04-21 devnull * at last, the magnitude
148 cd5bae78 2004-04-21 devnull sb = sin(b);
149 cd5bae78 2004-04-21 devnull mag = -8.68 +2.52*fabs(up+omg-u)-
150 cd5bae78 2004-04-21 devnull 2.60*fabs(sb) + 1.25*(sb*sb);
152 cd5bae78 2004-04-21 devnull helio();
