uran.c (3336B)
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 }