      SUBROUTINE nro_transf(tps,tu,xp,yp,dpsi0,deps0,itrans,xq)
C-------------------------------------------------------------------------C
C    Matrix transformation from the terrestrial reference system (TRS)    C
C    to the celestial reference system (CRS) at the date of observation   C
C                [CRS]= XQ [TRS]= xpn * xtheta * xyp [TRS]                C
C    Option (2), makes use of the nonrotating origin (Guinot 1979) to     C
C     realize the intermediate reference frame of date t,                 C
C      see IERS Conventions (1996), Chapter 5, p.34-39.                   C
C									  C
C    (a-h,o-z) REAL*8                                                     C
C                                                                         C
C      INPUT :                                                            C
C       tps= (JD - 2451545.0)/36525  (in Julian centuries)    and JD can  C
C           be considered as measured in TDB or TAI  (Julian days)        C
C       tu= Julian UT1 date - JD 2451545.0 UT1	(in days)		  C
C       xp, yp= "polar coordinates" of the CEP in the TRS  (in arcsecond) C
C                see IERS Conventions (1996), p.34                        C
C       dpsi0, deps0= celestial pole offsets (in arcsecond)               C
C                   see IERS Conventions (1996) p. 35                     C
C       itrans= 0  ==> nro using the complete equation of equinoxes       C
C             = 2  ==> degraded NRO, incomplete equation of equinoxes     C
C      (see "IERS Conventions (1996)", IERS Technical Note 21, pp. 34-35) C
C                                                                         C
C     OUTPUT :                                                            C
C       xq(3,3)= transformation matrix from the [TRS] to the [CRS]        C
C                dimension(3,3), (in radian)                              C
C                                                                         C
C     COMMON / conversion /                                               C
C       tpi= 2*pi                                                         C
C       secrad= arc second in radian                                      C
C       degrad= degrees in radian                                         C
C                                                                         C
C     SUBROUTINE INTERFACE :                                              C
C       called subroutines : mvtpole, angstellaire, precnut, mulmat       C
C-------------------------------------------------------------------------C	

       IMPLICIT NONE
       INTEGER itrans
       REAL*8 tps,tu,xp,yp,dpsi0,deps0,xq(3,3)
       REAL*8 tpi,secrad,degrad
       REAL*8 sp,ak,theta,x,y,stheta,E,d,s
       REAL*8 xpn(3,3),xyp(3,3),xtheta(3,3),x1(3,3) 

       COMMON/conversion/tpi,secrad,degrad


            ! Conversion factor
        tpi= 6.28318530717958648d0
        secrad= tpi / (360.d0*3600.d0)
        degrad= tpi / 360.d0

	CALL mvtpole(tps,xp,yp,sp,xyp)
	CALL angstellaire(tu,ak,theta,xtheta)
	CALL precnut(itrans,tps,dpsi0,deps0,x,y,stheta,E,d,s,xpn)

	CALL mulmat(xpn,xtheta,x1,3,3,3)
	CALL mulmat(x1,xyp,xq,3,3,3)
      END
C=======================================================================C


      SUBROUTINE mvtpole(tps,xp,yp,sp,xyp)
C-------------------------------------------------------------------------C
C      INPUT :                                                            C
C       tps= Julian date    (in Julian centuries)                         C
C       xp, yp= "polar coordinates" of the CEP in the TRS  (in arcsecond) C
C                see IERS Conventions (1996), p.34                        C
C                                                                         C
C     OUTPUT :                                                            C
C       sp=  (in radian)  see IERS Conventions (1996), p.34                 C
C       xyp(3,3)= R3(-sp)*R1(yp)*R2(xp) matrix   (in radian)              C
C                                                                         C
C     COMMON / conversion /                                               C
C       tpi= 2*pi                                                         C
C       secrad= arc second in radian                                      C
C       degrad= degrees in radian                                         C
C                                                                         C
C     SUBROUTINE INTERFACE :                                              C
C        caller subroutine : nro_transf                                   C
C        called subroutines: rota, mulmat                                 C
C-------------------------------------------------------------------------C

       IMPLICIT NONE
       REAL*8 tps,xp,yp,sp,xyp(3,3)
       REAL*8 tpi,secrad,degrad
       REAL*8 ac,am,ac2,am2
       REAL*8 xmat1(3,3),xmat2(3,3),xmat(3,3)

       COMMON/conversion/tpi,secrad,degrad


        ac=0.26d0        ! average amplitude of the Chandlerian wobble (in ")
        am=0.12d0        ! average amplitude of the annual wobble (in ")
        ac2=ac**2
        am2=am**2
        sp=.0015d0*(ac2/1.2d0+am2)*tps
        sp=sp*secrad

        CALL rota(1,yp*secrad,xmat1)
        CALL rota(2,xp*secrad,xmat2)
        CALL mulmat(xmat1,xmat2,xmat,3,3,3)
        CALL rota(3,-sp,xmat1)
        CALL mulmat(xmat1,xmat,xyp,3,3,3)

      END
C=======================================================================C


      SUBROUTINE angstellaire(tu,ak,theta,xtheta)
C-------------------------------------------------------------------------C
C      INPUT :                                                            C
C       tu= Julian UT1 date - JD 2451545.0 UT1  (in Julian days)          C
C                                                                         C
C     OUTPUT :                                                            C
C       ak=  constant                                                     C
C       theta= stellar angle   (in radian)                                C
C       xtheta(3,3)= R3(-theta) matrix   (in radian)                      C
C                                                                         C
C     COMMON / conversion /                                               C
C       tpi= 2*pi                                                         C
C       secrad= arc second in radian                                      C
C       degrad= degrees in radian                                         C
C                                                                         C
C     SUBROUTINE INTERFACE :                                              C
C        caller subroutine : nro_transf                                   C
C        called subroutines: rota                                         C
C-------------------------------------------------------------------------C

       IMPLICIT NONE
       REAL*8 tu,ak,theta,xtheta(3,3)
       REAL*8 tpi,secrad,degrad
       REAL*8 theta0 

       COMMON/conversion/tpi,secrad,degrad


        theta0=0.779057273264d0
        ak=1.00273781191135447d0

        theta=ak*tu+theta0
        theta=theta*tpi
        theta=DMOD(theta,tpi)

        CALL rota(3,-theta,xtheta)

      END
C=======================================================================C


      SUBROUTINE precnut(itrans,tps,dpsi0,deps0,xb,yb,stheta,E,d,s,xpn)
C-------------------------------------------------------------------------C
C      INPUT :                                                            C
C       itrans= 0  ==> nro using the complete equation of equinoxes       C
C             = 2  ==> degraded NRO, incomplete equation of equinoxes     C
C                see IERS Conventions (1996), pp.34-35                    C
C       tps= Julian date    (in Julian centuries)                         C
C       dpsi0, deps0= celestial pole offsets (in arcsecond)               C
C                   see IERS Conventions (1996), p. 35                    C
C                                                                         C
C     OUTPUT :                                                            C
C       xb, yb, stheta=           (in arcsecond)                          C 
C       E, d, s=  angles     (in radian)                                  C
C       xpn(3,3)= R3(-E)*R2(-d)*R3(E)*R3(s) matrix   (in radian)          C
C                                                                         C
C     COMMON / conversion /                                               C
C       tpi= 2*pi                                                         C
C       secrad= arc second in radian                                      C
C       degrad= degrees in radian                                         C
C                                                                         C
C     SUBROUTINE INTERFACE :                                              C
C        caller subroutine : nro_transf                                   C
C        called subroutines: coord_XY, rota, mulmat, arctangente          C
C-------------------------------------------------------------------------C

       IMPLICIT NONE
       INTEGER itrans
       REAL*8 tps,dpsi0,deps0,xb,yb,stheta,E,d,s,xpn(3,3)
       REAL*8 tpi,secrad,degrad
       REAL*8 arctangente
       REAL*8 sineps0,xa,ya,zb,v,cosE,sinE
       REAL*8 xmat1(3,3),xmat2(3,3),xmat(3,3),smat(3,3)

       COMMON/conversion/tpi,secrad,degrad


        CALL coord_XY(itrans,tps,xb,yb,stheta)
        sineps0=DSIN(84381.448D0*secrad)
        xb=xb+dpsi0*sineps0
        yb=yb+deps0
        xa=xb*secrad
        ya=yb*secrad
        s=stheta-(xb*ya)/2.d0
        s=s*secrad

        zb=DSQRT(1.d0-(xa**2+ya**2))
        v=-DSQRT(xb**2+yb**2)
        cosE=xb/v
        sinE=yb/v
        v=v*secrad
        E=arctangente(sinE,cosE)
        d=arctangente(v,zb)
        E=E*degrad
        d=d*degrad

        CALL rota(3,E,xmat)
        CALL rota(2,-d,xmat1)
        CALL mulmat(xmat1,xmat,xmat2,3,3,3)
        CALL rota(3,-E,xmat1)
        CALL mulmat(xmat1,xmat2,xmat,3,3,3)

        CALL rota(3,s,smat)
        CALL mulmat(xmat,smat,xpn,3,3,3)

      END
C=====================================================================C


	SUBROUTINE coord_XY(itrans,tps,xb,yb,sb)
C-------------------------------------------------------------------------C
C      INPUT :                                                            C
C       itrans= 0  ==> nro using the complete equation of equinoxes       C
C             = 2  ==> degraded NRO, incomplete equation of equinoxes     C
C                see IERS Conventions (1996), p.35-36                     C
C       tps= Julian date    (in Julian centuries)                         C
C                                                                         C
C     OUTPUT :                                                            C
C       xb, yb, sb=           (in arcsecond)                              C
C                                                                         C
C     COMMON / conversion /                                               C
C       tpi= 2*pi                                                         C
C       secrad= arc second in radian                                      C
C       degrad= degrees in radian                                         C
C                                                                         C
C     SUBROUTINE INTERFACE :                                              C
C        caller subroutine : precnut                                      C
C-------------------------------------------------------------------------C

       IMPLICIT NONE
       INTEGER itrans
       REAL*8 tps,xb,yb,sb
       REAL*8 tpi,secrad,degrad
       INTEGER matrix(7,106),mtrxsin(7,20),mtrxcos(7,16)
       INTEGER mtrx1(7,20),mtrx2(7,20),mtrx3(7,20)
       INTEGER mtrx4(7,20),mtrx5(7,20),mtrx6(7,6)
       INTEGER i,j
       REAL*8 tps2,tps3,r,dangle(5)
       REAL*8 deltapsi,deltaeps,angle,sineps0

       COMMON/conversion/tpi,secrad,degrad

       EQUIVALENCE (matrix(1,  1),mtrx1 (1, 1))
       EQUIVALENCE (matrix(1, 21),mtrx2 (1, 1))
       EQUIVALENCE (matrix(1, 41),mtrx3 (1, 1))
       EQUIVALENCE (matrix(1, 61),mtrx4 (1, 1))
       EQUIVALENCE (matrix(1, 81),mtrx5 (1, 1))
       EQUIVALENCE (matrix(1,101),mtrx6 (1, 1))

       	DATA mtrx1 /
     +     0, 0, 0, 0, 1,  -171996,      92025,
     +     0, 0, 0, 0, 2,     2062,       -895,
     +    -2, 0, 2, 0, 1,       46,        -24,
     +     2, 0,-2, 0, 0,       11,          0,
     +    -2, 0, 2, 0, 2,       -3,          1,
     +     1,-1, 0,-1, 0,       -3,          0,
     +     0,-2, 2,-2, 1,       -2,          1,
     +     2, 0,-2, 0, 1,        1,          0,
     +     0, 0, 2,-2, 2,   -13187,       5736,
     +     0, 1, 0, 0, 0,     1426,         54,
     +     0, 1, 2,-2, 2,     -517,        224,
     +     0,-1, 2,-2, 2,      217,        -95,
     +     0, 0, 2,-2, 1,      129,        -70,
     +     2, 0, 0,-2, 0,       48,          1,
     +     0, 0, 2,-2, 0,      -22,          0,
     +     0, 2, 0, 0, 0,       17,          0,
     +     0, 1, 0, 0, 1,      -15,          9,
     +     0, 2, 2,-2, 2,      -16,          7,
     +     0,-1, 0, 0, 1,      -12,          6,
     +    -2, 0, 0, 2, 1,       -6,          3/
       	DATA mtrx2 /
     +     0,-1, 2,-2, 1,       -5,          3,
     +     2, 0, 0,-2, 1,        4,         -2,
     +     0, 1, 2,-2, 1,        4,         -2,
     +     1, 0, 0,-1, 0,       -4,          0,
     +     2, 1, 0,-2, 0,        1,          0,
     +     0, 0,-2, 2, 1,        1,          0,
     +     0, 1,-2, 2, 0,       -1,          0,
     +     0, 1, 0, 0, 2,        1,          0,
     +    -1, 0, 0, 1, 1,        1,          0,
     +     0, 1, 2,-2, 0,       -1,          0,
     +     0, 0, 2, 0, 2,    -2274,        977,
     +     1, 0, 0, 0, 0,      712,         -7,
     +     0, 0, 2, 0, 1,     -386,        200,
     +     1, 0, 2, 0, 2,     -301,        129,
     +     1, 0, 0,-2, 0,     -158,         -1,
     +    -1, 0, 2, 0, 2,      123,        -53,
     +     0, 0, 0, 2, 0,       63,         -2,
     +     1, 0, 0, 0, 1,       63,        -33,
     +    -1, 0, 0, 0, 1,      -58,         32,
     +    -1, 0, 2, 2, 2,      -59,         26/
       	DATA mtrx3 /
     +     1, 0, 2, 0, 1,      -51,         27,
     +     0, 0, 2, 2, 2,      -38,         16,
     +     2, 0, 0, 0, 0,       29,         -1,
     +     1, 0, 2,-2, 2,       29,        -12,
     +     2, 0, 2, 0, 2,      -31,         13,
     +     0, 0, 2, 0, 0,       26,         -1,
     +    -1, 0, 2, 0, 1,       21,        -10,
     +    -1, 0, 0, 2, 1,       16,         -8,
     +     1, 0, 0,-2, 1,      -13,          7,
     +    -1, 0, 2, 2, 1,      -10,          5,
     +     1, 1, 0,-2, 0,       -7,          0,
     +     0, 1, 2, 0, 2,        7,         -3,
     +     0,-1, 2, 0, 2,       -7,          3,
     +     1, 0, 2, 2, 2,       -8,          3,
     +     1, 0, 0, 2, 0,        6,          0,
     +     2, 0, 2,-2, 2,        6,         -3,
     +     0, 0, 0, 2, 1,       -6,          3,
     +     0, 0, 2, 2, 1,       -7,          3,
     +     1, 0, 2,-2, 1,        6,         -3,
     +     0, 0, 0,-2, 1,       -5,          3/
       	DATA mtrx4 /
     +     1,-1, 0, 0, 0,        5,          0,
     +     2, 0, 2, 0, 1,       -5,          3,
     +     0, 1, 0,-2, 0,       -4,          0,
     +     1, 0,-2, 0, 0,        4,          0,
     +     0, 0, 0, 1, 0,       -4,          0,
     +     1, 1, 0, 0, 0,       -3,          0,
     +     1, 0, 2, 0, 0,        3,          0,
     +     1,-1, 2, 0, 2,       -3,          1,
     +    -1,-1, 2, 2, 2,       -3,          1,
     +    -2, 0, 0, 0, 1,       -2,          1,
     +     3, 0, 2, 0, 2,       -3,          1,
     +     0,-1 ,2, 2, 2,       -3,          1,
     +     1, 1, 2, 0, 2,        2,         -1,
     +    -1, 0, 2,-2, 1,       -2,          1,
     +     2, 0, 0, 0, 1,        2,         -1,
     +     1, 0, 0, 0, 2,       -2,          1,
     +     3, 0, 0, 0, 0,        2,          0,
     +     0, 0, 2, 1, 2,        2,         -1,
     +    -1, 0, 0, 0, 2,        1,         -1,
     +     1, 0, 0,-4, 0,       -1,          0/
       	DATA mtrx5 /
     +    -2, 0, 2, 2, 2,        1,         -1,
     +    -1, 0, 2, 4, 2,       -2,          1,
     +     2, 0, 0,-4, 0,       -1,          0,
     +     1, 1, 2,-2, 2,        1,         -1,
     +     1, 0, 2, 2, 1,       -1,          1,
     +    -2, 0, 2, 4, 2,       -1,          1,
     +    -1, 0, 4, 0, 2,        1,          0,
     +     1,-1, 0,-2, 0,        1,          0,
     +     2, 0, 2,-2, 1,        1,         -1,
     +     2, 0, 2, 2, 2,       -1,          0,
     +     1, 0, 0, 2, 1,       -1,          0,
     +     0, 0, 4,-2, 2,        1,          0,
     +     3, 0, 2,-2, 2,        1,          0,
     +     1, 0, 2,-2, 0,       -1,          0,
     +     0, 1, 2, 0, 1,        1,          0,
     +    -1,-1, 0, 2, 1,        1,          0,
     +     0, 0,-2, 0, 1,       -1,          0,
     +     0, 0, 2,-1, 2,       -1,          0,
     +     0, 1, 0, 2, 0,       -1,          0,
     +     1, 0,-2,-2, 0,       -1,          0/
       	DATA mtrx6 /
     +     0,-1, 2, 0, 1,       -1,          0,
     +     1, 1, 0,-2, 1,       -1,          0,
     +     1, 0,-2, 2, 0,       -1,          0,
     +     2, 0, 0, 2, 0,        1,          0,
     +     0, 0, 2, 4, 2,       -1,          0,
     +     0, 1, 0, 1, 0,        1,          0/
	DATA mtrxsin /
     +     0, 0, 0, 0, 1,      -335,     15299,
     +     0, 0, 0, 0, 2,        -4,      -183,
     +    -2, 0, 2, 0, 1,         0,        -4,
     +     0, 0, 2,-2, 2,        21,      1173,
     +     0, 1, 0, 0, 0,       -17,      -127,
     +     0, 1, 2,-2, 2,         6,        46,
     +     0,-1, 2,-2, 2,         0,       -19,
     +     0, 0, 2,-2, 1,         0,       -12,
     +     2, 0, 0,-2, 0,         0,        -4,
     +     0, 0, 2, 0, 2,         4,       202,
     +     1, 0, 0, 0, 0,         0,       -63,
     +     0, 0, 2, 0, 1,         0,        34,
     +     1, 0, 2, 0, 2,         0,        27,
     +     1, 0, 0,-2, 0,         0,        14,
     +    -1, 0, 2, 0, 2,         0,       -11,
     +     0, 0, 0, 2, 0,         0,        -6,
     +     1, 0, 0, 0, 1,         0,        -6,
     +    -1, 0, 0, 0, 1,         0,         5,
     +    -1, 0, 2, 2, 2,         0,         5,
     +     1, 0, 2, 0, 1,         0,         5/
	DATA mtrxcos /
     +     0, 0, 0, 0, 1,     20578,        89,
     +     0, 0, 0, 0, 2,      -201,         5,
     +    -2, 0, 2, 0, 1,        -5,         0,
     +     0, 0, 2,-2, 2,      1283,       -31,
     +     0, 1, 0, 0, 0,        12,         0,
     +     0, 1, 2,-2, 2,        50,        -6,
     +     0,-1, 2,-2, 2,       -21,         0,
     +     0, 0, 2,-2, 1,       -16,         0,
     +     0, 0, 2, 0, 2,       218,        -5,
     +     0, 0, 2, 0, 1,        45,         0,
     +     1, 0, 2, 0, 2,        29,         0,
     +    -1, 0, 2, 0, 2,       -12,         0,
     +     1, 0, 0, 0, 1,        -7,         0,
     +    -1, 0, 0, 0, 1,         7,         0,
     +    -1, 0, 2, 2, 2,         6,         0,
     +     1, 0, 2, 0, 1,         6,         0/


        tps2=tps**2
        tps3=tps2 * tps

        xb=2004.3109d0*tps-.42665d0*tps2-.198656d0*tps3+.000014d0*tps**4
        yb=-.00013d0-22.40992d0*tps2+.001836d0*tps3+.001113d0*tps**4
        sb=.00385d0*tps-.07259d0*tps3

        r= 1296000.d0
        dangle(1)= 485866.733d0 + (1325.d0*r + 715922.633d0)*tps +
     .            31.310d0*tps2 + 0.064d0*tps3
        dangle(2)= 1287099.804d0 + (99.d0*r + 1292581.224d0)*tps -
     .            0.577d0*tps2 - 0.012d0*tps3
        dangle(3)= 335778.877d0 + (1342.d0*r + 295263.137d0)*tps -
     .            13.257d0*tps2 + 0.011d0*tps3
        dangle(4)= 1072261.307d0 + (1236.d0*r + 1105601.328d0)*tps -
     .            6.891d0*tps2 + 0.019d0*tps3
        dangle(5)= 450160.280d0 - (5.d0*r + 482890.539d0)*tps +
     .            7.455d0*tps2 + 0.008d0*tps3

        DO i=1,5
          dangle(i)= dangle(i) * secrad
        ENDDO
        deltapsi=0.d0
        deltaeps=0.d0

        DO i=1,106
          angle=0.d0
          DO j=1,5
            angle=angle+ DFLOAT(matrix(j,i))*dangle(j)
            angle=dmod(angle,tpi)
          ENDDO
          deltapsi=deltapsi+ DFLOAT(matrix(6,i))*dsin(angle) / 1.d4
          deltaeps=deltaeps+ DFLOAT(matrix(7,i))*dcos(angle) / 1.d4
        ENDDO

C	  sineps0=0.39777716d0
        sineps0=DSIN(84381.448d0*secrad)
        deltapsi=deltapsi * sineps0

        xb=xb+deltapsi-.00035d0*dsin(2.d0*dangle(5))
        angle=2.d0*(dangle(3)-dangle(4))+3.d0*dangle(5)
        xb=xb-.000049d0*dsin(angle)
        yb=yb+deltaeps+.00013d0*dcos(2.d0*dangle(5))
        yb=yb+.00002d0*dcos(angle)
        yb=yb-.00002d0*dcos(2.d0*(dangle(3)-dangle(4))+dangle(5))

          ! NRO, complete equation of equinoxes (after 26 february 1997)
        IF (itrans.EQ.0) THEN 
          sb=sb-.00265d0*(dsin(dangle(5)))
          sb=sb-.00006d0*(dsin(2.d0*dangle(5)))
        ENDIF

          ! degraded NRO, incomplete equation of equinoxes
          !  = classical transformation (until 26 february 1997)
        IF (itrans.EQ.2) THEN 
          sb=sb+.000153D0*(dsin(2.D0*dangle(5)))
        ENDIF

        DO i=1,20
          angle=0.d0
          DO j=1,5
            angle=angle + DFLOAT(mtrxsin(j,i)) * dangle(j)
            angle=dmod(angle,tpi)
          ENDDO
          xb=xb + tps * dsin(angle) * DFLOAT(mtrxsin(6,i)) / 1.d5
          yb=yb + tps * dsin(angle) * DFLOAT(mtrxsin(7,i)) / 1.d5
        ENDDO
        DO i=1,16
          angle=0.d0
          DO j=1,5
            angle=angle + DFLOAT(mtrxcos(j,i)) * dangle(j)
            angle=dmod(angle,tpi)
          ENDDO
          xb=xb + tps * dcos(angle) * DFLOAT(mtrxcos(6,i)) / 1.d5
          yb=yb + tps * dcos(angle) * DFLOAT(mtrxcos(7,i)) / 1.d5
        ENDDO
        angle=2.d0 * (dangle(5)+dangle(3)-dangle(4))
        angle=dmod(angle,tpi)
        xb=xb+(.00204d0*dsin(dangle(5))+.00016d0*dsin(angle))*tps2
        xb=xb+(.00006d0*dcos(dangle(5)))*tps2
        yb=yb-.00231d0*dcos(dangle(5))*tps2
        yb=yb -.00014d0 *dcos(angle) * tps2
        yb=yb-.000034d0*dsin(dangle(5))*tps2
        sb=sb+(.00074d0*dsin(dangle(5))+.00006d0*dsin(angle))*tps2

C               geodesic nutation     Not yet introduce in the usual procedure
c         geonut=0.000153d0*DSIN(dangle(2))+0.0000002d0*DSIN(2.d0*dangle(2))
c         xb=xb-geonut*sineps0

      END
C=======================================================================C


      SUBROUTINE rota(iaxe,angl,rotmat)
C-----------------------------------------------------------------------C
C   rotation matrix about  iaxe-axis by an angle = angl 		C
C	iaxe=1 ==> axis X,  iaxe=2 ==> axis Y,  iaxe=3 ==> axis Z	C
C                                                                       C
C   (a-h,o-z) REAL*8                                                    C
C                                                                       C
C     SUBROUTINE INTERFACE :                                            C
C        caller subroutine : mvtpole, angstellaire, precnut             C
C-----------------------------------------------------------------------C

       IMPLICIT NONE
       INTEGER iaxe
       REAL*8 angl,rotmat(3,3)
       REAL*8 cangl,sangl


        IF ((iaxe.GE.1).AND.(iaxe.LE.3)) THEN
          cangl=dcos(angl)
          sangl=dsin(angl)
          IF ((ABS(cangl)).LE.1.d-15) THEN
            cangl=0.d0
          ELSE
            IF ((ABS(sangl)).LE.1.d-15) THEN
              sangl=0.d0
            ENDIF
          ENDIF

          IF (iaxe.EQ.1) THEN
            rotmat(1,1)=1.d0
            rotmat(1,2)=0.d0
            rotmat(1,3)=0.d0

            rotmat(2,1)=0.d0
            rotmat(2,2)=cangl
            rotmat(2,3)=sangl

            rotmat(3,1)=0.d0
            rotmat(3,2)=-sangl
            rotmat(3,3)=cangl
          ELSE
            IF (iaxe.EQ.2) THEN
              rotmat(1,1)=cangl
              rotmat(1,2)=0.d0
              rotmat(1,3)=-sangl

              rotmat(2,1)=0.d0
              rotmat(2,2)=1.d0
              rotmat(2,3)=0.d0

              rotmat(3,1)=sangl
              rotmat(3,2)=0.d0
              rotmat(3,3)=cangl
            ELSE
              rotmat(1,1)=cangl
              rotmat(1,2)=sangl
              rotmat(1,3)=0.d0

              rotmat(2,1)=-sangl
              rotmat(2,2)=cangl
              rotmat(2,3)=0.d0

              rotmat(3,1)=0.d0
              rotmat(3,2)=0.d0
              rotmat(3,3)=1.d0
            ENDIF
          ENDIF
        ELSE
          PRINT*,'the rotation about the axis ',iaxe,' is not defined'
        ENDIF

      END
C=======================================================================C


      SUBROUTINE mulmat(x1,x2,xr,ni,nc,nj)
C-----------------------------------------------------------------------C
C  Multiplication of matrices                                           C
C            xr(ni,nj) = x1(ni,nc) * x2(nc,nj)                          C
C                                                                       C
C     (a-h,o-z) REAL*8                                                  C
C                                                                       C
C     SUBROUTINE INTERFACE :                                            C
C        caller subroutine : nro_transf, mvtpole, precnut               C
C-----------------------------------------------------------------------C

       IMPLICIT NONE
       INTEGER ni,nc,nj
       REAL*8 x1(ni,nc),x2(nc,nj),xr(ni,nj)
       INTEGER i,j,k
       REAL*8 cal


        DO i=1,ni
          DO j=1,nj
            cal=0.d0
            DO k=1,nc
              cal= cal + x1(i,k) * x2(k,j)
            ENDDO
            xr(i,j)=cal
          ENDDO
        ENDDO 

      END
C=======================================================================C


      FUNCTION arctangente(a1,a2)
C-----------------------------------------------------------------------C
C   computes the angle in degree with the value of the sine,a1, and the C
C     cosine, a2.                                                       C
C                                                                       C
C     (a-h,o-z) REAL*8                                                  C
C                                                                       C
C     SUBROUTINE INTERFACE :                                            C
C        caller subroutine : precnut                                    C
C-----------------------------------------------------------------------C

        IMPLICIT NONE
        REAL*8 arctangente,a1,a2
        REAL*8 pi,radeg,angle

         pi=4.d0*DATAN(1.d0)
         radeg=180.d0/pi
         IF (ABS(a2).LT.1.d-17) THEN
           angle=DSIGN(pi/2.d0,a1)
         ELSE
           angle=DATAN(a1/a2)
           IF (a2.LT.0.d0) THEN
             angle=DSIGN(pi-ABS(angle),a1)
           ENDIF
         ENDIF
         IF (ABS(angle).LT.1.d-15) THEN
           angle=0.d0
         ENDIF

         arctangente=angle*radeg

       END
C=======================================================================C

