Initial commit

1 files changed, 62 insertions, 0 deletions

diff --git a/src/slalib/prebn.f b/src/slalib/prebn.f
new file mode 100644
index 0000000..a99dd51
--- /dev/null
+++ b/src/slalib/prebn.f
@@ -0,0 +1,62 @@
+      SUBROUTINE sla_PREBN (BEP0, BEP1, RMATP)
+*+
+*     - - - - - -
+*      P R E B N
+*     - - - - - -
+*
+*  Generate the matrix of precession between two epochs,
+*  using the old, pre-IAU1976, Bessel-Newcomb model, using
+*  Kinoshita's formulation (double precision)
+*
+*  Given:
+*     BEP0    dp         beginning Besselian epoch
+*     BEP1    dp         ending Besselian epoch
+*
+*  Returned:
+*     RMATP  dp(3,3)    precession matrix
+*
+*  The matrix is in the sense   V(BEP1)  =  RMATP * V(BEP0)
+*
+*  Reference:
+*     Kinoshita, H. (1975) 'Formulas for precession', SAO Special
+*     Report No. 364, Smithsonian Institution Astrophysical
+*     Observatory, Cambridge, Massachusetts.
+*
+*  Called:  sla_DEULER
+*
+*  P.T.Wallace   Starlink   23 August 1996
+*
+*  Copyright (C) 1996 Rutherford Appleton Laboratory
+*-
+
+      IMPLICIT NONE
+
+      DOUBLE PRECISION BEP0,BEP1,RMATP(3,3)
+
+*  Arc seconds to radians
+      DOUBLE PRECISION AS2R
+      PARAMETER (AS2R=0.484813681109535994D-5)
+
+      DOUBLE PRECISION BIGT,T,TAS2R,W,ZETA,Z,THETA
+
+
+
+*  Interval between basic epoch B1850.0 and beginning epoch in TC
+      BIGT = (BEP0-1850D0)/100D0
+
+*  Interval over which precession required, in tropical centuries
+      T = (BEP1-BEP0)/100D0
+
+*  Euler angles
+      TAS2R = T*AS2R
+      W = 2303.5548D0+(1.39720D0+0.000059D0*BIGT)*BIGT
+
+      ZETA = (W+(0.30242D0-0.000269D0*BIGT+0.017996D0*T)*T)*TAS2R
+      Z = (W+(1.09478D0+0.000387D0*BIGT+0.018324D0*T)*T)*TAS2R
+      THETA = (2005.1125D0+(-0.85294D0-0.000365D0*BIGT)*BIGT+
+     :        (-0.42647D0-0.000365D0*BIGT-0.041802D0*T)*T)*TAS2R
+
+*  Rotation matrix
+      CALL sla_DEULER('ZYZ',-ZETA,THETA,-Z,RMATP)
+
+      END