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