source: trunk/MagicSoft/slalib/prec.c

Last change on this file was 731, checked in by tbretz, 24 years ago
*** empty log message ***
  • Property svn:executable set to *
File size: 2.3 KB
Line 
1#include "slalib.h"
2#include "slamac.h"
3void slaPrec ( double ep0, double ep1, double rmatp[3][3] )
4/*
5** - - - - - - - -
6** s l a P r e c
7** - - - - - - - -
8**
9** Form the matrix of precession between two epochs (IAU 1976, FK5).
10**
11** (double precision)
12**
13** Given:
14** ep0 double beginning epoch
15** ep1 double ending epoch
16**
17** Returned:
18** rmatp double[3][3] precession matrix
19**
20** Notes:
21**
22** 1) The epochs are TDB (loosely ET) Julian epochs.
23**
24** 2) The matrix is in the sense v(ep1) = rmatp * v(ep0) .
25**
26** 3) Though the matrix method itself is rigorous, the precession
27** angles are expressed through canonical polynomials which are
28** valid only for a limited time span. There are also known
29** errors in the IAU precession rate. The absolute accuracy
30** of the present formulation is better than 0.1 arcsec from
31** 1960AD to 2040AD, better than 1 arcsec from 1640AD to 2360AD,
32** and remains below 3 arcsec for the whole of the period
33** 500BC to 3000AD. The errors exceed 10 arcsec outside the
34** range 1200BC to 3900AD, exceed 100 arcsec outside 4200BC to
35** 5600AD and exceed 1000 arcsec outside 6800BC to 8200AD.
36** The SLALIB routine slaPrecl implements a more elaborate
37** model which is suitable for problems spanning several
38** thousand years.
39**
40** References:
41** Lieske,J.H., 1979. Astron. Astrophys.,73,282.
42** equations (6) & (7), p283.
43** Kaplan,G.H., 1981. USNO circular no. 163, pa2.
44**
45** Called: slaDeuler
46**
47** Defined in slamac.h: DAS2R
48**
49** Last revision: 10 July 1994
50**
51** Copyright P.T.Wallace. All rights reserved.
52*/
53{
54 double t0, t, tas2r, w, zeta, z, theta;
55
56/* Interval between basic epoch J2000.0 and beginning epoch (JC) */
57 t0 = ( ep0 - 2000.0 ) / 100.0;
58
59/* Interval over which precession required (JC) */
60 t = ( ep1 - ep0 ) / 100.0;
61
62/* Euler angles */
63 tas2r = t * DAS2R;
64 w = 2306.2181 + ( ( 1.39656 - ( 0.000139 * t0 ) ) * t0 );
65 zeta = (w + ( ( 0.30188 - 0.000344 * t0 ) + 0.017998 * t ) * t ) * tas2r;
66 z = (w + ( ( 1.09468 + 0.000066 * t0 ) + 0.018203 * t ) * t ) * tas2r;
67 theta = ( ( 2004.3109 + ( - 0.85330 - 0.000217 * t0 ) * t0 )
68 + ( ( -0.42665 - 0.000217 * t0 ) - 0.041833 * t ) * t ) * tas2r;
69
70/* Rotation matrix */
71 slaDeuler ( "ZYZ", -zeta, theta, -z, rmatp );
72}
Note: See TracBrowser for help on using the repository browser.