sat.c (2562B)
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 }