13e12c5d1SDavid du Colombier #include "astro.h"
23e12c5d1SDavid du Colombier
33e12c5d1SDavid du Colombier void
mars(void)43e12c5d1SDavid du Colombier mars(void)
53e12c5d1SDavid du Colombier {
63e12c5d1SDavid du Colombier double pturbl, pturbb, pturbr;
73e12c5d1SDavid du Colombier double lograd;
83e12c5d1SDavid du Colombier double dele, enom, vnom, nd, sl;
93e12c5d1SDavid du Colombier double lsun, elong, ci, dlong;
103e12c5d1SDavid du Colombier
113e12c5d1SDavid du Colombier
123e12c5d1SDavid du Colombier ecc = .09331290 + .000092064*capt;
133e12c5d1SDavid du Colombier incl = 1.850333 - 6.75e-4*capt;
143e12c5d1SDavid du Colombier node = 48.786442 + .770992*capt;
153e12c5d1SDavid du Colombier argp = 334.218203 + 1.840758*capt + 1.30e-4*capt2;
163e12c5d1SDavid du Colombier mrad = 1.5236915;
173e12c5d1SDavid du Colombier anom = 319.529425 + .5240207666*eday + 1.808e-4*capt2;
183e12c5d1SDavid du Colombier motion = 0.5240711638;
193e12c5d1SDavid du Colombier
203e12c5d1SDavid du Colombier
213e12c5d1SDavid du Colombier incl = incl*radian;
223e12c5d1SDavid du Colombier node = node*radian;
233e12c5d1SDavid du Colombier argp = argp*radian;
243e12c5d1SDavid du Colombier anom = fmod(anom,360.)*radian;
253e12c5d1SDavid du Colombier
263e12c5d1SDavid du Colombier enom = anom + ecc*sin(anom);
273e12c5d1SDavid du Colombier do {
283e12c5d1SDavid du Colombier dele = (anom - enom + ecc * sin(enom)) /
293e12c5d1SDavid du Colombier (1. - ecc*cos(enom));
303e12c5d1SDavid du Colombier enom += dele;
31*7dd7cddfSDavid du Colombier } while(fabs(dele) > converge);
323e12c5d1SDavid du Colombier vnom = 2.*atan2(sqrt((1.+ecc)/(1.-ecc))*sin(enom/2.),cos(enom/2.));
333e12c5d1SDavid du Colombier rad = mrad*(1. - ecc*cos(enom));
343e12c5d1SDavid du Colombier
353e12c5d1SDavid du Colombier lambda = vnom + argp;
363e12c5d1SDavid du Colombier pturbl = 0.;
373e12c5d1SDavid du Colombier lambda = lambda + pturbl*radsec;
383e12c5d1SDavid du Colombier pturbb = 0.;
393e12c5d1SDavid du Colombier pturbr = 0.;
403e12c5d1SDavid du Colombier
413e12c5d1SDavid du Colombier /*
423e12c5d1SDavid du Colombier * reduce to the ecliptic
433e12c5d1SDavid du Colombier */
443e12c5d1SDavid du Colombier
453e12c5d1SDavid du Colombier nd = lambda - node;
463e12c5d1SDavid du Colombier lambda = node + atan2(sin(nd)*cos(incl),cos(nd));
473e12c5d1SDavid du Colombier
483e12c5d1SDavid du Colombier sl = sin(incl)*sin(nd) + pturbb*radsec;
493e12c5d1SDavid du Colombier beta = atan2(sl, pyth(sl));
503e12c5d1SDavid du Colombier
513e12c5d1SDavid du Colombier lograd = pturbr*2.30258509;
523e12c5d1SDavid du Colombier rad *= 1. + lograd;
533e12c5d1SDavid du Colombier
543e12c5d1SDavid du Colombier
553e12c5d1SDavid du Colombier motion *= radian*mrad*mrad/(rad*rad);
563e12c5d1SDavid du Colombier semi = 4.68;
573e12c5d1SDavid du Colombier
583e12c5d1SDavid du Colombier lsun = 99.696678 + 0.9856473354*eday;
593e12c5d1SDavid du Colombier lsun *= radian;
603e12c5d1SDavid du Colombier elong = lambda - lsun;
613e12c5d1SDavid du Colombier ci = (rad - cos(elong))/sqrt(1. + rad*rad - 2.*rad*cos(elong));
623e12c5d1SDavid du Colombier dlong = atan2(pyth(ci), ci)/radian;
633e12c5d1SDavid du Colombier mag = -1.30 + .01486*dlong;
643e12c5d1SDavid du Colombier
653e12c5d1SDavid du Colombier helio();
663e12c5d1SDavid du Colombier geo();
673e12c5d1SDavid du Colombier }
68