Subversion Repositories planix.SVN

Rev

Rev 2 | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
2 - 1 
#include "astro.h"
2 
 
3 
#define	MAXE	(.999)	/* cant do hyperbolas */
4 
 
5 
void
6 
comet(void)
7 
{
8 
	double pturbl, pturbb, pturbr;
9 
	double lograd;
10 
	double dele, enom, vnom, nd, sl;
11 
 
12 
	struct	elem
13 
	{
14 
		double	t;	/* time of perihelion */
15 
		double	q;	/* perihelion distance */
16 
		double	e;	/* eccentricity */
17 
		double	i;	/* inclination */
18 
		double	w;	/* argument of perihelion */
19 
		double	o;	/* longitude of ascending node */
20 
	} elem;
21 
 
22 
/*	elem = (struct elem)
23 
	{
24 
		etdate(1990, 5, 19.293),
25 
		0.9362731,
26 
		0.6940149,
27 
		11.41096,
28 
		198.77059,
29 
		69.27130,
30 
	};	/* p/schwassmann-wachmann 3, 1989d */
31 
/*	elem = (struct elem)
32 
	{
33 
		etdate(1990, 4, 9.9761),
34 
		0.349957,
35 
		1.00038,
36 
		58.9596,
37 
		61.5546,
38 
		75.2132,
39 
	};	/* austin 3, 1989c */
40 
/*	elem = (struct elem)
41 
	{
42 
		etdate(1990, 10, 24.36),
43 
		0.9385,
44 
		1.00038,
45 
		131.62,
46 
		242.58,
47 
		138.57,
48 
	};	/* levy 6 , 1990c */
49 
/*	elem=(struct elem)
50 
	{
51 
		etdate(1996, 5, 1.3965),
52 
		0.230035,
53 
		0.999662,
54 
		124.9098,
55 
		130.2102,
56 
		188.0430,
57 
	};	/* C/1996 B2 (Hyakutake) */
58 
/*	elem=(struct elem)
59 
	{
60 
		etdate(1997, 4, 1.13413),
61 
		0.9141047,
62 
		0.9950989,
63 
		89.42932,
64 
		130.59066,
65 
		282.47069,
66 
	};	/*C/1995 O1 (Hale-Bopp) */
67 
/*	elem=(struct elem)
68 
	{
69 
		etdate(2000, 7, 26.1754),
70 
		0.765126,
71 
		0.999356,
72 
		149.3904,
73 
		151.0510,
74 
		83.1909,
75 
	};	/*C/1999 S4 (Linear) */
76 
	elem=(struct elem)
77 
	{
78 
		etdate(2002, 3, 18.9784),
79 
		0.5070601,
80 
		0.990111,
81 
		28.12106,
82 
		34.6666,
83 
		93.1206,
84 
	};	/*C/2002 C1 (Ikeya-Zhang) */
85 
 
86 
	ecc = elem.e;
87 
	if(ecc > MAXE)
88 
		ecc = MAXE;
89 
	incl = elem.i * radian;
90 
	node = (elem.o + 0.4593) * radian;
91 
	argp = (elem.w + elem.o + 0.4066) * radian;
92 
	mrad = elem.q / (1-ecc);
93 
        motion = .01720209895 * sqrt(1/(mrad*mrad*mrad))/radian;
94 
	anom = (eday - (elem.t - 2415020)) * motion * radian;
95 
	enom = anom + ecc*sin(anom);
96 
 
97 
	do {
98 
		dele = (anom - enom + ecc * sin(enom)) /
99 
			(1 - ecc*cos(enom));
100 
		enom += dele;
101 
	} while(fabs(dele) > converge);
102 
 
103 
	vnom = 2*atan2(
104 
		sqrt((1+ecc)/(1-ecc))*sin(enom/2),
105 
		cos(enom/2));
106 
	rad = mrad*(1-ecc*cos(enom));
107 
	lambda = vnom + argp;
108 
	pturbl = 0;
109 
	lambda += pturbl*radsec;
110 
	pturbb = 0;
111 
	pturbr = 0;
112 
 
113 
/*
114 
 *	reduce to the ecliptic
115 
 */
116 
	nd = lambda - node;
117 
	lambda = node + atan2(sin(nd)*cos(incl),cos(nd));
118 
 
119 
	sl = sin(incl)*sin(nd) + pturbb*radsec;
120 
	beta = atan2(sl, sqrt(1-sl*sl));
121 
 
122 
	lograd = pturbr*2.30258509;
123 
	rad *= 1 + lograd;
124 
 
125 
	motion *= radian*mrad*mrad/(rad*rad);
126 
	semi = 0;
127 
 
128 
	mag = 5.47 + 6.1/2.303*log(rad);
129 
 
130 
	helio();
131 
	geo();
132 
}