/* compute the position of the m`planets and return the epli value */
#include <stdio.h>
#include <math.h>

/* crude magnitudes of planets (x100) for output filtering by brightness */

#define MAGSOL -9.9
#define MAGMER -1.50
#define MAGVEN -4.0
#define MAGMAR -1.00
#define MAGJUP -2.00
#define MAGSAT 1.00
#define MAGURA 5.90
#define MAGNEP 8.00

extern double	rad;
double	htod(), atof(), kepler(), truean();
double	longi(), lati(), poly(), aint(), range();
int	cls(), helio_trans(), geo_trans(), speak();

double
planet_pos(jd, T0)
double	jd, T0;
{
	double	l, a, e, i, w1, w2, M, M0, M1, M2, M4, M5, M6, M7, M8;
	double	a0, a1, a2, a3, RA, DEC;
	double	ECC, nu, r, u, ll, b, lonpert, radpert;
	double	esun, Lsun, Cen, Stheta, Snu, Sr;
	double	N, D, epli, thapp, omeg;
	double	nu2, P, Q, S, V, W, ze, l1pert, epert, w1pert, apert;
	double	psi, H, G, eta, th;


	/* Calculate orbital elements for Mercury */
	a0 = 178.179078;
	a1 = 149474.07078;
	a2 = 0.0003011;
	a3 = 0.0;
	l = poly(a0, a1, a2, a3, T0);
	l = range(l);
	a = 0.3870986;
	a0 = 0.20561421;
	a1 = 0.00002046;
	a2 = -0.000000030;
	e = poly(a0, a1, a2, a3, T0);
	a0 = 7.002881;
	a1 = 0.0018608;
	a2 = -0.0000183;
	i = poly(a0, a1, a2, a3, T0);
	a0 = 28.753753;
	a1 = 0.3702806;
	a2 = 0.0001208;
	w1 = poly(a0, a1, a2, a3, T0);
	w1 = range(w1);
	a0 = 47.145944;
	a1 = 1.1852083;
	a2 = 0.0001739;
	w2 = poly(a0, a1, a2, a3, T0);
	w2 = range(w2);
	M1 = 102.27938 + 149472.51529 * T0 + 0.000007 * T0 * T0;
	M1 = range(M1);

	/* solving Kepler find the eccentric anomaly ECC    */

	ECC = kepler(e, M1);
	nu = truean(e, ECC);
	r = a * (1.0 - e * cos(ECC * rad));
	u = l + nu - M1 - w2;
	ll = longi(w2, i, u);
	b = lati(u, i);

	/* Now make corrections due to perturbations */
	M2 = 212.60322 + 58517.80387 * T0 + 0.001286 * T0 * T0;
	M2 = range(M2);
	M4 = 319.51913 + 19139.85475 * T0 + 0.000181 * T0 * T0;
	M4 = range(M4);
	M5 = 225.32833 + 3034.69202 * T0 - 0.000722 * T0 * T0;
	M5 = range(M5);
	M6 = 175.46622 + 1221.55147 * T0 - 0.000502 * T0 * T0;
	M6 = range(M6);
	lonpert = 0.00204 * cos((5 * M2 - 2 * M1 + 12.220) * rad)
	 + 0.00103 * cos((2 * M2 - M1 - 160.692) * rad)
	 + 0.00091 * cos((2 * M5 - M1 - 37.003) * rad)
	 + 0.00078 * cos((5 * M2 - 3 * M1 + 10.137) * rad);

	radpert = 0.000007525 * cos((2 * M5 - M1 + 53.013) * rad)
	 + 0.000006802 * cos((5 * M2 - 3 * M1 - 259.918) * rad)
	 + 0.000005457 * cos((2 * M2 - 2 * M1 - 71.188) * rad)
	 + 0.000003569 * cos((5 * M2 - M1 - 77.75) * rad);

	r = r + radpert;
	ll = ll + lonpert;

	/* Now find the Longi and radius vector of the sun */

	M = 358.47583 + 35999.04975 * T0 - 0.000150 * T0 * T0
	    -0.0000033 * T0 * T0 * T0;
	M = range(M);
	esun = 0.01675104 - 0.0000418 * T0 - 0.000000126 * T0 * T0;
	Lsun = 279.69668 + 36000.76892 * T0 + 0.0003025 * T0 * T0;
	Lsun = range(Lsun);
	Cen = (1.919460 - 0.004789 * T0 - 0.000014 * T0 * T0) * sin(M * rad)
	 + (0.020094 - 0.000100 * T0) * sin(2 * M * rad)
	 + 0.000293 * sin(3 * M * rad);

	Stheta = Lsun + Cen;
	Stheta = range(Stheta);
	Snu = M + Cen;
	Sr = 1.0000002 * (1.0 - esun * esun) / (1.0 + esun * cos(Snu * rad));

	omeg = 259.18 - 1934.142 * T0;
	thapp = Stheta - 0.00569 - 0.00479 * sin(omeg * rad);
	epli = 23.452294 - 0.0130125 * T0
	    -0.00000164 * T0 * T0 + 0.000000503 * T0 * T0 * T0;

	N = cos(epli * rad) * sin(thapp * rad);
	D = cos(thapp * rad);
	RA = atan2(N, D) / rad;
	DEC = asin(sin(epli * rad) * sin(thapp * rad)) / rad;
	speak(RA, DEC, Sr, MAGSOL, "PS", "Sol");
	/* tansformation of coordinates on Mercury and output */
	helio_trans(r, b, ll, Stheta, Sr, epli, MAGMER, "PM", "Mercury");

	/* Now start on Venus */
	a0 = 342.767053;
	a1 = 58519.21191;
	a2 = 0.0003097;
	a3 = 0.0;
	l = poly(a0, a1, a2, a3, T0);
	l = range(l);
	a = 0.7233316;
	a0 = 0.00682069;
	a1 = -0.00004774;
	a2 = 0.000000091;
	e = poly(a0, a1, a2, a3, T0);
	a0 = 3.393631;
	a1 = 0.0010058;
	a2 = -0.0000010;
	i = poly(a0, a1, a2, a3, T0);
	a0 = 54.384186;
	a1 = 0.5081861;
	a2 = -0.0013864;
	w1 = poly(a0, a1, a2, a3, T0);
	w1 = range(w1);
	a0 = 75.779647;
	a1 = 0.8998500;
	a2 = 0.0004100;
	w2 = poly(a0, a1, a2, a3, T0);
	w2 = range(w2);
	/* Venus has a long period pert that needs to be added before Kelper */
	lonpert = 0.00077 * sin((237.24 + 150.27 * T0) * rad);
	l = l + lonpert;
	M0 = M2 + lonpert;
	ECC = kepler(e, M0);
	nu = truean(e, ECC);
	r = a * (1.0 - e * cos(ECC * rad));
	u = l + nu - M0 - w2;
	ll = longi(w2, i, u);
	b = lati(u, i);

	/* now Venus perturbations */

	lonpert = 0.00313 * cos((2 * M - 2 * M2 - 148.225) * rad)
	 + 0.00198 * cos((3 * M - 3 * M2 + 2.565) * rad)
	 + 0.00136 * cos((M - M2 - 119.107) * rad)
	 + 0.00096 * cos((3 * M - 2 * M2 - 135.912) * rad)
	 + 0.00082 * cos((M5 - M2 - 208.087) * rad);
	ll = ll + lonpert;

	radpert = 0.000022501 * cos((2 * M - 2 * M2 - 58.208) * rad)
	 + 0.000019045 * cos((3 * M - 3 * M2 + 92.577) * rad)
	 + 0.000006887 * cos((M5 - M2 - 118.090) * rad)
	 + 0.000005172 * cos((M - M2 - 29.110) * rad)
	 + 0.000003620 * cos((5 * M - 4 * M2 - 104.208) * rad)
	 + 0.000003283 * cos((4 * M - 4 * M2 + 63.513) * rad)
	 + 0.000003074 * cos((2 * M5 - 2 * M2 - 55.167) * rad);
	r = r + radpert;
	helio_trans(r, b, ll, Stheta, Sr, epli, MAGVEN, "PV", "Venus");

	/* Now  start the planet Mars */
	a0 = 293.737334;
	a1 = 19141.69551;
	a2 = 0.0003107;
	a3 = 0.0;
	l = poly(a0, a1, a2, a3, T0);
	l = range(l);
	a = 1.5236883;
	a0 = 0.09331290;
	a1 = 0.000092064;
	a2 = -0.000000077;
	e = poly(a0, a1, a2, a3, T0);
	a0 = 1.850333;
	a1 = -0.0006750;
	a2 = 0.0000126;
	i = poly(a0, a1, a2, a3, T0);
	a0 = 285.431761;
	a1 = 1.0697667;
	a2 = 0.0001313;
	a3 = 0.00000414;
	w1 = poly(a0, a1, a2, a3, T0);
	w1 = range(w1);
	a0 = 48.786442;
	a1 = 0.7709917;
	a2 = -0.0000014;
	a3 = -0.00000533;
	w2 = poly(a0, a1, a2, a3, T0);
	w2 = range(w2);

	/* Mars has a long period perturbation */
	lonpert = -0.01133 * sin((3 * M5 - 8 * M4 + 4 * M) * rad)
	-0.00933 * cos((3 * M5 - 8 * M4 + 4 * M) * rad);
	l = l + lonpert;
	M0 = M4 + lonpert;
	ECC = kepler(e, M0);
	nu = truean(e, ECC);
	r = a * (1.0 - e * cos(ECC * rad));
	u = l + nu - M0 - w2;
	ll = longi(w2, i, u);
	b = lati(u, i);

	/* Now Mars Perturbations */
	lonpert = 0.00705 * cos((M5 - M4 - 48.958) * rad)
	 + 0.00607 * cos((2 * M5 - M4 - 188.350) * rad)
	 + 0.00445 * cos((2 * M5 - 2 * M4 - 191.897) * rad)
	 + 0.00388 * cos((M - 2 * M4 + 20.495) * rad)
	 + 0.00238 * cos((M - M4 + 35.097) * rad)
	 + 0.00204 * cos((2 * M - 3 * M4 + 158.638) * rad)
	 + 0.00177 * cos((3 * M4 - M2 - 57.602) * rad)
	 + 0.00136 * cos((2 * M - 4 * M4 + 154.093) * rad)
	 + 0.00104 * cos((M5 + 17.618) * rad);
	ll = ll + lonpert;

	radpert = 0.000053227 * cos((M5 - M4 + 41.1306) * rad)
	 + 0.000050989 * cos((2 * M5 - 2 * M4 - 101.9847) * rad)
	 + 0.000038278 * cos((2 * M5 - M4 - 98.3292) * rad)
	 + 0.000015996 * cos((M - M4 - 55.555) * rad)
	 + 0.000014764 * cos((2 * M - 3 * M4 + 68.622) * rad)
	 + 0.000008966 * cos((M5 - 2 * M4 + 43.615) * rad);
	/*
     * broken into two pieces for wimpy C compilers
     */
	radpert += 0.000007914 * cos((3 * M5 - 2 * M4 - 139.737) * rad)
	 + 0.000007004 * cos((2 * M5 - 3 * M4 - 102.888) * rad)
	 + 0.000006620 * cos((M - 2 * M4 + 113.202) * rad)
	 + 0.000004930 * cos((3 * M5 - 3 * M4 - 76.243) * rad)
	 + 0.000004693 * cos((3 * M - 5 * M4 + 190.603) * rad)
	 + 0.000004571 * cos((2 * M - 4 * M4 + 244.702) * rad)
	 + 0.000004409 * cos((3 * M5 - M4 - 115.828) * rad);
	r = r + radpert;
	helio_trans(r, b, ll, Stheta, Sr, epli, MAGMAR, "Pm", "Mars");

	/* Now start Jupiter */
	a0 = 238.049257;
	a1 = 3036.301986;
	a2 = 0.0003347;
	a3 = -0.00000165;
	l = poly(a0, a1, a2, a3, T0);
	l = range(l);
	a = 5.202561;
	a0 = 0.04833475;
	a1 = 0.000164180;
	a2 = -0.0000004676;
	a3 = -0.0000000017;
	e = poly(a0, a1, a2, a3, T0);
	a0 = 1.308736;
	a1 = -0.0056961;
	a2 = 0.0000039;
	a3 = 0.0;
	i = poly(a0, a1, a2, a3, T0);
	a0 = 273.277558;
	a1 = 0.5994317;
	a2 = 0.00070405;
	a3 = 0.00000508;
	w1 = poly(a0, a1, a2, a3, T0);
	w1 = range(w1);
	a0 = 99.443414;
	a1 = 1.0105300;
	a2 = 0.00035222;
	a3 = -0.00000851;
	w2 = poly(a0, a1, a2, a3, T0);
	w2 = range(w2);
	/*
     * Now start perturbation calclations
     */

	nu2 = T0 / 5.0 + 0.1;
	P = 237.47555 + 3034.9061 * T0;
	Q = 265.91650 + 1222.1139 * T0;
	S = 243.51721 + 428.4677 * T0;
	V = 5.0 * Q - 2.0 * P;
	W = 2.0 * P - 6.0 * Q + 3.0 * S;
	ze = Q - P;
	psi = S - Q;

	P = range(P) * rad;
	Q = range(Q) * rad;
	S = range(S) * rad;
	V = range(V) * rad;
	W = range(W) * rad;
	ze = range(ze) * rad;
	psi = range(psi) * rad;

	l1pert = (0.331364 - 0.010281 * nu2 - 0.004692 * nu2 * nu2) * sin(V)
	 + (0.003228 - 0.064436 * nu2 + 0.002075 * nu2 * nu2) * cos(V)
	-(0.003083 + 0.000275 * nu2 - 0.000489 * nu2 * nu2) * sin(2 * V)
	 + 0.002472 * sin(W)
	 + 0.013619 * sin(ze)
	 + 0.018472 * sin(2 * ze)
	 + 0.006717 * sin(3 * ze)
	 + 0.002775 * sin(4 * ze)
	 + (0.007275 - 0.001253 * nu2) * sin(ze) * sin(Q)
	 + 0.006417 * sin(2 * ze) * sin(Q)
	 + 0.002439 * sin(3 * ze) * sin(Q);

	l1pert = l1pert - (0.033839 + 0.001125 * nu2) * cos(ze) * sin(Q)
	-0.003767 * cos(2 * ze) * sin(Q)
	-(0.035681 + 0.001208 * nu2) * sin(ze) * cos(Q)
	-0.004261 * sin(2 * ze) * cos(Q)
	 + 0.002178 * cos(Q)
	 + (-0.006333 + 0.001161 * nu2) * cos(ze) * cos(Q)
	-0.006675 * cos(2 * ze) * cos(Q)
	-0.002664 * cos(3 * ze) * cos(Q)
	-0.002572 * sin(ze) * sin(2 * Q)
	-0.003567 * sin(2 * ze) * sin(2 * Q)
	 + 0.002094 * cos(ze) * cos(2 * Q)
	 + 0.003342 * cos(2 * ze) * cos(2 * Q);

	epert = (.0003606 + .0000130 * nu2 - .0000043 * nu2 * nu2) * sin(V)
	 + (.0001289 - .0000580 * nu2) * cos(V)
	-.0006764 * sin(ze) * sin(Q)
	-.0001110 * sin(2 * ze) * sin(Q)
	-.0000224 * sin(3 * ze) * sin(Q)
	-.0000204 * sin(Q)
	 + (.0001284 + .0000116 * nu2) * cos(ze) * sin(Q)
	 + .0000188 * cos(2 * ze) * sin(Q)
	 + (.0001460 + .0000130 * nu2) * sin(ze) * cos(Q)
	 + .0000224 * sin(2 * ze) * cos(Q)
	-.0000817 * cos(Q);

	epert = epert + .0006074 * cos(ze) * cos(Q)
	 + .0000992 * cos(2 * ze) * cos(Q)
	 + .0000508 * cos(3 * ze) * cos(Q)
	 + .0000230 * cos(4 * ze) * cos(Q)
	 + .0000108 * cos(5 * ze) * cos(Q)
	-(.0000956 + .0000073 * nu2) * sin(ze) * sin(2 * Q)
	 + .0000448 * sin(2 * ze) * sin(2 * Q)
	 + .0000137 * sin(3 * ze) * sin(2 * Q)
	 + (-.0000997 + .0000108 * nu2) * cos(ze) * sin(2 * Q)
	 + .0000480 * cos(2 * ze) * sin(2 * Q);

	epert = epert + .0000148 * cos(3 * ze) * sin(2 * Q)
	 + (-.0000956 + .0000099 * nu2) * sin(ze) * cos(2 * Q)
	 + .0000490 * sin(2 * ze) * cos(2 * Q)
	 + .0000158 * sin(3 * ze) * cos(2 * Q)
	 + .0000179 * cos(2 * Q)
	 + (.0001024 + .0000075 * nu2) * cos(ze) * cos(2 * Q)
	-.0000437 * cos(2 * ze) * cos(2 * Q)
	-.0000132 * cos(3 * ze) * cos(2 * Q);

	w1pert = (0.007192 - 0.003147 * nu2) * sin(V)
	 + (-0.020428 - 0.000675 * nu2 + 0.000197 * nu2 * nu2) * cos(V)
	 + (0.007269 + 0.000672 * nu2) * sin(ze) * sin(Q)
	-0.004344 * sin(Q)
	 + 0.034036 * cos(ze) * sin(Q)
	 + 0.005614 * cos(2 * ze) * sin(Q)
	 + 0.002964 * cos(3 * ze) * sin(Q)
	 + 0.037761 * sin(ze) * cos(Q);

	w1pert = w1pert
	     + 0.006158 * sin(2 * ze) * cos(Q)
	-0.006603 * cos(ze) * cos(Q)
	-0.005356 * sin(ze) * sin(2 * Q)
	 + 0.002722 * sin(2 * ze) * sin(2 * Q)
	 + 0.004483 * cos(ze) * sin(2 * Q)
	-0.002642 * cos(2 * ze) * sin(2 * Q)
	 + 0.004403 * sin(ze) * cos(2 * Q)
	-0.002536 * sin(2 * ze) * cos(2 * Q)
	 + 0.005547 * cos(ze) * cos(2 * Q)
	-0.002689 * cos(2 * ze) * cos(2 * Q);

	lonpert = l1pert - (w1pert / e);
	l = l + l1pert;
	M0 = M5 + lonpert;
	e = e + epert;
	w1 = w1 + w1pert;

	apert = -.000263 * cos(V)
	 + .000205 * cos(ze)
	 + .000693 * cos(2 * ze)
	 + .000312 * cos(3 * ze)
	 + .000147 * cos(4 * ze)
	 + .000299 * sin(ze) * sin(Q)
	 + .000181 * cos(2 * ze) * sin(Q)
	 + .000204 * sin(2 * ze) * cos(Q)
	 + .000111 * sin(3 * ze) * cos(Q)
	-.000337 * cos(ze) * cos(Q)
	-.000111 * cos(2 * ze) * cos(Q);

	a = a + apert;
	ECC = kepler(e, M0);
	nu = truean(e, ECC);
	r = a * (1.0 - e * cos(ECC * rad));
	u = l + nu - M0 - w2;
	ll = longi(w2, i, u);
	b = lati(u, i);

	helio_trans(r, b, ll, Stheta, Sr, epli, MAGJUP, "PJ", "Jupiter");

	/* Now start Saturn */

	a0 = 266.564377;
	a1 = 1223.509884;
	a2 = 0.0003245;
	a3 = -0.0000058;
	l = poly(a0, a1, a2, a3, T0);
	l = range(l);
	a = 9.554747;
	a0 = 0.05589232;
	a1 = -0.00034550;
	a2 = -0.000000728;
	a3 = 0.00000000074;
	e = poly(a0, a1, a2, a3, T0);
	a0 = 2.492519;
	a1 = -0.0039189;
	a2 = -0.00001549;
	a3 = 0.00000004;
	i = poly(a0, a1, a2, a3, T0);
	a0 = 338.307800;
	a1 = 1.0852207;
	a2 = 0.00097854;
	a3 = 0.00000992;
	w1 = poly(a0, a1, a2, a3, T0);
	w1 = range(w1);
	a0 = 112.790414;
	a1 = 0.8731951;
	a2 = -0.00015218;
	a3 = -0.00000531;
	w2 = poly(a0, a1, a2, a3, T0);
	w2 = range(w2);

	/* Now Saturn's perturbations */

	l1pert = (-0.814181 + 0.018150 * nu2 + 0.016714 * nu2 * nu2) * sin(V)
	 + (-0.010497 + 0.160906 * nu2 - 0.004100 * nu2 * nu2) * cos(V)
	 + 0.007581 * sin(2 * V)
	-0.007986 * sin(W)
	-0.148811 * sin(ze)
	-0.040786 * sin(2 * ze)
	-0.015208 * sin(3 * ze)
	-0.006339 * sin(4 * ze)
	-0.006244 * sin(Q);
	l1pert = l1pert
	     + (0.008931 + 0.002728 * nu2) * sin(ze) * sin(Q)
	-0.016500 * sin(2 * ze) * sin(Q)
	-0.005775 * sin(3 * ze) * sin(Q)
	 + (0.081344 + 0.003206 * nu2) * cos(ze) * sin(Q)
	 + 0.015019 * cos(2 * ze) * sin(Q)
	 + (0.085581 + 0.002494 * nu2) * sin(ze) * cos(Q)
	 + (0.025328 - 0.003117 * nu2) * cos(ze) * cos(Q);
	l1pert = l1pert
	     + 0.014394 * cos(2 * ze) * cos(Q)
	 + 0.006319 * cos(3 * ze) * cos(Q)
	 + 0.006369 * sin(ze) * sin(2 * Q)
	 + 0.009156 * sin(2 * ze) * sin(2 * Q)
	 + 0.007525 * sin(3 * psi) * sin(2 * Q)
	-0.005236 * cos(ze) * cos(2 * Q)
	-0.007736 * cos(2 * ze) * cos(2 * Q)
	-0.007528 * cos(3 * psi) * cos(2 * Q);

	epert = (-.0007927 + .0002548 * nu2 + .0000091 * nu2 * nu2) * sin(V)
	 + (.0013381 + .0001226 * nu2 - .0000253 * nu2 * nu2) * cos(V)
	 + (.0000248 - .0000121 * nu2) * sin(2 * V)
	-(.0000305 + .0000091 * nu2) * cos(2 * V)
	 + .0000412 * sin(2 * ze)
	 + .0012415 * sin(Q)
	 + (.0000390 - .0000617 * nu2) * sin(ze) * sin(Q)
	 + (.0000165 - .0000204 * nu2) * sin(2 * ze) * sin(Q)
	 + .0026599 * cos(ze) * sin(Q)
	-.0004687 * cos(2 * ze) * sin(Q);
	epert = epert
	    -.0001870 * cos(3 * ze) * sin(Q)
	-.0000821 * cos(4 * ze) * sin(Q)
	-.0000377 * cos(5 * ze) * sin(Q)
	 + .0000497 * cos(2 * psi) * sin(Q)
	 + (.0000163 - .0000611 * nu2) * cos(Q)
	-.0012696 * sin(ze) * cos(Q)
	-.0004200 * sin(2 * ze) * cos(Q)
	-.0001503 * sin(3 * ze) * cos(Q)
	-.0000619 * sin(4 * ze) * cos(Q)
	-.0000268 * sin(5 * ze) * cos(Q);
	epert = epert
	    -(.0000282 + .0001306 * nu2) * cos(ze) * cos(Q)
	 + (-.0000086 + .0000230 * nu2) * cos(2 * ze) * cos(Q)
	 + .0000461 * sin(2 * psi) * cos(Q)
	-.0000350 * sin(2 * Q)
	 + (.0002211 - .0000286 * nu2) * sin(ze) * sin(2 * Q)
	-.0002208 * sin(2 * ze) * sin(2 * Q)
	-.0000568 * sin(3 * ze) * sin(2 * Q)
	-.0000346 * sin(4 * ze) * sin(2 * Q)
	-(.0002780 + .0000222 * nu2) * cos(ze) * sin(2 * Q)
	 + (.0002022 + .0000263 * nu2) * cos(2 * ze) * sin(2 * Q);
	epert = epert
	     + .0000248 * cos(3 * ze) * sin(2 * Q)
	 + .0000242 * sin(3 * psi) * sin(2 * Q)
	 + .0000467 * cos(3 * psi) * sin(2 * Q)
	-.0000490 * cos(2 * Q)
	-(.0002842 + .0000279 * nu2) * sin(ze) * cos(2 * Q)
	 + (.0000128 + .0000226 * nu2) * sin(2 * ze) * cos(2 * Q)
	 + .0000224 * sin(3 * ze) * cos(2 * Q)
	 + (-.0001594 + .0000282 * nu2) * cos(ze) * cos(2 * Q)
	 + (.0002162 - .0000207 * nu2) * cos(2 * ze) * cos(2 * Q)
	 + .0000561 * cos(3 * ze) * cos(2 * Q);
	epert = epert
	     + .0000343 * cos(4 * ze) * cos(2 * Q)
	 + .0000469 * sin(3 * psi) * cos(2 * Q)
	-.0000242 * cos(3 * psi) * cos(2 * Q)
	-.0000205 * sin(ze) * sin(3 * Q)
	 + .0000262 * sin(3 * ze) * sin(3 * Q)
	 + .0000208 * cos(ze) * cos(3 * Q)
	-.0000271 * cos(3 * ze) * cos(3 * Q)
	-.0000382 * cos(3 * ze) * sin(4 * Q)
	-.0000376 * sin(3 * ze) * cos(4 * Q);

	w1pert = (0.077108 + 0.007186 * nu2 - 0.001533 * nu2 * nu2) * sin(V)
	 + (0.045803 - 0.014766 * nu2 - 0.000536 * nu2 * nu2) * cos(V)
	-0.007075 * sin(ze)
	-0.075825 * sin(ze) * sin(Q)
	-0.024839 * sin(2 * ze) * sin(Q)
	-0.008631 * sin(3 * ze) * sin(Q)
	-0.072586 * cos(Q)
	-0.150383 * cos(ze) * cos(Q)
	 + 0.026897 * cos(2 * ze) * cos(Q)
	 + 0.010053 * cos(3 * ze) * cos(Q);
	w1pert = w1pert
	    -(0.013597 + 0.001719 * nu2) * sin(ze) * sin(2 * Q)
	 + (-0.007742 + 0.001517 * nu2) * cos(ze) * sin(2 * Q)
	 + (0.013586 - 0.001375 * nu2) * cos(2 * ze) * sin(2 * Q)
	 + (-0.013667 + 0.001239 * nu2) * sin(ze) * cos(2 * Q)
	 + 0.011981 * sin(2 * ze) * cos(2 * Q)
	 + (0.014861 + 0.001136 * nu2) * cos(ze) * cos(2 * Q)
	-(0.013064 + 0.001628 * nu2) * cos(2 * ze) * cos(2 * Q);

	lonpert = l1pert - (w1pert / e);
	l = l + l1pert;
	M0 = M6 + lonpert;
	e = e + epert;
	w1 = w1 + w1pert;

	apert = .000572 * sin(V) - .001590 * sin(2 * ze) * cos(Q)
	 + .002933 * cos(V) - .000647 * sin(3 * ze) * cos(Q)
	 + .033629 * cos(ze) - .000344 * sin(4 * ze) * cos(Q)
	-.003081 * cos(2 * ze) + .002885 * cos(ze) * cos(Q)
	-.001423 * cos(3 * ze) + (.002172 + .000102 * nu2) * cos(2 * ze) * cos(Q)
	-.000671 * cos(4 * ze) + .000296 * cos(3 * ze) * cos(Q)
	-.000320 * cos(5 * ze) - .000267 * sin(2 * ze) * sin(2 * Q);
	apert = apert
	     + .001098 * sin(Q) - .000778 * cos(ze) * sin(2 * Q)
	-.002812 * sin(ze) * sin(Q) + .000495 * cos(2 * ze) * sin(2 * Q)
	 + .000688 * sin(2 * ze) * sin(Q) + .000250 * cos(3 * ze) * sin(2 * Q);
	apert = apert
	    -.000393 * sin(3 * ze) * sin(Q)
	-.000228 * sin(4 * ze) * sin(Q)
	 + .002138 * cos(ze) * sin(Q)
	-.000999 * cos(2 * ze) * sin(Q)
	-.000642 * cos(3 * ze) * sin(Q)
	-.000325 * cos(4 * ze) * sin(Q)
	-.000890 * cos(Q)
	 + .002206 * sin(ze) * cos(Q);
	apert = apert
	    -.000856 * sin(ze) * cos(2 * Q)
	 + .000441 * sin(2 * ze) * cos(2 * Q)
	 + .000296 * cos(2 * ze) * cos(2 * Q)
	 + .000211 * cos(3 * ze) * cos(2 * Q)
	-.000427 * sin(ze) * sin(3 * Q)
	 + .000398 * sin(3 * ze) * sin(3 * Q)
	 + .000344 * cos(ze) * cos(3 * Q)
	-.000427 * cos(3 * ze) * cos(3 * Q);

	a = a + apert;
	ECC = kepler(e, M0);
	nu = truean(e, ECC);
	r = a * (1.0 - e * cos(ECC * rad));
	u = l + nu - M0 - w2;
	ll = longi(w2, i, u);
	b = lati(u, i);

	b = b
	     + 0.000747 * cos(ze) * sin(Q)
	 + 0.001069 * cos(ze) * cos(Q)
	 + 0.002108 * sin(2 * ze) * sin(2 * Q)
	 + 0.001261 * cos(2 * ze) * sin(2 * Q)
	 + 0.001236 * sin(2 * ze) * cos(2 * Q)
	-0.002075 * cos(2 * ze) * cos(2 * Q);

	helio_trans(r, b, ll, Stheta, Sr, epli, MAGSAT, "Ps", "Saturn");

	/* Now Start on Uranus */
	a0 = 244.197470;
	a1 = 429.863546;
	a2 = 0.0003160;
	a3 = -0.00000060;
	l = poly(a0, a1, a2, a3, T0);
	l = range(l);
	a = 19.21814;
	a0 = 0.0463444;
	a1 = -0.00002658;
	a2 = 0.000000077;
	a3 = 0.0;
	e = poly(a0, a1, a2, a3, T0);
	a0 = 0.772464;
	a1 = 0.0006253;
	a2 = 0.0000395;
	i = poly(a0, a1, a2, a3, T0);
	a0 = 98.071581;
	a1 = 0.9857650;
	a2 = -0.0010745;
	a3 = -0.00000061;
	w1 = poly(a0, a1, a2, a3, T0);
	w1 = range(w1);
	a0 = 73.477111;
	a1 = 0.4986678;
	a2 = 0.0013117;
	a3 = 0.0;
	w2 = poly(a0, a1, a2, a3, T0);
	w2 = range(w2);
	M7 = 72.64878 + 428.37911 * T0 + 0.000079 * T0 * T0;
	M7 = range(M7);
	/* now perturbation corrections for Uranus */
	G = 83.76922 + 218.4901 * T0;
	S = S / rad;
	P = P / rad;
	Q = Q / rad;
	H = 2.0 * G - S;
	ze = S - P;
	eta = S - Q;
	th = G - S;
	S = S * rad;
	G = range(G) * rad;
	P = P * rad;
	Q = Q * rad;
	H = range(H) * rad;
	ze = range(ze) * rad;
	eta = range(eta) * rad;
	th = range(th) * rad;

	l1pert = (0.864319 - 0.001583 * nu2) * sin(H)
	 + (0.082222 - 0.006833 * nu2) * cos(H)
	 + 0.036017 * sin(2 * H)
	-0.003019 * cos(2 * H)
	 + 0.008122 * sin(W);

	epert = (-.0003349 + .0000163 * nu2) * sin(H)
	 + .0020981 * cos(H)
	 + .0001311 * cos(H);

	w1pert = 0.120303 * sin(H)
	 + (0.019472 - 0.000947 * nu2) * cos(H)
	 + 0.006197 * sin(2 * H);

	lonpert = l1pert - (w1pert / e);
	l = l + l1pert;
	M0 = M7 + lonpert;
	e = e + epert;
	w1 = w1 + w1pert;

	a = a - 0.003825 * cos(H);
	ECC = kepler(e, M0);
	nu = truean(e, ECC);
	r = a * (1.0 - e * cos(ECC * rad));
	u = l + nu - M0 - w2;
	ll = longi(w2, i, u);
	b = lati(u, i);

	ll = ll
	     + (0.010122 - 0.000988 * nu2) * sin(S + eta)
	 + (-0.038581 + 0.002031 * nu2 - 0.001910 * nu2 * nu2) * cos(S + eta)
	 + (0.034964 - 0.001038 * nu2 + 0.000868 * nu2 * nu2) * cos(2 * S + eta)
	 + 0.005594 * sin(S + 3 * th);
	ll = ll
	    -0.014808 * sin(ze)
	-0.005794 * sin(eta)
	 + 0.002347 * cos(eta)
	 + 0.009872 * sin(th)
	 + 0.008803 * sin(2 * th)
	-0.004308 * sin(3 * th);

	b = b
	     + (0.000458 * sin(eta) - 0.000642 * cos(eta) - 0.000517 * cos(4 * th))
	*sin(S)
	-(0.000347 * sin(eta) + 0.000853 * cos(eta) + 0.000517 * sin(4 * eta))
	*cos(S)
	 + 0.000403 * (cos(2 * th) * sin(2 * S) + sin(2 * th) * cos(2 * S));

	r = r
	    -.025948
	     + .004985 * cos(ze)
	-.001230 * cos(S)
	 + .003354 * cos(eta)
	 + (.005795 * cos(S) - .001165 * sin(S) + .001388 * cos(2 * S)) * sin(eta)
	 + (.001351 * cos(S) + .005702 * sin(S) + .001388 * sin(2 * S)) * cos(eta)
	 + .000904 * cos(2 * th)
	 + .000894 * (cos(th) - cos(3 * th));
	helio_trans(r, b, ll, Stheta, Sr, epli, MAGURA, "PU", "Uranus");

	/* Now start Neptune */
	a0 = 84.457994;
	a1 = 219.885914;
	a2 = 0.0003205;
	a3 = -0.00000060;
	l = poly(a0, a1, a2, a3, T0);
	l = range(l);
	a = 30.10957;
	a0 = 0.00899704;
	a1 = 0.000006330;
	a2 = -0.000000002;
	a3 = 0.0;
	e = poly(a0, a1, a2, a3, T0);
	a0 = 1.779242;
	a1 = -0.0095436;
	a2 = -0.0000091;
	i = poly(a0, a1, a2, a3, T0);
	a0 = 276.045975;
	a1 = 0.3256394;
	a2 = 0.00014095;
	a3 = 0.000004113;
	w1 = poly(a0, a1, a2, a3, T0);
	w1 = range(w1);
	a0 = 130.681389;
	a1 = 1.0989350;
	a2 = 0.00024987;
	a3 = -0.000004718;
	w2 = poly(a0, a1, a2, a3, T0);
	w2 = range(w2);
	M8 = 37.73063 + 218.46134 * T0 - 0.000070 * T0 * T0;

	/* now perturbation corrections for neptune */
	G = G / rad;
	P = P / rad;
	Q = Q / rad;
	S = S / rad;
	ze = G - P;
	eta = G - Q;
	th = G - S;
	G = G * rad;
	P = P * rad;
	Q = Q * rad;
	S = S * rad;
	ze = range(ze) * rad;
	eta = range(eta) * rad;
	th = range(th) * rad;

	l1pert = (-0.589833 + 0.001089 * nu2) * sin(H)
	 + (-0.056094 + 0.004658 * nu2) * cos(H)
	-0.024286 * sin(2 * H);

	epert = .0004389 * sin(H)
	 + .0004262 * cos(H)
	 + .0001129 * sin(2 * H)
	 + .0001089 * cos(2 * H);

	w1pert = 0.024039 * sin(H)
	-0.025303 * cos(H)
	 + 0.006206 * sin(2 * H)
	-0.005992 * cos(2 * H);


	lonpert = l1pert - (w1pert / e);
	l = l + l1pert;
	M0 = M8 + lonpert;
	e = e + epert;
	w1 = w1 + w1pert;

	a = a - 0.000817 * sin(H)
	 + 0.008189 * cos(H)
	 + 0.000781 * cos(2 * H);

	ECC = kepler(e, M0);
	nu = truean(e, ECC);
	r = a * (1.0 - e * cos(ECC * rad));
	u = l + nu - M0 - w2;
	ll = longi(w2, i, u);
	b = lati(u, i);

	ll = ll
	    -0.009556 * sin(ze)
	-0.005178 * sin(eta)
	 + 0.002572 * sin(2 * th)
	-0.002972 * cos(2 * th) * sin(G)
	-0.002833 * sin(2 * th) * cos(G);

	b = b
	     + 0.000336 * cos(2 * th) * sin(G)
	 + 0.000364 * sin(2 * th) * cos(G);

	r = r
	    -.040596
	     + .004992 * cos(ze)
	 + .002744 * cos(eta)
	 + .002044 * cos(th)
	 + .001051 * cos(2 * th);
	helio_trans(r, b, ll, Stheta, Sr, epli, MAGNEP, "PN", "Neptune");

	return epli;
}