Subversion Repositories planix.SVN

#include "astro.h"

#define MAXE    (.999)  /* cant do hyperbolas */

void
comet(void)
{
        double pturbl, pturbb, pturbr;
        double lograd;
        double dele, enom, vnom, nd, sl;

        struct  elem
        {
                double  t;      /* time of perihelion */
                double  q;      /* perihelion distance */
                double  e;      /* eccentricity */
                double  i;      /* inclination */
                double  w;      /* argument of perihelion */
                double  o;      /* longitude of ascending node */
        } elem;

/*      elem = (struct elem)
        {
                etdate(1990, 5, 19.293),
                0.9362731,
                0.6940149,
                11.41096,
                198.77059,
                69.27130,
        };      /* p/schwassmann-wachmann 3, 1989d */
/*      elem = (struct elem)
        {
                etdate(1990, 4, 9.9761),
                0.349957,
                1.00038,
                58.9596,
                61.5546,
                75.2132,
        };      /* austin 3, 1989c */
/*      elem = (struct elem)
        {
                etdate(1990, 10, 24.36),
                0.9385,
                1.00038,
                131.62,
                242.58,
                138.57,
        };      /* levy 6 , 1990c */
/*      elem=(struct elem)
        {
                etdate(1996, 5, 1.3965),
                0.230035,
                0.999662,
                124.9098,
                130.2102,
                188.0430,
        };      /* C/1996 B2 (Hyakutake) */
/*      elem=(struct elem)
        {
                etdate(1997, 4, 1.13413),
                0.9141047,
                0.9950989,
                89.42932,
                130.59066,
                282.47069,
        };      /*C/1995 O1 (Hale-Bopp) */
/*      elem=(struct elem)
        {
                etdate(2000, 7, 26.1754),
                0.765126,
                0.999356,
                149.3904,
                151.0510,
                83.1909,
        };      /*C/1999 S4 (Linear) */
        elem=(struct elem)
        {
                etdate(2002, 3, 18.9784),
                0.5070601,
                0.990111,
                28.12106,
                34.6666,
                93.1206,
        };      /*C/2002 C1 (Ikeya-Zhang) */

        ecc = elem.e;
        if(ecc > MAXE)
                ecc = MAXE;
        incl = elem.i * radian;
        node = (elem.o + 0.4593) * radian;
        argp = (elem.w + elem.o + 0.4066) * radian;
        mrad = elem.q / (1-ecc);
        motion = .01720209895 * sqrt(1/(mrad*mrad*mrad))/radian;
        anom = (eday - (elem.t - 2415020)) * motion * radian;
        enom = anom + ecc*sin(anom);

        do {
                dele = (anom - enom + ecc * sin(enom)) /
                        (1 - ecc*cos(enom));
                enom += dele;
        } while(fabs(dele) > converge);

        vnom = 2*atan2(
                sqrt((1+ecc)/(1-ecc))*sin(enom/2),
                cos(enom/2));
        rad = mrad*(1-ecc*cos(enom));
        lambda = vnom + argp;
        pturbl = 0;
        lambda += pturbl*radsec;
        pturbb = 0;
        pturbr = 0;

/*
 *      reduce to the ecliptic
 */
        nd = lambda - node;
        lambda = node + atan2(sin(nd)*cos(incl),cos(nd));

        sl = sin(incl)*sin(nd) + pturbb*radsec;
        beta = atan2(sl, sqrt(1-sl*sl));

        lograd = pturbr*2.30258509;
        rad *= 1 + lograd;

        motion *= radian*mrad*mrad/(rad*rad);
        semi = 0;

        mag = 5.47 + 6.1/2.303*log(rad);

        helio();
        geo();
}