source: trunk/MagicSoft/slalib/evp.c@ 8098

Last change on this file since 8098 was 731, checked in by tbretz, 25 years ago
*** empty log message ***
  • Property svn:executable set to *
File size: 22.4 KB
Line 
1#include "slalib.h"
2#include "slamac.h"
3void slaEvp ( double date, double deqx, double dvb[3], double dpb[3],
4 double dvh[3], double dph[3] )
5/*
6** - - - - - - -
7** s l a E v p
8** - - - - - - -
9**
10** Barycentric and heliocentric velocity and position of the Earth.
11**
12** Given:
13**
14** date double TDB (loosely ET) as a Modified Julian Date
15** (JD-2400000.5)
16**
17** deqx double Julian epoch (e.g. 2000.0) of mean equator and
18** equinox of the vectors returned. If deqx <= 0.0,
19** all vectors are referred to the mean equator and
20** equinox (FK5) of epoch date.
21**
22** Returned (all 3D Cartesian vectors):
23**
24** dvb,dpb double[3] barycentric velocity, position
25**
26** dvh,dph double[3] heliocentric velocity, position
27**
28** (Units are AU/s for velocity and AU for position)
29**
30** Called: slaEpj, slaPrec
31**
32** Accuracy:
33**
34** The maximum deviations from the JPL DE96 ephemeris are as
35** follows:
36**
37** barycentric velocity 42 cm/s
38** barycentric position 6900 km
39**
40** heliocentric velocity 42 cm/s
41** heliocentric position 1600 km
42**
43** This routine is adapted from the BARVEL and BARCOR Fortran
44** subroutines of P.Stumpff, which are described in
45** Astron. Astrophys. Suppl. Ser. 41, 1-8 (1980). The present
46** routine uses double precision throughout; most of the other
47** changes are essentially cosmetic and do not affect the
48** results. However, some adjustments have been made so as to
49** give results that refer to the new (IAU 1976 "FK5") equinox
50** and precession, although the differences these changes make
51** relative to the results from Stumpff's original "FK4" version
52** are smaller than the inherent accuracy of the algorithm. One
53** minor shortcoming in the original routines that has not been
54** corrected is that better numerical accuracy could be achieved
55** if the various polynomial evaluations were nested. Note also
56** that one of Stumpff's precession constants differs by 0.001 arcsec
57** from the value given in the Explanatory Supplement to the A.E.
58**
59** Defined in slamac.h: D2PI, DS2R, dmod
60**
61** Last revision: 21 March 1999
62**
63** Copyright P.T.Wallace. All rights reserved.
64*/
65{
66 int ideq, i, j, k;
67
68 double a, pertl,
69 pertld, pertr, pertrd, cosa, sina, e, twoe, esq, g, twog,
70 phi, f, sf, cf, phid, psid, pertp, pertpd, tl, sinlm, coslm,
71 sigma, b, plon, pomg, pecc, flatm, flat;
72
73 double dt, dlocal, dml,
74 deps, dparam, dpsi, d1pdro, drd, drld, dtl, dsinls,
75 dcosls, dxhd, dyhd, dzhd, dxbd, dybd, dzbd, dcosep,
76 dsinep, dyahd, dzahd, dyabd, dzabd, dr,
77 dxh, dyh, dzh, dxb, dyb, dzb, dyah, dzah, dyab,
78 dzab, depj, deqcor;
79
80 double sn[4], forbel[7], sorbel[17], sinlp[4], coslp[4];
81
82 double dprema[3][3], w, vw[3];
83
84/* Sidereal rate dcsld in longitude, rate ccsgd in mean anomaly */
85 static double dcsld = 1.990987e-7;
86 static double ccsgd = 1.990969e-7;
87
88/* Some constants used in the calculation of the lunar contribution */
89 static double cckm = 3.122140e-5;
90 static double ccmld = 2.661699e-6;
91 static double ccfdi = 2.399485e-7;
92
93/* Besselian epoch 1950.0 expressed as a Julian epoch */
94 static double b1950 = 1949.9997904423;
95
96/*
97** ccpamv(k)=a*m*dl/dt (planets), dc1mme=1-mass(Earth+Moon)
98*/
99 static double ccpamv[4] = {
100 8.326827e-11,
101 1.843484e-11,
102 1.988712e-12,
103 1.881276e-12
104 };
105 static double dc1mme = 0.99999696;
106
107/*
108** ccpam(k)=a*m(planets)
109** ccim=inclination(Moon)
110*/
111 static double ccpam[4] = {
112 4.960906e-3,
113 2.727436e-3,
114 8.392311e-4,
115 1.556861e-3
116 };
117 static double ccim = 8.978749e-2;
118
119/*
120** Constants dcfel(i,k) of fast changing elements
121*/
122 static double dcfel[3][8] = {
123 { 1.7400353, /* dcfel[0][0] */
124 6.2565836, /* dcfel[0][1] */
125 4.7199666, /* dcfel[0][2] */
126 1.9636505e-1, /* dcfel[0][3] */
127 4.1547339, /* dcfel[0][4] */
128 4.6524223, /* dcfel[0][5] */
129 4.2620486, /* dcfel[0][6] */
130 1.4740694 }, /* dcfel[0][7] */
131 { 6.2833195099091e+2, /* dcfel[1][0] */
132 6.2830194572674e+2, /* dcfel[1][1] */
133 8.3997091449254e+3, /* dcfel[1][2] */
134 8.4334662911720e+3, /* dcfel[1][3] */
135 5.2993466764997e+1, /* dcfel[1][4] */
136 2.1354275911213e+1, /* dcfel[1][5] */
137 7.5025342197656, /* dcfel[1][6] */
138 3.8377331909193 }, /* dcfel[1][7] */
139 { 5.2796e-6, /* dcfel[2][0] */
140 -2.6180e-6, /* dcfel[2][1] */
141 -1.9780e-5, /* dcfel[2][2] */
142 -5.6044e-5, /* dcfel[2][3] */
143 5.8845e-6, /* dcfel[2][4] */
144 5.6797e-6, /* dcfel[2][5] */
145 5.5317e-6, /* dcfel[2][6] */
146 5.6093e-6 } /* dcfel[2][7] */
147 };
148
149/*
150** Constants dceps and ccsel(i,k) of slowly changing elements
151*/
152 static double dceps[3] = {
153 4.093198e-1,
154 -2.271110e-4,
155 -2.860401e-8
156 };
157 static double ccsel[3][17] = {
158 { 1.675104e-2, /* ccsel[0][0] */
159 2.220221e-1, /* ccsel[0][1] */
160 1.589963, /* ccsel[0][2] */
161 2.994089, /* ccsel[0][3] */
162 8.155457e-1, /* ccsel[0][4] */
163 1.735614, /* ccsel[0][5] */
164 1.968564, /* ccsel[0][6] */
165 1.282417, /* ccsel[0][7] */
166 2.280820, /* ccsel[0][8] */
167 4.833473e-2, /* ccsel[0][9] */
168 5.589232e-2, /* ccsel[0][10] */
169 4.634443e-2, /* ccsel[0][11] */
170 8.997041e-3, /* ccsel[0][12] */
171 2.284178e-2, /* ccsel[0][13] */
172 4.350267e-2, /* ccsel[0][14] */
173 1.348204e-2, /* ccsel[0][15] */
174 3.106570e-2 }, /* ccsel[0][16] */
175 { -4.179579e-5, /* ccsel[1][0] */
176 2.809917e-2, /* ccsel[1][1] */
177 3.418075e-2, /* ccsel[1][2] */
178 2.590824e-2, /* ccsel[1][3] */
179 2.486352e-2, /* ccsel[1][4] */
180 1.763719e-2, /* ccsel[1][5] */
181 1.524020e-2, /* ccsel[1][6] */
182 8.703393e-3, /* ccsel[1][7] */
183 1.918010e-2, /* ccsel[1][8] */
184 1.641773e-4, /* ccsel[1][9] */
185 -3.455092e-4, /* ccsel[1][10] */
186 -2.658234e-5, /* ccsel[1][11] */
187 6.329728e-6, /* ccsel[1][12] */
188 -9.941590e-5, /* ccsel[1][13] */
189 -6.839749e-5, /* ccsel[1][14] */
190 1.091504e-5, /* ccsel[1][15] */
191 -1.665665e-4 }, /* ccsel[1][16] */
192 { -1.260516e-7, /* ccsel[2][0] */
193 1.852532e-5, /* ccsel[2][1] */
194 1.430200e-5, /* ccsel[2][2] */
195 4.155840e-6, /* ccsel[2][3] */
196 6.836840e-6, /* ccsel[2][4] */
197 6.370440e-6, /* ccsel[2][5] */
198 -2.517152e-6, /* ccsel[2][6] */
199 2.289292e-5, /* ccsel[2][7] */
200 4.484520e-6, /* ccsel[2][8] */
201 -4.654200e-7, /* ccsel[2][9] */
202 -7.388560e-7, /* ccsel[2][10] */
203 7.757000e-8, /* ccsel[2][11] */
204 -1.939256e-9, /* ccsel[2][12] */
205 6.787400e-8, /* ccsel[2][13] */
206 -2.714956e-7, /* ccsel[2][14] */
207 6.903760e-7, /* ccsel[2][15] */
208 -1.590188e-7 } /* ccsel[2][16] */
209 };
210
211/*
212** Constants of the arguments of the short-period perturbations
213** by the planets: dcargs(i,k)
214*/
215 static double dcargs[2][15] = {
216 { 5.0974222, /* dcargs[0][0] */
217 3.9584962, /* dcargs[0][1] */
218 1.6338070, /* dcargs[0][2] */
219 2.5487111, /* dcargs[0][3] */
220 4.9255514, /* dcargs[0][4] */
221 1.3363463, /* dcargs[0][5] */
222 1.6072053, /* dcargs[0][6] */
223 1.3629480, /* dcargs[0][7] */
224 5.5657014, /* dcargs[0][8] */
225 5.0708205, /* dcargs[0][9] */
226 3.9318944, /* dcargs[0][10] */
227 4.8989497, /* dcargs[0][11] */
228 1.3097446, /* dcargs[0][12] */
229 3.5147141, /* dcargs[0][13] */
230 3.5413158 }, /* dcargs[0][14] */
231 { -7.8604195454652e+2, /* dcargs[1][0] */
232 -5.7533848094674e+2, /* dcargs[1][1] */
233 -1.1506769618935e+3, /* dcargs[1][2] */
234 -3.9302097727326e+2, /* dcargs[1][3] */
235 -5.8849265665348e+2, /* dcargs[1][4] */
236 -5.5076098609303e+2, /* dcargs[1][5] */
237 -5.2237501616674e+2, /* dcargs[1][6] */
238 -1.1790629318198e+3, /* dcargs[1][7] */
239 -1.0977134971135e+3, /* dcargs[1][8] */
240 -1.5774000881978e+2, /* dcargs[1][9] */
241 5.2963464780000e+1, /* dcargs[1][10] */
242 3.9809289073258e+1, /* dcargs[1][11] */
243 7.7540959633708e+1, /* dcargs[1][12] */
244 7.9618578146517e+1, /* dcargs[1][13] */
245 -5.4868336758022e+2 } /* dcargs[1][14] */
246 };
247
248/*
249** Amplitudes ccamps(n,k) of the short-period perturbations
250*/
251 static double ccamps[5][15] = {
252 { -2.279594e-5, /* ccamps[0][0] */
253 -3.494537e-5, /* ccamps[0][1] */
254 6.593466e-7, /* ccamps[0][2] */
255 1.140767e-5, /* ccamps[0][3] */
256 9.516893e-6, /* ccamps[0][4] */
257 7.310990e-6, /* ccamps[0][5] */
258 -2.603449e-6, /* ccamps[0][6] */
259 -3.228859e-6, /* ccamps[0][7] */
260 3.442177e-7, /* ccamps[0][8] */
261 8.702406e-6, /* ccamps[0][9] */
262 -1.488378e-6, /* ccamps[0][10] */
263 -8.043059e-6, /* ccamps[0][11] */
264 3.699128e-6, /* ccamps[0][12] */
265 2.550120e-6, /* ccamps[0][13] */
266 -6.351059e-7 }, /* ccamps[0][14] */
267 { 1.407414e-5, /* ccamps[1][0] */
268 2.860401e-7, /* ccamps[1][1] */
269 1.322572e-5, /* ccamps[1][2] */
270 -2.049792e-5, /* ccamps[1][3] */
271 -2.748894e-6, /* ccamps[1][4] */
272 -1.924710e-6, /* ccamps[1][5] */
273 7.359472e-6, /* ccamps[1][6] */
274 1.308997e-7, /* ccamps[1][7] */
275 2.671323e-6, /* ccamps[1][8] */
276 -8.421214e-6, /* ccamps[1][9] */
277 -1.251789e-5, /* ccamps[1][10] */
278 -2.991300e-6, /* ccamps[1][11] */
279 -3.316126e-6, /* ccamps[1][12] */
280 -1.241123e-6, /* ccamps[1][13] */
281 2.341650e-6 }, /* ccamps[1][14] */
282 { 8.273188e-6, /* ccamps[2][0] */
283 1.289448e-7, /* ccamps[2][1] */
284 9.258695e-6, /* ccamps[2][2] */
285 -4.747930e-6, /* ccamps[2][3] */
286 -1.319381e-6, /* ccamps[2][4] */
287 -8.772849e-7, /* ccamps[2][5] */
288 3.168357e-6, /* ccamps[2][6] */
289 1.013137e-7, /* ccamps[2][7] */
290 1.832858e-6, /* ccamps[2][8] */
291 -1.372341e-6, /* ccamps[2][9] */
292 5.226868e-7, /* ccamps[2][10] */
293 1.473654e-7, /* ccamps[2][11] */
294 2.901257e-7, /* ccamps[2][12] */
295 9.901116e-8, /* ccamps[2][13] */
296 1.061492e-6 }, /* ccamps[2][14] */
297 { 1.340565e-5, /* ccamps[3][0] */
298 1.627237e-5, /* ccamps[3][1] */
299 -4.674248e-7, /* ccamps[3][2] */
300 -2.638763e-6, /* ccamps[3][3] */
301 -4.549908e-6, /* ccamps[3][4] */
302 -3.334143e-6, /* ccamps[3][5] */
303 1.119056e-6, /* ccamps[3][6] */
304 2.403899e-6, /* ccamps[3][7] */
305 -2.394688e-7, /* ccamps[3][8] */
306 -1.455234e-6, /* ccamps[3][9] */
307 -2.049301e-7, /* ccamps[3][10] */
308 -3.154542e-7, /* ccamps[3][11] */
309 3.407826e-7, /* ccamps[3][12] */
310 2.210482e-7, /* ccamps[3][13] */
311 2.878231e-7 }, /* ccamps[3][14] */
312 { -2.490817e-7, /* ccamps[4][0] */
313 -1.823138e-7, /* ccamps[4][1] */
314 -3.646275e-7, /* ccamps[4][2] */
315 -1.245408e-7, /* ccamps[4][3] */
316 -1.864821e-7, /* ccamps[4][4] */
317 -1.745256e-7, /* ccamps[4][5] */
318 -1.655307e-7, /* ccamps[4][6] */
319 -3.736225e-7, /* ccamps[4][7] */
320 -3.478444e-7, /* ccamps[4][8] */
321 -4.998479e-8, /* ccamps[4][9] */
322 0.0, /* ccamps[4][10] */
323 0.0, /* ccamps[4][11] */
324 0.0, /* ccamps[4][12] */
325 0.0, /* ccamps[4][13] */
326 0.0 } /* ccamps[4][14] */
327 };
328
329/*
330** Constants of the secular perturbations in longitude
331** ccsec3 and ccsec(n,k)
332*/
333 static double ccsec3 = -7.757020e-8;
334 static double ccsec[3][4] = {
335 { 1.289600e-6, /* ccsec[0][0] */
336 3.102810e-5, /* ccsec[0][1] */
337 9.124190e-6, /* ccsec[0][2] */
338 9.793240e-7 }, /* ccsec[0][3] */
339 { 5.550147e-1, /* ccsec[1][0] */
340 4.035027, /* ccsec[1][1] */
341 9.990265e-1, /* ccsec[1][2] */
342 5.508259 }, /* ccsec[1][3] */
343 { 2.076942, /* ccsec[2][0] */
344 3.525565e-1, /* ccsec[2][1] */
345 2.622706, /* ccsec[2][2] */
346 1.559103e+1 } /* ccsec[2][3] */
347 };
348
349/*
350** Constants dcargm(i,k) of the arguments of the perturbations
351** of the motion of the Moon
352*/
353 static double dcargm[2][3] = {
354 { 5.167983, /* dcargm[0][0] */
355 5.491315, /* dcargm[0][1] */
356 5.959853 }, /* dcargm[0][2] */
357 { 8.3286911095275e+3, /* dcargm[1][0] */
358 -7.2140632838100e+3, /* dcargm[1][1] */
359 1.5542754389685e+4 } /* dcargm[1][2] */
360 };
361
362/*
363** Amplitudes ccampm(n,k) of the perturbations of the Moon
364*/
365 static double ccampm[4][3] = {
366 { 1.097594e-1, /* ccampm[0][0] */
367 -2.223581e-2, /* ccampm[0][1] */
368 1.148966e-2 }, /* ccampm[0][2] */
369 { 2.896773e-7, /* ccampm[1][0] */
370 5.083103e-8, /* ccampm[1][1] */
371 5.658888e-8 }, /* ccampm[1][2] */
372 { 5.450474e-2, /* ccampm[2][0] */
373 1.002548e-2, /* ccampm[2][1] */
374 8.249439e-3 }, /* ccampm[2][2] */
375 { 1.438491e-7, /* ccampm[3][0] */
376 -2.291823e-8, /* ccampm[3][1] */
377 4.063015e-8 } /* ccampm[3][2] */
378 };
379
380/*
381**
382** Execution
383** ---------
384**
385** Control parameter ideq, and time arguments
386*/
387 ideq = ( deqx <= 0.0 ) ? 0 : 1;
388 dt = ( date - 15019.5 ) / 36525.0;
389
390/* Values of all elements for the instant date */
391 for ( k = 0; k < 8; k++ ) {
392 dlocal = dmod ( dcfel[0][k]
393 + dt * ( dcfel[1][k]
394 + dt * dcfel[2][k] ), D2PI );
395 if ( k == 0 ) {
396 dml = dlocal;
397 } else {
398 forbel[k-1] = dlocal;
399 }
400 }
401 deps = dmod ( dceps[0]
402 + dt * ( dceps[1]
403 + dt * dceps[2] ) , D2PI );
404 for ( k = 0; k < 17; k++ ) {
405 sorbel[k] = dmod ( ccsel[0][k]
406 + dt * ( ccsel[1][k]
407 + dt * ccsel[2][k] ), D2PI );
408 }
409
410/* Secular perturbations in longitude */
411 for ( k = 0; k < 4; k++ ) {
412 a = dmod ( ccsec[1][k] + dt * ccsec[2][k] , D2PI );
413 sn[k] = sin ( a );
414 }
415
416/* Periodic perturbations of the EMB (Earth-Moon barycentre) */
417 pertl = ccsec[0][0] * sn[0]
418 + ccsec[0][1] * sn[1]
419 + ( ccsec[0][2] + dt * ccsec3 ) * sn[2]
420 + ccsec[0][3] * sn[3];
421 pertld = 0.0;
422 pertr = 0.0;
423 pertrd = 0.0;
424 for ( k = 0; k < 15; k++ ) {
425 a = dmod ( dcargs[0][k] + dt * dcargs[1][k] , D2PI );
426 cosa = cos ( a );
427 sina = sin ( a );
428 pertl += ccamps[0][k] * cosa + ccamps[1][k] * sina;
429 pertr += ccamps[2][k] * cosa + ccamps[3][k] * sina;
430 if ( k < 10 ) {
431 pertld += ( ccamps[1][k] * cosa
432 - ccamps[0][k] * sina ) * ccamps[4][k];
433 pertrd += ( ccamps[3][k] * cosa
434 - ccamps[2][k] * sina ) * ccamps[4][k];
435 }
436 }
437
438/* Elliptic part of the motion of the EMB */
439 e = sorbel[0];
440 twoe = e + e;
441 esq = e * e;
442 dparam = 1.0 - esq;
443 g = forbel[0];
444 twog = g + g;
445 phi = twoe * ( ( 1.0 - esq / 8.0 ) * sin ( g )
446 + 5.0 * e * sin ( twog ) / 8.0
447 + 13.0 * esq * sin ( g + twog ) / 24.0 );
448 f = forbel[0] + phi;
449 sf = sin ( f );
450 cf = cos ( f );
451 dpsi = dparam / ( 1.0 + e * cf );
452 phid = twoe * ccsgd * ( ( 1.0 + esq * 1.5 ) * cf
453 + e * ( 1.25 - sf * sf / 2.0 ) );
454 psid = ccsgd * e * sf / sqrt ( dparam );
455
456/* Perturbed heliocentric motion of the EMB */
457 d1pdro = 1.0 + pertr;
458 drd = d1pdro * ( psid + dpsi * pertrd );
459 drld = d1pdro * dpsi * ( dcsld + phid + pertld );
460 dtl = dmod ( dml + phi + pertl , D2PI );
461 dsinls = sin ( dtl );
462 dcosls = cos ( dtl );
463 dxhd = drd * dcosls - drld * dsinls;
464 dyhd = drd * dsinls + drld * dcosls;
465
466/* Influence of eccentricity, evection and variation on the
467** geocentric motion of the Moon */
468 pertl = 0.0;
469 pertld = 0.0;
470 pertp = 0.0;
471 pertpd = 0.0;
472 for ( k = 0; k < 3; k++ ) {
473 a = dmod ( dcargm[0][k] + dt * dcargm[1][k] , D2PI );
474 sina = sin ( a );
475 cosa = cos ( a );
476 pertl += ccampm[0][k] * sina;
477 pertld += ccampm[1][k] * cosa;
478 pertp += ccampm[2][k] * cosa;
479 pertpd += - ccampm[3][k] * sina;
480 }
481
482/* Heliocentric motion of the Earth */
483 tl = forbel[1] + pertl;
484 sinlm = sin ( tl );
485 coslm = cos ( tl );
486 sigma = cckm / ( 1.0 + pertp );
487 a = sigma * ( ccmld + pertld );
488 b = sigma * pertpd;
489 dxhd += a * sinlm + b * coslm;
490 dyhd += - a * coslm + b * sinlm;
491 dzhd = - sigma * ccfdi * cos ( forbel[2] );
492
493/* Barycentric motion of the Earth */
494 dxbd = dxhd * dc1mme;
495 dybd = dyhd * dc1mme;
496 dzbd = dzhd * dc1mme;
497 for ( k = 0; k < 4; k++ ) {
498 plon = forbel[k+3];
499 pomg = sorbel[k+1];
500 pecc = sorbel[k+9];
501 tl = dmod( plon + 2.0 * pecc * sin ( plon - pomg ) , D2PI );
502 sinlp[k] = sin ( tl );
503 coslp[k] = cos ( tl );
504 dxbd += ccpamv[k] * ( sinlp[k] + pecc * sin ( pomg ) );
505 dybd += - ccpamv[k] * ( coslp[k] + pecc * cos ( pomg ) );
506 dzbd += - ccpamv[k] * sorbel[k+13] * cos ( plon - sorbel[k+5] );
507 }
508
509/* Transition to mean equator of date */
510 dcosep = cos ( deps );
511 dsinep = sin ( deps );
512 dyahd = dcosep * dyhd - dsinep * dzhd;
513 dzahd = dsinep * dyhd + dcosep * dzhd;
514 dyabd = dcosep * dybd - dsinep * dzbd;
515 dzabd = dsinep * dybd + dcosep * dzbd;
516
517/* Heliocentric coordinates of the Earth */
518 dr = dpsi * d1pdro;
519 flatm = ccim * sin ( forbel[2] );
520 a = sigma * cos ( flatm );
521 dxh = dr * dcosls - a * coslm;
522 dyh = dr * dsinls - a * sinlm;
523 dzh = - sigma * sin ( flatm );
524
525/* Barycentric coordinates of the Earth */
526 dxb = dxh * dc1mme;
527 dyb = dyh * dc1mme;
528 dzb = dzh * dc1mme;
529 for ( k = 0; k < 4; k++ ) {
530 flat = sorbel[k+13] * sin ( forbel[k+3] - sorbel[k+5] );
531 a = ccpam[k] * (1.0 - sorbel[k+9] * cos ( forbel[k+3] - sorbel[k+1]));
532 b = a * cos(flat);
533 dxb -= b * coslp[k];
534 dyb -= b * sinlp[k];
535 dzb -= a * sin ( flat );
536 }
537
538/* Transition to mean equator of date */
539 dyah = dcosep * dyh - dsinep * dzh;
540 dzah = dsinep * dyh + dcosep * dzh;
541 dyab = dcosep * dyb - dsinep * dzb;
542 dzab = dsinep * dyb + dcosep * dzb;
543
544/* Copy result components into vectors, correcting for FK4 equinox */
545 depj = slaEpj ( date );
546 deqcor = DS2R * ( 0.035 + ( 0.00085 * ( depj - b1950 ) ) );
547 dvh[0] = dxhd - deqcor * dyahd;
548 dvh[1] = dyahd + deqcor * dxhd;
549 dvh[2] = dzahd;
550 dvb[0] = dxbd - deqcor * dyabd;
551 dvb[1] = dyabd + deqcor * dxbd;
552 dvb[2] = dzabd;
553 dph[0] = dxh - deqcor * dyah;
554 dph[1] = dyah + deqcor * dxh;
555 dph[2] = dzah;
556 dpb[0] = dxb - deqcor * dyab;
557 dpb[1] = dyab + deqcor * dxb;
558 dpb[2] = dzab;
559
560/* Was precession to another equinox requested? */
561 if ( ideq != 0 ) {
562
563 /* Yes: compute precession matrix from MJD date to Julian Epoch deqx */
564 slaPrec ( depj, deqx, dprema );
565
566 /* Rotate dvh */
567 for ( j = 0; j < 3; j++ ) {
568 w = 0.0;
569 for ( i = 0; i < 3; i++ ) {
570 w += dprema[j][i] * dvh[i];
571 }
572 vw[j] = w;
573 }
574 for ( j = 0; j < 3; j++ ) {
575 dvh[j] = vw[j];
576 }
577
578 /* Rotate dvb */
579 for ( j = 0; j < 3; j++ ) {
580 w = 0.0;
581 for ( i = 0; i < 3; i++ ) {
582 w += dprema[j][i] * dvb[i];
583 }
584 vw[j] = w;
585 }
586 for ( j = 0; j < 3; j++ ) {
587 dvb[j] = vw[j];
588 }
589
590 /* Rotate dph */
591 for ( j = 0; j < 3; j++ ) {
592 w = 0.0;
593 for ( i = 0; i < 3; i++ ) {
594 w += dprema[j][i] * dph[i];
595 }
596 vw[j] = w;
597 }
598 for ( j = 0; j < 3; j++ ) {
599 dph[j] = vw[j];
600 }
601
602 /* Rotate dpb */
603 for ( j = 0; j < 3; j++ ) {
604 w = 0.0;
605 for ( i = 0; i < 3; i++ ) {
606 w += dprema[j][i] * dpb[i];
607 }
608 vw[j] = w;
609 }
610 for ( j = 0; j < 3; j++ ) {
611 dpb[j] = vw[j];
612 }
613 }
614}
Note: See TracBrowser for help on using the repository browser.