1 cd5bae78 2004-04-21 devnull #include "astro.h"
  2 cd5bae78 2004-04-21 devnull 
  3 cd5bae78 2004-04-21 devnull void
  4 cd5bae78 2004-04-21 devnull sat(void)
  5 cd5bae78 2004-04-21 devnull {
  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;
  9 cd5bae78 2004-04-21 devnull 
 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;
 13 cd5bae78 2004-04-21 devnull 
 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.;
 21 cd5bae78 2004-04-21 devnull 
 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;
 26 cd5bae78 2004-04-21 devnull 
 27 cd5bae78 2004-04-21 devnull 	enom = anom + ecc*sin(anom);
 28 cd5bae78 2004-04-21 devnull 	do {
 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));
 36 cd5bae78 2004-04-21 devnull 
 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.;
 42 cd5bae78 2004-04-21 devnull 
 43 cd5bae78 2004-04-21 devnull /*
 44 cd5bae78 2004-04-21 devnull  *	reduce to the ecliptic
 45 cd5bae78 2004-04-21 devnull  */
 46 cd5bae78 2004-04-21 devnull 
 47 cd5bae78 2004-04-21 devnull 	nd = lambda - node;
 48 cd5bae78 2004-04-21 devnull 	lambda = node + atan2(sin(nd)*cos(incl),cos(nd));
 49 cd5bae78 2004-04-21 devnull 
 50 cd5bae78 2004-04-21 devnull 	sl = sin(incl)*sin(nd) + pturbb*radsec;
 51 cd5bae78 2004-04-21 devnull 	beta = atan2(sl, pyth(sl));
 52 cd5bae78 2004-04-21 devnull 
 53 cd5bae78 2004-04-21 devnull 	lograd = pturbr*2.30258509;
 54 cd5bae78 2004-04-21 devnull 	rad *= 1. + lograd;
 55 cd5bae78 2004-04-21 devnull 
 56 cd5bae78 2004-04-21 devnull 
 57 cd5bae78 2004-04-21 devnull 	lambda -= 1185.*radsec;
 58 cd5bae78 2004-04-21 devnull 	beta -= 51.*radsec;
 59 cd5bae78 2004-04-21 devnull 
 60 cd5bae78 2004-04-21 devnull 	motion *= radian*mrad*mrad/(rad*rad);
 61 cd5bae78 2004-04-21 devnull 	semi = 83.33;
 62 cd5bae78 2004-04-21 devnull 
 63 cd5bae78 2004-04-21 devnull /*
 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
 66 cd5bae78 2004-04-21 devnull  */
 67 cd5bae78 2004-04-21 devnull 
 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));
 78 cd5bae78 2004-04-21 devnull 
 79 cd5bae78 2004-04-21 devnull /*
 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.
 82 cd5bae78 2004-04-21 devnull  */
 83 cd5bae78 2004-04-21 devnull 
 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;
 89 cd5bae78 2004-04-21 devnull 
 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;
 95 cd5bae78 2004-04-21 devnull 
 96 cd5bae78 2004-04-21 devnull /*
 97 cd5bae78 2004-04-21 devnull  *	now find saturnicentric ring-plane coords of the earth
 98 cd5bae78 2004-04-21 devnull  */
 99 cd5bae78 2004-04-21 devnull 
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));
107 cd5bae78 2004-04-21 devnull 
108 cd5bae78 2004-04-21 devnull /*
109 cd5bae78 2004-04-21 devnull  *	and then the saturnicentric ring-plane coords of the sun
110 cd5bae78 2004-04-21 devnull  */
111 cd5bae78 2004-04-21 devnull 
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);
116 cd5bae78 2004-04-21 devnull 
117 cd5bae78 2004-04-21 devnull /*
118 cd5bae78 2004-04-21 devnull  *	at last, the magnitude
119 cd5bae78 2004-04-21 devnull  */
120 cd5bae78 2004-04-21 devnull 
121 cd5bae78 2004-04-21 devnull 
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);
125 cd5bae78 2004-04-21 devnull 
126 cd5bae78 2004-04-21 devnull 	helio();
127 cd5bae78 2004-04-21 devnull 	geo();
128 cd5bae78 2004-04-21 devnull }