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