| 1 | #include "slalib.h" | 
|---|
| 2 | #include "slamac.h" | 
|---|
| 3 | void slaPrebn ( double bep0, double bep1, double rmatp[3][3] ) | 
|---|
| 4 | /* | 
|---|
| 5 | **  - - - - - - - - - | 
|---|
| 6 | **   s l a P r e b n | 
|---|
| 7 | **  - - - - - - - - - | 
|---|
| 8 | ** | 
|---|
| 9 | **  Generate the matrix of precession between two epochs, | 
|---|
| 10 | **  using the old, pre-IAU1976, Bessel-Newcomb model, using | 
|---|
| 11 | **  Kinoshita's formulation (double precision) | 
|---|
| 12 | ** | 
|---|
| 13 | **  Given: | 
|---|
| 14 | **     BEP0    double        beginning Besselian epoch | 
|---|
| 15 | **     BEP1    double        ending Besselian epoch | 
|---|
| 16 | ** | 
|---|
| 17 | **  Returned: | 
|---|
| 18 | **     RMATP   double[3][3]  precession matrix | 
|---|
| 19 | ** | 
|---|
| 20 | **  The matrix is in the sense   v(bep1)  =  rmatp * v(bep0) | 
|---|
| 21 | ** | 
|---|
| 22 | **  Reference: | 
|---|
| 23 | **     Kinoshita, H. (1975) 'Formulas for precession', SAO Special | 
|---|
| 24 | **     Report No. 364, Smithsonian Institution Astrophysical | 
|---|
| 25 | **     Observatory, Cambridge, Massachusetts. | 
|---|
| 26 | ** | 
|---|
| 27 | **  Called:  slaDeuler | 
|---|
| 28 | ** | 
|---|
| 29 | **  Defined in slamac.h:  DAS2R | 
|---|
| 30 | ** | 
|---|
| 31 | **  Last revision:   30 October 1993 | 
|---|
| 32 | ** | 
|---|
| 33 | **  Copyright P.T.Wallace.  All rights reserved. | 
|---|
| 34 | */ | 
|---|
| 35 | { | 
|---|
| 36 | double bigt, t, tas2r, w, zeta, z, theta; | 
|---|
| 37 |  | 
|---|
| 38 | /* Interval between basic epoch B1850.0 and beginning epoch in TC */ | 
|---|
| 39 | bigt  = ( bep0 - 1850.0 ) / 100.0; | 
|---|
| 40 |  | 
|---|
| 41 | /* Interval over which precession required, in tropical centuries */ | 
|---|
| 42 | t = ( bep1 - bep0 ) / 100.0; | 
|---|
| 43 |  | 
|---|
| 44 | /* Euler angles */ | 
|---|
| 45 | tas2r = t * DAS2R; | 
|---|
| 46 | w = 2303.5548 + ( 1.39720 + 0.000059 * bigt ) * bigt; | 
|---|
| 47 | zeta = (w + ( 0.30242 - 0.000269 * bigt + 0.017996 * t ) * t ) * tas2r; | 
|---|
| 48 | z = (w + ( 1.09478 + 0.000387 * bigt + 0.018324 * t ) * t ) * tas2r; | 
|---|
| 49 | theta = ( 2005.1125 + ( - 0.85294 - 0.000365* bigt ) * bigt + | 
|---|
| 50 | ( - 0.42647 - 0.000365 * bigt - 0.041802 * t ) * t ) * tas2r; | 
|---|
| 51 |  | 
|---|
| 52 | /* Rotation matrix */ | 
|---|
| 53 | slaDeuler ( "ZYZ", -zeta, theta, -z, rmatp ); | 
|---|
| 54 | } | 
|---|