1 cd5bae78 2004-04-21 devnull #include "astro.h"
  2 cd5bae78 2004-04-21 devnull 
  3 cd5bae78 2004-04-21 devnull double k1, k2, k3, k4;
  4 cd5bae78 2004-04-21 devnull double mnom, msun, noded, dmoon;
  5 cd5bae78 2004-04-21 devnull 
  6 cd5bae78 2004-04-21 devnull void
  7 cd5bae78 2004-04-21 devnull moon(void)
  8 cd5bae78 2004-04-21 devnull {
  9 cd5bae78 2004-04-21 devnull 	Moontab *mp;
 10 cd5bae78 2004-04-21 devnull 	double dlong, lsun, psun;
 11 cd5bae78 2004-04-21 devnull 	double eccm, eccs, chp, cpe;
 12 cd5bae78 2004-04-21 devnull 	double v0, t0, m0, j0;
 13 cd5bae78 2004-04-21 devnull 	double arg1, arg2, arg3, arg4, arg5, arg6, arg7;
 14 cd5bae78 2004-04-21 devnull 	double arg8, arg9, arg10;
 15 cd5bae78 2004-04-21 devnull 	double dgamma, k5, k6;
 16 cd5bae78 2004-04-21 devnull 	double lterms, sterms, cterms, nterms, pterms, spterms;
 17 cd5bae78 2004-04-21 devnull 	double gamma1, gamma2, gamma3, arglat;
 18 cd5bae78 2004-04-21 devnull 	double xmp, ymp, zmp;
 19 cd5bae78 2004-04-21 devnull 	double obl2;
 20 cd5bae78 2004-04-21 devnull 
 21 cd5bae78 2004-04-21 devnull /*
 22 cd5bae78 2004-04-21 devnull  *	the fundamental elements - all referred to the epoch of
 23 cd5bae78 2004-04-21 devnull  *	Jan 0.5, 1900 and to the mean equinox of date.
 24 cd5bae78 2004-04-21 devnull  */
 25 cd5bae78 2004-04-21 devnull 
 26 cd5bae78 2004-04-21 devnull 	dlong = 270.434164 + 13.1763965268*eday - .001133*capt2
 27 cd5bae78 2004-04-21 devnull 		 + 2.e-6*capt3;
 28 cd5bae78 2004-04-21 devnull 	argp = 334.329556 + .1114040803*eday - .010325*capt2
 29 cd5bae78 2004-04-21 devnull 		 - 12.e-6*capt3;
 30 cd5bae78 2004-04-21 devnull 	node = 259.183275 - .0529539222*eday + .002078*capt2
 31 cd5bae78 2004-04-21 devnull 		 + 2.e-6*capt3;
 32 cd5bae78 2004-04-21 devnull 	lsun = 279.696678 + .9856473354*eday + .000303*capt2;
 33 cd5bae78 2004-04-21 devnull 	psun = 281.220833 + .0000470684*eday + .000453*capt2
 34 cd5bae78 2004-04-21 devnull 		 + 3.e-6*capt3;
 35 cd5bae78 2004-04-21 devnull 
 36 cd5bae78 2004-04-21 devnull 	dlong = fmod(dlong, 360.);
 37 cd5bae78 2004-04-21 devnull 	argp = fmod(argp, 360.);
 38 cd5bae78 2004-04-21 devnull 	node = fmod(node, 360.);
 39 cd5bae78 2004-04-21 devnull 	lsun = fmod(lsun, 360.);
 40 cd5bae78 2004-04-21 devnull 	psun = fmod(psun, 360.);
 41 cd5bae78 2004-04-21 devnull 
 42 cd5bae78 2004-04-21 devnull 	eccm = 22639.550;
 43 cd5bae78 2004-04-21 devnull 	eccs = .01675104 - .00004180*capt;
 44 cd5bae78 2004-04-21 devnull 	incl = 18461.400;
 45 cd5bae78 2004-04-21 devnull 	cpe = 124.986;
 46 cd5bae78 2004-04-21 devnull 	chp = 3422.451;
 47 cd5bae78 2004-04-21 devnull 
 48 cd5bae78 2004-04-21 devnull /*
 49 cd5bae78 2004-04-21 devnull  *	some subsidiary elements - they are all longitudes
 50 cd5bae78 2004-04-21 devnull  *	and they are referred to the epoch 1/0.5 1900 and
 51 cd5bae78 2004-04-21 devnull  *	to the fixed mean equinox of 1850.0.
 52 cd5bae78 2004-04-21 devnull  */
 53 cd5bae78 2004-04-21 devnull 
 54 cd5bae78 2004-04-21 devnull 	v0 = 342.069128 + 1.6021304820*eday;
 55 cd5bae78 2004-04-21 devnull 	t0 =  98.998753 + 0.9856091138*eday;
 56 cd5bae78 2004-04-21 devnull 	m0 = 293.049675 + 0.5240329445*eday;
 57 cd5bae78 2004-04-21 devnull 	j0 = 237.352319 + 0.0830912295*eday;
 58 cd5bae78 2004-04-21 devnull 
 59 cd5bae78 2004-04-21 devnull /*
 60 cd5bae78 2004-04-21 devnull  *	the following are periodic corrections to the
 61 cd5bae78 2004-04-21 devnull  *	fundamental elements and constants.
 62 cd5bae78 2004-04-21 devnull  *	arg3 is the "Great Venus Inequality".
 63 cd5bae78 2004-04-21 devnull  */
 64 cd5bae78 2004-04-21 devnull 
 65 cd5bae78 2004-04-21 devnull 	arg1 = 41.1 + 20.2*(capt+.5);
 66 cd5bae78 2004-04-21 devnull 	arg2 = dlong - argp + 33. + 3.*t0 - 10.*v0 - 2.6*(capt+.5);
 67 cd5bae78 2004-04-21 devnull 	arg3 = dlong - argp + 151.1 + 16.*t0 - 18.*v0 - (capt+.5);
 68 cd5bae78 2004-04-21 devnull 	arg4 = node;
 69 cd5bae78 2004-04-21 devnull 	arg5 = node + 276.2 - 2.3*(capt+.5);
 70 cd5bae78 2004-04-21 devnull 	arg6 = 313.9 + 13.*t0 - 8.*v0;
 71 cd5bae78 2004-04-21 devnull 	arg7 = dlong - argp + 112.0 + 29.*t0 - 26.*v0;
 72 cd5bae78 2004-04-21 devnull 	arg8 = dlong + argp - 2.*lsun + 273. + 21.*t0 - 20.*v0;
 73 cd5bae78 2004-04-21 devnull 	arg9 = node + 290.1 - 0.9*(capt+.5);
 74 cd5bae78 2004-04-21 devnull 	arg10 = 115. + 38.5*(capt+.5);
 75 cd5bae78 2004-04-21 devnull 	arg1 *= radian;
 76 cd5bae78 2004-04-21 devnull 	arg2 *= radian;
 77 cd5bae78 2004-04-21 devnull 	arg3 *= radian;
 78 cd5bae78 2004-04-21 devnull 	arg4 *= radian;
 79 cd5bae78 2004-04-21 devnull 	arg5 *= radian;
 80 cd5bae78 2004-04-21 devnull 	arg6 *= radian;
 81 cd5bae78 2004-04-21 devnull 	arg7 *= radian;
 82 cd5bae78 2004-04-21 devnull 	arg8 *= radian;
 83 cd5bae78 2004-04-21 devnull 	arg9 *= radian;
 84 cd5bae78 2004-04-21 devnull 	arg10 *= radian;
 85 cd5bae78 2004-04-21 devnull 
 86 cd5bae78 2004-04-21 devnull 	dlong +=
 87 cd5bae78 2004-04-21 devnull 		   (0.84 *sin(arg1)
 88 cd5bae78 2004-04-21 devnull 		 +  0.31 *sin(arg2)
 89 cd5bae78 2004-04-21 devnull 		 + 14.27 *sin(arg3)
 90 cd5bae78 2004-04-21 devnull 		 +  7.261*sin(arg4)
 91 cd5bae78 2004-04-21 devnull 		 +  0.282*sin(arg5)
 92 cd5bae78 2004-04-21 devnull 		 +  0.237*sin(arg6)
 93 cd5bae78 2004-04-21 devnull 		 +  0.108*sin(arg7)
 94 cd5bae78 2004-04-21 devnull 		 +  0.126*sin(arg8))/3600.;
 95 cd5bae78 2004-04-21 devnull 
 96 cd5bae78 2004-04-21 devnull 	argp +=
 97 cd5bae78 2004-04-21 devnull 		 (- 2.10 *sin(arg1)
 98 cd5bae78 2004-04-21 devnull 		 -  0.118*sin(arg3)
 99 cd5bae78 2004-04-21 devnull 		 -  2.076*sin(arg4)
100 cd5bae78 2004-04-21 devnull 		 -  0.840*sin(arg5)
101 cd5bae78 2004-04-21 devnull 		 -  0.593*sin(arg6))/3600.;
102 cd5bae78 2004-04-21 devnull 
103 cd5bae78 2004-04-21 devnull 	node +=
104 cd5bae78 2004-04-21 devnull 		   (0.63*sin(arg1)
105 cd5bae78 2004-04-21 devnull 		 +  0.17*sin(arg3)
106 cd5bae78 2004-04-21 devnull 		 + 95.96*sin(arg4)
107 cd5bae78 2004-04-21 devnull 		 + 15.58*sin(arg5)
108 cd5bae78 2004-04-21 devnull 		 +  1.86*sin(arg9))/3600.;
109 cd5bae78 2004-04-21 devnull 
110 cd5bae78 2004-04-21 devnull 	t0 +=
111 cd5bae78 2004-04-21 devnull 		 (- 6.40*sin(arg1)
112 cd5bae78 2004-04-21 devnull 		 -  1.89*sin(arg6))/3600.;
113 cd5bae78 2004-04-21 devnull 
114 cd5bae78 2004-04-21 devnull 	psun +=
115 cd5bae78 2004-04-21 devnull 		   (6.40*sin(arg1)
116 cd5bae78 2004-04-21 devnull 		 +  1.89*sin(arg6))/3600.;
117 cd5bae78 2004-04-21 devnull 
118 cd5bae78 2004-04-21 devnull 	dgamma = -  4.318*cos(arg4)
119 cd5bae78 2004-04-21 devnull 		 -  0.698*cos(arg5)
120 cd5bae78 2004-04-21 devnull 		 -  0.083*cos(arg9);
121 cd5bae78 2004-04-21 devnull 
122 cd5bae78 2004-04-21 devnull 	j0 +=
123 cd5bae78 2004-04-21 devnull 		   0.33*sin(arg10);
124 cd5bae78 2004-04-21 devnull 
125 cd5bae78 2004-04-21 devnull /*
126 cd5bae78 2004-04-21 devnull  *	the following factors account for the fact that the
127 cd5bae78 2004-04-21 devnull  *	eccentricity, solar eccentricity, inclination and
128 cd5bae78 2004-04-21 devnull  *	parallax used by Brown to make up his coefficients
129 cd5bae78 2004-04-21 devnull  *	are both wrong and out of date.  Brown did the same
130 cd5bae78 2004-04-21 devnull  *	thing in a different way.
131 cd5bae78 2004-04-21 devnull  */
132 cd5bae78 2004-04-21 devnull 
133 cd5bae78 2004-04-21 devnull 	k1 = eccm/22639.500;
134 cd5bae78 2004-04-21 devnull 	k2 = eccs/.01675104;
135 cd5bae78 2004-04-21 devnull 	k3 = 1. + 2.708e-6 + .000108008*dgamma;
136 cd5bae78 2004-04-21 devnull 	k4 = cpe/125.154;
137 cd5bae78 2004-04-21 devnull 	k5 = chp/3422.700;
138 cd5bae78 2004-04-21 devnull 
139 cd5bae78 2004-04-21 devnull /*
140 cd5bae78 2004-04-21 devnull  *	the principal arguments that are used to compute
141 cd5bae78 2004-04-21 devnull  *	perturbations are the following differences of the
142 cd5bae78 2004-04-21 devnull  *	fundamental elements.
143 cd5bae78 2004-04-21 devnull  */
144 cd5bae78 2004-04-21 devnull 
145 cd5bae78 2004-04-21 devnull 	mnom = dlong - argp;
146 cd5bae78 2004-04-21 devnull 	msun = lsun - psun;
147 cd5bae78 2004-04-21 devnull 	noded = dlong - node;
148 cd5bae78 2004-04-21 devnull 	dmoon = dlong - lsun;
149 cd5bae78 2004-04-21 devnull 
150 cd5bae78 2004-04-21 devnull /*
151 cd5bae78 2004-04-21 devnull  *	solar terms in longitude
152 cd5bae78 2004-04-21 devnull  */
153 cd5bae78 2004-04-21 devnull 
154 cd5bae78 2004-04-21 devnull 	lterms = 0.0;
155 cd5bae78 2004-04-21 devnull 	mp = moontab;
156 cd5bae78 2004-04-21 devnull 	for(;;) {
157 cd5bae78 2004-04-21 devnull 		if(mp->f == 0.0)
158 cd5bae78 2004-04-21 devnull 			break;
159 cd5bae78 2004-04-21 devnull 		lterms += sinx(mp->f,
160 cd5bae78 2004-04-21 devnull 			mp->c[0], mp->c[1],
161 cd5bae78 2004-04-21 devnull 			mp->c[2], mp->c[3], 0.0);
162 cd5bae78 2004-04-21 devnull 		mp++;
163 cd5bae78 2004-04-21 devnull 	}
164 cd5bae78 2004-04-21 devnull 	mp++;
165 cd5bae78 2004-04-21 devnull 
166 cd5bae78 2004-04-21 devnull /*
167 cd5bae78 2004-04-21 devnull  *	planetary terms in longitude
168 cd5bae78 2004-04-21 devnull  */
169 cd5bae78 2004-04-21 devnull 
170 cd5bae78 2004-04-21 devnull 	lterms += sinx(0.822, 0,0,0,0, t0-v0);
171 cd5bae78 2004-04-21 devnull 	lterms += sinx(0.307, 0,0,0,0, 2.*t0-2.*v0+179.8);
172 cd5bae78 2004-04-21 devnull 	lterms += sinx(0.348, 0,0,0,0, 3.*t0-2.*v0+272.9);
173 cd5bae78 2004-04-21 devnull 	lterms += sinx(0.176, 0,0,0,0, 4.*t0-3.*v0+271.7);
174 cd5bae78 2004-04-21 devnull 	lterms += sinx(0.092, 0,0,0,0, 5.*t0-3.*v0+199.);
175 cd5bae78 2004-04-21 devnull 	lterms += sinx(0.129, 1,0,0,0, -t0+v0+180.);
176 cd5bae78 2004-04-21 devnull 	lterms += sinx(0.152, 1,0,0,0, t0-v0);
177 cd5bae78 2004-04-21 devnull 	lterms += sinx(0.127, 1,0,0,0, 3.*t0-3.*v0+180.);
178 cd5bae78 2004-04-21 devnull 	lterms += sinx(0.099, 0,0,0,2, t0-v0);
179 cd5bae78 2004-04-21 devnull 	lterms += sinx(0.136, 0,0,0,2, 2.*t0-2.*v0+179.5);
180 cd5bae78 2004-04-21 devnull 	lterms += sinx(0.083, -1,0,0,2, -4.*t0+4.*v0+180.);
181 cd5bae78 2004-04-21 devnull 	lterms += sinx(0.662, -1,0,0,2, -3.*t0+3.*v0+180.0);
182 cd5bae78 2004-04-21 devnull 	lterms += sinx(0.137, -1,0,0,2, -2.*t0+2.*v0);
183 cd5bae78 2004-04-21 devnull 	lterms += sinx(0.133, -1,0,0,2, t0-v0);
184 cd5bae78 2004-04-21 devnull 	lterms += sinx(0.157, -1,0,0,2, 2.*t0-2.*v0+179.6);
185 cd5bae78 2004-04-21 devnull 	lterms += sinx(0.079, -1,0,0,2, -8.*t0+6.*v0+162.6);
186 cd5bae78 2004-04-21 devnull 	lterms += sinx(0.073, 2,0,0,-2, 3.*t0-3.*v0+180.);
187 cd5bae78 2004-04-21 devnull 	lterms += sinx(0.643, 0,0,0,0, -t0+j0+178.8);
188 cd5bae78 2004-04-21 devnull 	lterms += sinx(0.187, 0,0,0,0, -2.*t0+2.*j0+359.6);
189 cd5bae78 2004-04-21 devnull 	lterms += sinx(0.087, 0,0,0,0, j0+289.9);
190 cd5bae78 2004-04-21 devnull 	lterms += sinx(0.165, 0,0,0,0, -t0+2.*j0+241.5);
191 cd5bae78 2004-04-21 devnull 	lterms += sinx(0.144, 1,0,0,0, t0-j0+1.0);
192 cd5bae78 2004-04-21 devnull 	lterms += sinx(0.158, 1,0,0,0, -t0+j0+179.0);
193 cd5bae78 2004-04-21 devnull 	lterms += sinx(0.190, 1,0,0,0, -2.*t0+2.*j0+180.0);
194 cd5bae78 2004-04-21 devnull 	lterms += sinx(0.096, 1,0,0,0, -2.*t0+3.*j0+352.5);
195 cd5bae78 2004-04-21 devnull 	lterms += sinx(0.070, 0,0,0,2, 2.*t0-2.*j0+180.);
196 cd5bae78 2004-04-21 devnull 	lterms += sinx(0.167, 0,0,0,2, -t0+j0+178.5);
197 cd5bae78 2004-04-21 devnull 	lterms += sinx(0.085, 0,0,0,2, -2.*t0+2.*j0+359.2);
198 cd5bae78 2004-04-21 devnull 	lterms += sinx(1.137, -1,0,0,2, 2.*t0-2.*j0+180.3);
199 cd5bae78 2004-04-21 devnull 	lterms += sinx(0.211, -1,0,0,2, -t0+j0+178.4);
200 cd5bae78 2004-04-21 devnull 	lterms += sinx(0.089, -1,0,0,2, -2.*t0+2.*j0+359.2);
201 cd5bae78 2004-04-21 devnull 	lterms += sinx(0.436, -1,0,0,2, 2.*t0-3.*j0+7.5);
202 cd5bae78 2004-04-21 devnull 	lterms += sinx(0.240, 2,0,0,-2, -2.*t0+2.*j0+179.9);
203 cd5bae78 2004-04-21 devnull 	lterms += sinx(0.284, 2,0,0,-2, -2.*t0+3.*j0+172.5);
204 cd5bae78 2004-04-21 devnull 	lterms += sinx(0.195, 0,0,0,0, -2.*t0+2.*m0+180.2);
205 cd5bae78 2004-04-21 devnull 	lterms += sinx(0.327, 0,0,0,0, -t0+2.*m0+224.4);
206 cd5bae78 2004-04-21 devnull 	lterms += sinx(0.093, 0,0,0,0, -2.*t0+4.*m0+244.8);
207 cd5bae78 2004-04-21 devnull 	lterms += sinx(0.073, 1,0,0,0, -t0+2.*m0+223.3);
208 cd5bae78 2004-04-21 devnull 	lterms += sinx(0.074, 1,0,0,0, t0-2.*m0+306.3);
209 cd5bae78 2004-04-21 devnull 	lterms += sinx(0.189, 0,0,0,0, node+180.);
210 cd5bae78 2004-04-21 devnull 
211 cd5bae78 2004-04-21 devnull /*
212 cd5bae78 2004-04-21 devnull  *	solar terms in latitude
213 cd5bae78 2004-04-21 devnull  */
214 cd5bae78 2004-04-21 devnull 
215 cd5bae78 2004-04-21 devnull 	sterms = 0;
216 cd5bae78 2004-04-21 devnull 	for(;;) {
217 cd5bae78 2004-04-21 devnull 		if(mp->f == 0)
218 cd5bae78 2004-04-21 devnull 			break;
219 cd5bae78 2004-04-21 devnull 		sterms += sinx(mp->f,
220 cd5bae78 2004-04-21 devnull 			mp->c[0], mp->c[1],
221 cd5bae78 2004-04-21 devnull 			mp->c[2], mp->c[3], 0);
222 cd5bae78 2004-04-21 devnull 		mp++;
223 cd5bae78 2004-04-21 devnull 	}
224 cd5bae78 2004-04-21 devnull 	mp++;
225 cd5bae78 2004-04-21 devnull 
226 cd5bae78 2004-04-21 devnull 	cterms = 0;
227 cd5bae78 2004-04-21 devnull 	for(;;) {
228 cd5bae78 2004-04-21 devnull 		if(mp->f == 0)
229 cd5bae78 2004-04-21 devnull 			break;
230 cd5bae78 2004-04-21 devnull 		cterms += cosx(mp->f,
231 cd5bae78 2004-04-21 devnull 			mp->c[0], mp->c[1],
232 cd5bae78 2004-04-21 devnull 			mp->c[2], mp->c[3], 0);
233 cd5bae78 2004-04-21 devnull 		mp++;
234 cd5bae78 2004-04-21 devnull 	}
235 cd5bae78 2004-04-21 devnull 	mp++;
236 cd5bae78 2004-04-21 devnull 
237 cd5bae78 2004-04-21 devnull 	nterms = 0;
238 cd5bae78 2004-04-21 devnull 	for(;;) {
239 cd5bae78 2004-04-21 devnull 		if(mp->f == 0)
240 cd5bae78 2004-04-21 devnull 			break;
241 cd5bae78 2004-04-21 devnull 		nterms += sinx(mp->f,
242 cd5bae78 2004-04-21 devnull 			mp->c[0], mp->c[1],
243 cd5bae78 2004-04-21 devnull 			mp->c[2], mp->c[3], 0);
244 cd5bae78 2004-04-21 devnull 		mp++;
245 cd5bae78 2004-04-21 devnull 	}
246 cd5bae78 2004-04-21 devnull 	mp++;
247 cd5bae78 2004-04-21 devnull 
248 cd5bae78 2004-04-21 devnull /*
249 cd5bae78 2004-04-21 devnull  *	planetary terms in latitude
250 cd5bae78 2004-04-21 devnull  */
251 cd5bae78 2004-04-21 devnull 
252 cd5bae78 2004-04-21 devnull 	pterms =
253 cd5bae78 2004-04-21 devnull 		   sinx(0.215, 0,0,0,0, dlong);
254 cd5bae78 2004-04-21 devnull 
255 cd5bae78 2004-04-21 devnull /*
256 cd5bae78 2004-04-21 devnull  *	solar terms in parallax
257 cd5bae78 2004-04-21 devnull  */
258 cd5bae78 2004-04-21 devnull 
259 cd5bae78 2004-04-21 devnull 	spterms = 3422.700;
260 cd5bae78 2004-04-21 devnull 	for(;;) {
261 cd5bae78 2004-04-21 devnull 		if(mp->f == 0)
262 cd5bae78 2004-04-21 devnull 			break;
263 cd5bae78 2004-04-21 devnull 		spterms += cosx(mp->f,
264 cd5bae78 2004-04-21 devnull 			mp->c[0], mp->c[1],
265 cd5bae78 2004-04-21 devnull 			mp->c[2], mp->c[3], 0);
266 cd5bae78 2004-04-21 devnull 		mp++;
267 cd5bae78 2004-04-21 devnull 	}
268 cd5bae78 2004-04-21 devnull 
269 cd5bae78 2004-04-21 devnull /*
270 cd5bae78 2004-04-21 devnull  *	planetary terms in parallax
271 cd5bae78 2004-04-21 devnull  */
272 cd5bae78 2004-04-21 devnull 
273 c5bfba48 2012-10-21 rsc 	//spterms = spterms;
274 cd5bae78 2004-04-21 devnull 
275 cd5bae78 2004-04-21 devnull /*
276 cd5bae78 2004-04-21 devnull  *	computation of longitude
277 cd5bae78 2004-04-21 devnull  */
278 cd5bae78 2004-04-21 devnull 
279 cd5bae78 2004-04-21 devnull 	lambda = (dlong + lterms/3600.)*radian;
280 cd5bae78 2004-04-21 devnull 
281 cd5bae78 2004-04-21 devnull /*
282 cd5bae78 2004-04-21 devnull  *	computation of latitude
283 cd5bae78 2004-04-21 devnull  */
284 cd5bae78 2004-04-21 devnull 
285 cd5bae78 2004-04-21 devnull 	arglat = (noded + sterms/3600.)*radian;
286 cd5bae78 2004-04-21 devnull 	gamma1 = 18519.700 * k3;
287 cd5bae78 2004-04-21 devnull 	gamma2 = -6.241 * k3*k3*k3;
288 cd5bae78 2004-04-21 devnull 	gamma3 = 0.004 * k3*k3*k3*k3*k3;
289 cd5bae78 2004-04-21 devnull 
290 cd5bae78 2004-04-21 devnull 	k6 = (gamma1 + cterms) / gamma1;
291 cd5bae78 2004-04-21 devnull 
292 cd5bae78 2004-04-21 devnull 	beta = k6 * (gamma1*sin(arglat) + gamma2*sin(3.*arglat)
293 cd5bae78 2004-04-21 devnull 		 + gamma3*sin(5.*arglat) + nterms)
294 cd5bae78 2004-04-21 devnull 		 + pterms;
295 cd5bae78 2004-04-21 devnull 	if(flags['o'])
296 cd5bae78 2004-04-21 devnull 		beta -= 0.6;
297 cd5bae78 2004-04-21 devnull 	beta *= radsec;
298 cd5bae78 2004-04-21 devnull 
299 cd5bae78 2004-04-21 devnull /*
300 cd5bae78 2004-04-21 devnull  *	computation of parallax
301 cd5bae78 2004-04-21 devnull  */
302 cd5bae78 2004-04-21 devnull 
303 cd5bae78 2004-04-21 devnull 	spterms = k5 * spterms *radsec;
304 cd5bae78 2004-04-21 devnull 	hp = spterms + (spterms*spterms*spterms)/6.;
305 cd5bae78 2004-04-21 devnull 
306 cd5bae78 2004-04-21 devnull 	rad = hp/radsec;
307 cd5bae78 2004-04-21 devnull 	rp = 1.;
308 cd5bae78 2004-04-21 devnull 	semi = .0799 + .272453*(hp/radsec);
309 cd5bae78 2004-04-21 devnull 	if(dmoon < 0.)
310 cd5bae78 2004-04-21 devnull 		dmoon += 360.;
311 cd5bae78 2004-04-21 devnull 	mag = dmoon/360.;
312 cd5bae78 2004-04-21 devnull 
313 cd5bae78 2004-04-21 devnull /*
314 cd5bae78 2004-04-21 devnull  *	change to equatorial coordinates
315 cd5bae78 2004-04-21 devnull  */
316 cd5bae78 2004-04-21 devnull 
317 cd5bae78 2004-04-21 devnull 	lambda += phi;
318 cd5bae78 2004-04-21 devnull 	obl2 = obliq + eps;
319 cd5bae78 2004-04-21 devnull 	xmp = rp*cos(lambda)*cos(beta);
320 cd5bae78 2004-04-21 devnull 	ymp = rp*(sin(lambda)*cos(beta)*cos(obl2) - sin(obl2)*sin(beta));
321 cd5bae78 2004-04-21 devnull 	zmp = rp*(sin(lambda)*cos(beta)*sin(obl2) + cos(obl2)*sin(beta));
322 cd5bae78 2004-04-21 devnull 
323 cd5bae78 2004-04-21 devnull 	alpha = atan2(ymp, xmp);
324 cd5bae78 2004-04-21 devnull 	delta = atan2(zmp, sqrt(xmp*xmp+ymp*ymp));
325 cd5bae78 2004-04-21 devnull 	meday = eday;
326 cd5bae78 2004-04-21 devnull 	mhp = hp;
327 cd5bae78 2004-04-21 devnull 
328 cd5bae78 2004-04-21 devnull 	geo();
329 cd5bae78 2004-04-21 devnull }
330 cd5bae78 2004-04-21 devnull 
331 cd5bae78 2004-04-21 devnull double
332 cd5bae78 2004-04-21 devnull sinx(double coef, int i, int j, int k, int m, double angle)
333 cd5bae78 2004-04-21 devnull {
334 cd5bae78 2004-04-21 devnull 	double x;
335 cd5bae78 2004-04-21 devnull 
336 cd5bae78 2004-04-21 devnull 	x = i*mnom + j*msun + k*noded + m*dmoon + angle;
337 cd5bae78 2004-04-21 devnull 	x = coef*sin(x*radian);
338 cd5bae78 2004-04-21 devnull 	if(i < 0)
339 cd5bae78 2004-04-21 devnull 		i = -i;
340 cd5bae78 2004-04-21 devnull 	for(; i>0; i--)
341 cd5bae78 2004-04-21 devnull 		x *= k1;
342 cd5bae78 2004-04-21 devnull 	if(j < 0)
343 cd5bae78 2004-04-21 devnull 		j = -j;
344 cd5bae78 2004-04-21 devnull 	for(; j>0; j--)
345 cd5bae78 2004-04-21 devnull 		x *= k2;
346 cd5bae78 2004-04-21 devnull 	if(k < 0)
347 cd5bae78 2004-04-21 devnull 		k = -k;
348 cd5bae78 2004-04-21 devnull 	for(; k>0; k--)
349 cd5bae78 2004-04-21 devnull 		x *= k3;
350 cd5bae78 2004-04-21 devnull 	if(m & 1)
351 cd5bae78 2004-04-21 devnull 		x *= k4;
352 cd5bae78 2004-04-21 devnull 
353 cd5bae78 2004-04-21 devnull 	return x;
354 cd5bae78 2004-04-21 devnull }
355 cd5bae78 2004-04-21 devnull 
356 cd5bae78 2004-04-21 devnull double
357 cd5bae78 2004-04-21 devnull cosx(double coef, int i, int j, int k, int m, double angle)
358 cd5bae78 2004-04-21 devnull {
359 cd5bae78 2004-04-21 devnull 	double x;
360 cd5bae78 2004-04-21 devnull 
361 cd5bae78 2004-04-21 devnull 	x = i*mnom + j*msun + k*noded + m*dmoon + angle;
362 cd5bae78 2004-04-21 devnull 	x = coef*cos(x*radian);
363 cd5bae78 2004-04-21 devnull 	if(i < 0)
364 cd5bae78 2004-04-21 devnull 		i = -i;
365 cd5bae78 2004-04-21 devnull 	for(; i>0; i--)
366 cd5bae78 2004-04-21 devnull 		x *= k1;
367 cd5bae78 2004-04-21 devnull 	if(j < 0)
368 cd5bae78 2004-04-21 devnull 		j = -j;
369 cd5bae78 2004-04-21 devnull 	for(; j>0; j--)
370 cd5bae78 2004-04-21 devnull 		x *= k2;
371 cd5bae78 2004-04-21 devnull 	if(k < 0)
372 cd5bae78 2004-04-21 devnull 		k = -k;
373 cd5bae78 2004-04-21 devnull 	for(; k>0; k--)
374 cd5bae78 2004-04-21 devnull 		x *= k3;
375 cd5bae78 2004-04-21 devnull 	if(m & 1)
376 cd5bae78 2004-04-21 devnull 		x *= k4;
377 cd5bae78 2004-04-21 devnull 
378 cd5bae78 2004-04-21 devnull 	return x;
379 cd5bae78 2004-04-21 devnull }