13e12c5d1SDavid du Colombier #include "astro.h"
23e12c5d1SDavid du Colombier
33e12c5d1SDavid du Colombier void
merc(void)43e12c5d1SDavid du Colombier merc(void)
53e12c5d1SDavid du Colombier {
63e12c5d1SDavid du Colombier double pturbl, pturbr;
73e12c5d1SDavid du Colombier double lograd;
83e12c5d1SDavid du Colombier double dele, enom, vnom, nd, sl;
93e12c5d1SDavid du Colombier double q0, v0, t0, j0 , s0;
103e12c5d1SDavid du Colombier double lsun, elong, ci, dlong;
113e12c5d1SDavid du Colombier
123e12c5d1SDavid du Colombier ecc = .20561421 + .00002046*capt - 0.03e-6*capt2;
133e12c5d1SDavid du Colombier incl = 7.0028806 + .0018608*capt - 18.3e-6*capt2;
143e12c5d1SDavid du Colombier node = 47.145944 + 1.185208*capt + .0001739*capt2;
153e12c5d1SDavid du Colombier argp = 75.899697 + 1.555490*capt + .0002947*capt2;
163e12c5d1SDavid du Colombier mrad = .3870986;
173e12c5d1SDavid du Colombier anom = 102.279381 + 4.0923344364*eday + 6.7e-6*capt2;
183e12c5d1SDavid du Colombier motion = 4.0923770233;
193e12c5d1SDavid du Colombier
203e12c5d1SDavid du Colombier q0 = 102.28 + 4.092334429*eday;
213e12c5d1SDavid du Colombier v0 = 212.536 + 1.602126105*eday;
223e12c5d1SDavid du Colombier t0 = -1.45 + .985604737*eday;
233e12c5d1SDavid du Colombier j0 = 225.36 + .083086735*eday;
243e12c5d1SDavid du Colombier s0 = 175.68 + .033455441*eday;
253e12c5d1SDavid du Colombier
263e12c5d1SDavid du Colombier q0 *= radian;
273e12c5d1SDavid du Colombier v0 *= radian;
283e12c5d1SDavid du Colombier t0 *= radian;
293e12c5d1SDavid du Colombier j0 *= radian;
303e12c5d1SDavid du Colombier s0 *= radian;
313e12c5d1SDavid du Colombier
323e12c5d1SDavid du Colombier incl *= radian;
333e12c5d1SDavid du Colombier node *= radian;
343e12c5d1SDavid du Colombier argp *= radian;
353e12c5d1SDavid du Colombier anom = fmod(anom, 360.)*radian;
363e12c5d1SDavid du Colombier
373e12c5d1SDavid du Colombier
383e12c5d1SDavid du Colombier enom = anom + ecc*sin(anom);
393e12c5d1SDavid du Colombier do {
403e12c5d1SDavid du Colombier dele = (anom - enom + ecc * sin(enom)) /
413e12c5d1SDavid du Colombier (1. - ecc*cos(enom));
423e12c5d1SDavid du Colombier enom += dele;
43*7dd7cddfSDavid du Colombier } while(fabs(dele) > converge);
443e12c5d1SDavid du Colombier vnom = 2.*atan2(sqrt((1.+ecc)/(1.-ecc))*sin(enom/2.),
453e12c5d1SDavid du Colombier cos(enom/2.));
463e12c5d1SDavid du Colombier rad = mrad*(1. - ecc*cos(enom));
473e12c5d1SDavid du Colombier
483e12c5d1SDavid du Colombier icosadd(mercfp, merccp);
493e12c5d1SDavid du Colombier pturbl = cosadd(2, q0, -v0);
503e12c5d1SDavid du Colombier pturbl += cosadd(2, q0, -t0);
513e12c5d1SDavid du Colombier pturbl += cosadd(2, q0, -j0);
523e12c5d1SDavid du Colombier pturbl += cosadd(2, q0, -s0);
533e12c5d1SDavid du Colombier
543e12c5d1SDavid du Colombier pturbr = cosadd(2, q0, -v0);
553e12c5d1SDavid du Colombier pturbr += cosadd(2, q0, -t0);
563e12c5d1SDavid du Colombier pturbr += cosadd(2, q0, -j0);
573e12c5d1SDavid du Colombier
583e12c5d1SDavid du Colombier /*
593e12c5d1SDavid du Colombier * reduce to the ecliptic
603e12c5d1SDavid du Colombier */
613e12c5d1SDavid du Colombier
623e12c5d1SDavid du Colombier lambda = vnom + argp + pturbl*radsec;
633e12c5d1SDavid du Colombier nd = lambda - node;
643e12c5d1SDavid du Colombier lambda = node + atan2(sin(nd)*cos(incl), cos(nd));
653e12c5d1SDavid du Colombier
663e12c5d1SDavid du Colombier sl = sin(incl)*sin(nd);
673e12c5d1SDavid du Colombier beta = atan2(sl, pyth(sl));
683e12c5d1SDavid du Colombier
693e12c5d1SDavid du Colombier lograd = pturbr*2.30258509;
703e12c5d1SDavid du Colombier rad *= 1. + lograd;
713e12c5d1SDavid du Colombier
723e12c5d1SDavid du Colombier motion *= radian*mrad*mrad/(rad*rad);
733e12c5d1SDavid du Colombier semi = 3.34;
743e12c5d1SDavid du Colombier
753e12c5d1SDavid du Colombier lsun = 99.696678 + 0.9856473354*eday;
763e12c5d1SDavid du Colombier lsun *= radian;
773e12c5d1SDavid du Colombier elong = lambda - lsun;
783e12c5d1SDavid du Colombier ci = (rad - cos(elong))/sqrt(1. + rad*rad - 2.*rad*cos(elong));
793e12c5d1SDavid du Colombier dlong = atan2(pyth(ci), ci)/radian;
803e12c5d1SDavid du Colombier mag = -.003 + .01815*dlong + .0001023*dlong*dlong;
813e12c5d1SDavid du Colombier
823e12c5d1SDavid du Colombier helio();
833e12c5d1SDavid du Colombier geo();
843e12c5d1SDavid du Colombier }
85