1 #include "astro.h"
2 
3 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 }