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 23 incl *= radian; 24 node *= radian; 25 argp *= radian; 26 anom = fmod(anom,360.)*radian; 27 28 enom = anom + ecc*sin(anom); 29 do { 30 dele = (anom - enom + ecc * sin(enom)) / 31 (1. - ecc*cos(enom)); 32 enom += dele; 33 } while(fabs(dele) > 1.e-8); 34 vnom = 2.*atan2(sqrt((1.+ecc)/(1.-ecc))*sin(enom/2.), 35 cos(enom/2.)); 36 rad = mrad*(1. - ecc*cos(enom)); 37 38 lambda = vnom + argp; 39 pturbl = 0.; 40 lambda += pturbl*radsec; 41 pturbb = 0.; 42 pturbr = 0.; 43 44 /* 45 * reduce to the ecliptic 46 */ 47 48 nd = lambda - node; 49 lambda = node + atan2(sin(nd)*cos(incl),cos(nd)); 50 51 sl = sin(incl)*sin(nd) + pturbb*radsec; 52 beta = atan2(sl, pyth(sl)); 53 54 lograd = pturbr*2.30258509; 55 rad *= 1. + lograd; 56 57 58 lambda -= 1185.*radsec; 59 beta -= 51.*radsec; 60 61 motion *= radian*mrad*mrad/(rad*rad); 62 semi = 83.33; 63 64 /* 65 * here begins the computation of magnitude 66 * first find the geocentric equatorial coordinates of Saturn 67 */ 68 69 sd = rad*(cos(beta)*sin(lambda)*sin(obliq) + 70 sin(beta)*cos(obliq)); 71 sa = rad*(cos(beta)*sin(lambda)*cos(obliq) - 72 sin(beta)*sin(obliq)); 73 ca = rad*cos(beta)*cos(lambda); 74 sd += zms; 75 sa += yms; 76 ca += xms; 77 alpha = atan2(sa,ca); 78 delta = atan2(sd,sqrt(sa*sa+ca*ca)); 79 80 /* 81 * here are the necessary elements of Saturn's rings 82 * cf. Exp. Supp. p. 363ff. 83 */ 84 85 capj = 6.9056 - 0.4322*capt; 86 capn = 126.3615 + 3.9894*capt + 0.2403*capt2; 87 eye = 28.0743 - 0.0128*capt; 88 comg = 168.1179 + 1.3936*capt; 89 omg = 42.9236 - 2.7390*capt - 0.2344*capt2; 90 91 capj *= radian; 92 capn *= radian; 93 eye *= radian; 94 comg *= radian; 95 omg *= radian; 96 97 /* 98 * now find saturnicentric ring-plane coords of the earth 99 */ 100 101 sb = sin(capj)*cos(delta)*sin(alpha-capn) - 102 cos(capj)*sin(delta); 103 su = cos(capj)*cos(delta)*sin(alpha-capn) + 104 sin(capj)*sin(delta); 105 cu = cos(delta)*cos(alpha-capn); 106 u = atan2(su,cu); 107 b = atan2(sb,sqrt(su*su+cu*cu)); 108 109 /* 110 * and then the saturnicentric ring-plane coords of the sun 111 */ 112 113 su = sin(eye)*sin(beta) + 114 cos(eye)*cos(beta)*sin(lambda-comg); 115 cu = cos(beta)*cos(lambda-comg); 116 up = atan2(su,cu); 117 118 /* 119 * at last, the magnitude 120 */ 121 122 123 sb = sin(b); 124 mag = -8.68 +2.52*fabs(up+omg-u)- 125 2.60*fabs(sb) + 1.25*(sb*sb); 126 127 helio(); 128 geo(); 129 } 130