merc.c (1873B)
1 #include "astro.h" 2 3 void 4 merc(void) 5 { 6 double pturbl, pturbr; 7 double lograd; 8 double dele, enom, vnom, nd, sl; 9 double q0, v0, t0, j0 , s0; 10 double lsun, elong, ci, dlong; 11 12 ecc = .20561421 + .00002046*capt - 0.03e-6*capt2; 13 incl = 7.0028806 + .0018608*capt - 18.3e-6*capt2; 14 node = 47.145944 + 1.185208*capt + .0001739*capt2; 15 argp = 75.899697 + 1.555490*capt + .0002947*capt2; 16 mrad = .3870986; 17 anom = 102.279381 + 4.0923344364*eday + 6.7e-6*capt2; 18 motion = 4.0923770233; 19 20 q0 = 102.28 + 4.092334429*eday; 21 v0 = 212.536 + 1.602126105*eday; 22 t0 = -1.45 + .985604737*eday; 23 j0 = 225.36 + .083086735*eday; 24 s0 = 175.68 + .033455441*eday; 25 26 q0 *= radian; 27 v0 *= radian; 28 t0 *= radian; 29 j0 *= radian; 30 s0 *= radian; 31 32 incl *= radian; 33 node *= radian; 34 argp *= radian; 35 anom = fmod(anom, 360.)*radian; 36 37 38 enom = anom + ecc*sin(anom); 39 do { 40 dele = (anom - enom + ecc * sin(enom)) / 41 (1. - ecc*cos(enom)); 42 enom += dele; 43 } while(fabs(dele) > converge); 44 vnom = 2.*atan2(sqrt((1.+ecc)/(1.-ecc))*sin(enom/2.), 45 cos(enom/2.)); 46 rad = mrad*(1. - ecc*cos(enom)); 47 48 icosadd(mercfp, merccp); 49 pturbl = cosadd(2, q0, -v0); 50 pturbl += cosadd(2, q0, -t0); 51 pturbl += cosadd(2, q0, -j0); 52 pturbl += cosadd(2, q0, -s0); 53 54 pturbr = cosadd(2, q0, -v0); 55 pturbr += cosadd(2, q0, -t0); 56 pturbr += cosadd(2, q0, -j0); 57 58 /* 59 * reduce to the ecliptic 60 */ 61 62 lambda = vnom + argp + pturbl*radsec; 63 nd = lambda - node; 64 lambda = node + atan2(sin(nd)*cos(incl), cos(nd)); 65 66 sl = sin(incl)*sin(nd); 67 beta = atan2(sl, pyth(sl)); 68 69 lograd = pturbr*2.30258509; 70 rad *= 1. + lograd; 71 72 motion *= radian*mrad*mrad/(rad*rad); 73 semi = 3.34; 74 75 lsun = 99.696678 + 0.9856473354*eday; 76 lsun *= radian; 77 elong = lambda - lsun; 78 ci = (rad - cos(elong))/sqrt(1. + rad*rad - 2.*rad*cos(elong)); 79 dlong = atan2(pyth(ci), ci)/radian; 80 mag = -.003 + .01815*dlong + .0001023*dlong*dlong; 81 82 helio(); 83 geo(); 84 }