C
C David J. Heisterberg
C The Ohio Supercomputer Center
C 1224 Kinnear Rd.
C Columbus, OH  43212-1163
C (614)292-6036
C djh@ccl.net    djh@ohstpy.bitnet    ohstpy::djh
C
C ANALYZE
C analyze the x to y fit and the fitting matrix
C
      SUBROUTINE analyz (n, x, y, w, u)
      IMPLICIT CHARACTER (A-Z)
C
      INTEGER n
      DOUBLEPRECISION x (3, n)
      DOUBLEPRECISION y (3, n)
      DOUBLEPRECISION w (n)
      DOUBLEPRECISION u (3, 3)
C
      DOUBLEPRECISION err
      DOUBLEPRECISION wnorm
      DOUBLEPRECISION urnorm (3)
      DOUBLEPRECISION ucnorm (3)
      INTEGER i, j
C
C find the mean sum of squares error
C
      err = 0.0D0
      wnorm = 0.0D0
      DO 10000 i = 1, n
        err = err + w (i) * ((y (1, i) - x (1, i)) ** 2 +
     &                       (y (2, i) - x (2, i)) ** 2 +
     &                       (y (3, i) - x (3, i)) ** 2)
        wnorm = wnorm + w (i)
10000   CONTINUE
      err = SQRT (err / wnorm)
C
C find the row and column norms of u
C
      ucnorm (1) = u (1, 1) ** 2 + u (2, 1) ** 2 + u (3, 1) ** 2
      ucnorm (2) = u (1, 2) ** 2 + u (2, 2) ** 2 + u (3, 2) ** 2
      ucnorm (3) = u (1, 3) ** 2 + u (2, 3) ** 2 + u (3, 3) ** 2
      urnorm (1) = u (1, 1) ** 2 + u (1, 2) ** 2 + u (1, 3) ** 2
      urnorm (2) = u (2, 1) ** 2 + u (2, 2) ** 2 + u (2, 3) ** 2
      urnorm (3) = u (3, 1) ** 2 + u (3, 2) ** 2 + u (3, 3) ** 2
C
C write the error and u norms
C
      WRITE (*, *)
      DO 11000 i = 1, 3
        WRITE (*, 1000) (u (i, j), j = 1, 3), urnorm (i)
11000   CONTINUE
      WRITE (*, *)
      WRITE (*, 1000) (ucnorm (j), j = 1, 3), err
C
      RETURN
C
 1000 FORMAT (1X, 3E18.10, 4X, E18.10)
C
      END
C
C CENTER
C center a molecule about its weighted centroid or other origin
C
      SUBROUTINE center (n, x, w, io, o, t)
      IMPLICIT CHARACTER (A-Z)
C
      INTEGER n
      DOUBLEPRECISION x (3, n)
      DOUBLEPRECISION w (n)
      LOGICAL io
      DOUBLEPRECISION o (3)
      DOUBLEPRECISION t (3)
C
      DOUBLEPRECISION wnorm
      INTEGER i
C
      IF (io) THEN
        t (1) = o (1)
        t (2) = o (2)
        t (3) = o (3)
      ELSE
        t (1) = 0.0D0
        t (2) = 0.0D0
        t (3) = 0.0D0
        wnorm = 0.0D0
        DO 10000 i = 1, n
          t (1) = t (1) + x (1, i) * SQRT (w (i))
          t (2) = t (2) + x (2, i) * SQRT (w (i))
          t (3) = t (3) + x (3, i) * SQRT (w (i))
          wnorm = wnorm + SQRT (w (i))
10000     CONTINUE
        t (1) = t (1) / wnorm
        t (2) = t (2) / wnorm
        t (3) = t (3) / wnorm
      ENDIF
C
      DO 10100 i = 1, n
        x (1, i) = x (1, i) - t (1)
        x (2, i) = x (2, i) - t (2)
        x (3, i) = x (3, i) - t (3)
10100   CONTINUE
C
      RETURN
      END
C
C FITEST
C rigid fit test driver
C
      PROGRAM fitest
      IMPLICIT CHARACTER (A-Z)
C
      INTEGER n
      DOUBLEPRECISION angle
      DOUBLEPRECISION pert
      PARAMETER (n = 401)
      PARAMETER (angle = 1.0D0)
      PARAMETER (pert = 0.2D0)
C
      DOUBLEPRECISION y (3, n)
      DOUBLEPRECISION x (3, n)
      DOUBLEPRECISION z (3, n)
      DOUBLEPRECISION w (n)
      DOUBLEPRECISION o (3)
      DOUBLEPRECISION q (4)
      DOUBLEPRECISION u (3, 3)
      INTEGER ierr
C
C seed the random number generator
C call ranini for vax
C
      CALL ranset (123571)
*     CALL ranini ()
C
C generate reference and test molecules and weights
C
      CALL genxyw (angle, pert, n, x, y, w)
      CALL center (n, y, w, .FALSE., o, o)
      CALL center (n, x, w, .FALSE., o, o)
      CALL identm (3, 3, u)
      CALL analyz (n, x, y, w, u)
C
C fit x to y
C
      WRITE (*, *)
      WRITE (*, *) 'fitting x to y'
      CALL qtrfit (n, x, y, w, q, u, ierr)
      CALL rotmol (n, x, z, u)
      CALL analyz (n, z, y, w, u)
C
C fit y to x
C
      WRITE (*, *)
      WRITE (*, *) 'fitting y to x'
      CALL qtrfit (n, y, x, w, q, u, ierr)
      CALL rotmol (n, y, z, u)
      CALL analyz (n, z, x, w, u)
C
      END
C
C GENROT
C Generate a left rotation matrix from a normalized rotation axis
C and cosine and sine of the rotation angle.
C
C INPUT
C   a      - rotation axis
C   cost   - cosine of rotation angle
C   sint   - sine of rotation angle
C
C OUTPUT
C   u      - the rotation matrix
C
      SUBROUTINE genrot (a, cost, sint, u)
      IMPLICIT CHARACTER (A-Z)
C
      DOUBLEPRECISION a (3)
      DOUBLEPRECISION cost, sint
      DOUBLEPRECISION u (3, 3)
C
      u (1, 1) = a (1) ** 2 + (1.0D0 - a (1) ** 2) * cost
      u (2, 1) = a (1) * a (2) * (1.0D0 - cost) + a (3) * sint
      u (3, 1) = a (1) * a (3) * (1.0D0 - cost) - a (2) * sint
C
      u (1, 2) = a (1) * a (2) * (1.0D0 - cost) - a (3) * sint
      u (2, 2) = a (2) ** 2 + (1.0D0 - a (2) ** 2) * cost
      u (3, 2) = a (2) * a (3) * (1.0D0 - cost) + a (1) * sint
C
      u (1, 3) = a (1) * a (3) * (1.0D0 - cost) + a (2) * sint
      u (2, 3) = a (2) * a (3) * (1.0D0 - cost) - a (1) * sint
      u (3, 3) = a (3) ** 2 + (1.0D0 - a (3) ** 2) * cost
C
      RETURN
      END
C
C GENXYW
C generate random reference and test molecules and weights
C
      SUBROUTINE genxyw (angle, pert, n, x, y, w)
      IMPLICIT CHARACTER (A-Z)
C
      DOUBLEPRECISION angle
      DOUBLEPRECISION pert
      INTEGER n
      DOUBLEPRECISION x (3, n)
      DOUBLEPRECISION y (3, n)
      DOUBLEPRECISION w (n)
C
      DOUBLEPRECISION u (3, 3)            
      INTEGER i
C
      DOUBLEPRECISION random
      EXTERNAL random
C
      CALL ranmol (n, y)
      CALL ranrot (angle, u)
      CALL rotmol (n, y, x, u)
      CALL pertur (pert, n, x)
      CALL ranrot (angle, u)
      CALL rotmol (n, x, x, u)
C
      DO 10000 i = 1, n
        w (i) = random ()
10000   CONTINUE
C
      RETURN
      END
C
C IDENTM
C generate an identity matrix
C
      SUBROUTINE identm (n, np, u)
      IMPLICIT CHARACTER (A-Z)
C
      INTEGER n, np
      DOUBLEPRECISION u (np, n)
C
      INTEGER i, j
C
      DO 10000 j = 1, n
        DO 10010 i = 1, n
          u (i, j) = 0.0D0
10010     CONTINUE
        u (j, j) = 1.0D0
10000   CONTINUE
C
      RETURN
      END
C
C JACOBI
C Jacobi diagonalizer with sorted output.  Same calling sequence as
C EISPACK routine, but must specify nrot!
C
      SUBROUTINE jacobi (a, n, np, d, v, nrot)
      IMPLICIT CHARACTER (A-Z)
C
      INTEGER n, np, nrot
      DOUBLEPRECISION a (np, n)
      DOUBLEPRECISION d (n)
      DOUBLEPRECISION v (np, n)
C
      DOUBLEPRECISION onorm, dnorm
      DOUBLEPRECISION b, dma, q, t, c, s
      DOUBLEPRECISION atemp, vtemp, dtemp
      INTEGER i, j, k, l
C
      DO 10000 j = 1, n
        DO 10010 i = 1, n
          v (i, j) = 0.0D0
10010     CONTINUE
        v (j, j) = 1.0D0
        d (j) = a (j, j)
10000   CONTINUE
C
      DO 20000 l = 1, nrot
        dnorm = 0.0D0
        onorm = 0.0D0
        DO 20100 j = 1, n
          dnorm = dnorm + ABS (d (j))
          DO 20110 i = 1, j - 1
            onorm = onorm + ABS (a (i, j))
20110       CONTINUE
20100     CONTINUE
        IF (onorm / dnorm .LE. 0.0D0) GOTO 19999
        DO 21000 j = 2, n
          DO 21010 i = 1, j - 1
            b = a (i, j)
            IF (ABS (b) .GT. 0.0D0) THEN
              dma = d (j) - d (i)
              IF (ABS (dma) + ABS (b) .LE. ABS (dma)) THEN
                t = b / dma
              ELSE
                q = 0.5D0 * dma / b
                t = SIGN (1.0D0 / (ABS (q) + SQRT (1.0D0 + q * q)), q)
              ENDIF
              c = 1.0D0 / SQRT (t * t + 1.0D0)
              s = t * c
              a (i, j) = 0.0D0
              DO 21110 k = 1, i - 1
                atemp    = c * a (k, i) - s * a (k, j)
                a (k, j) = s * a (k, i) + c * a (k, j)
                a (k, i) = atemp
21110           CONTINUE
              DO 21120 k = i + 1, j - 1
                atemp    = c * a (i, k) - s * a (k, j)
                a (k, j) = s * a (i, k) + c * a (k, j)
                a (i, k) = atemp
21120           CONTINUE
              DO 21130 k = j + 1, n
                atemp    = c * a (i, k) - s * a (j, k)
                a (j, k) = s * a (i, k) + c * a (j, k)
                a (i, k) = atemp
21130           CONTINUE
              DO 21140 k = 1, n
                vtemp    = c * v (k, i) - s * v (k, j)
                v (k, j) = s * v (k, i) + c * v (k, j)
                v (k, i) = vtemp
21140           CONTINUE
              dtemp = c * c * d (i) + s * s * d (j) -
     &                2.0D0 * c * s * b
              d (j) = s * s * d (i) + c * c * d (j) +
     &                2.0D0 * c * s * b
              d (i) = dtemp
              ENDIF
21010       CONTINUE
21000     CONTINUE
20000   CONTINUE
19999 CONTINUE
      nrot = l
C
      DO 30000 j = 1, n - 1
        k = j
        dtemp = d (k)
        DO 30100 i = j + 1, n
          IF (d (i) .LT. dtemp) THEN
            k = i
            dtemp = d (k)
            ENDIF
30100     CONTINUE
        IF (k .GT. j) THEN
          d (k) = d (j)
          d (j) = dtemp
          DO 30200 i = 1, n
            dtemp    = v (i, k)
            v (i, k) = v (i, j)
            v (i, j) = dtemp
30200       CONTINUE
          ENDIF
30000   CONTINUE
C
      RETURN
      END
C
C PERTUR
C apply a random perturbation
C
      SUBROUTINE pertur (pert, n, x)
      IMPLICIT CHARACTER (A-Z)
C
      DOUBLEPRECISION pert
      INTEGER n
      DOUBLEPRECISION x (3, n)
C
      INTEGER i, j
C
      DOUBLEPRECISION random
      EXTERNAL random
C
      DO 10000 j = 1, 3
        DO 10010 i = 1, n
          x (j, i) = x (j, i) *
     &               (1.0D0 + 2.0D0 * pert * (0.5D0 - random ()))
10010     CONTINUE
10000   CONTINUE
C
      RETURN
      END
C
C Q2MAT
C Generate a left rotation matrix from a normalized quaternion
C
C INPUT
C   q      - normalized quaternion
C
C OUTPUT
C   u      - the rotation matrix
C
      SUBROUTINE q2mat (q, u)
      IMPLICIT NONE
C
      DOUBLEPRECISION q (0 : 3)
      DOUBLEPRECISION u (3, 3)
C
      u (1, 1) = q (0) ** 2 + q (1) ** 2 - q (2) ** 2 - q (3) ** 2
      u (2, 1) = 2.0D0 * (q (1) * q (2) - q (0) * q (3))
      u (3, 1) = 2.0D0 * (q (1) * q (3) + q (0) * q (2))
C
      u (1, 2) = 2.0D0 * (q (2) * q (1) + q (0) * q (3))
      u (2, 2) = q (0) ** 2 - q (1) ** 2 + q (2) ** 2 - q (3) ** 2
      u (3, 2) = 2.0D0 * (q (2) * q (3) - q (0) * q (1))
C
      u (1, 3) = 2.0D0 * (q (3) * q (1) - q (0) * q (2))
      u (2, 3) = 2.0D0 * (q (3) * q (2) + q (0) * q (1))
      u (3, 3) = q (0) ** 2 - q (1) ** 2 - q (2) ** 2 + q (3) ** 2
C
      RETURN
      END
C
C QTRFIT
C Find the quaternion, q, [and left rotation matrix, u] that minimizes
C
C   |qTXq - Y| ^ 2    [|uX - Y| ^ 2]
C
C This is equivalent to maximizing Re (qTXTqY).
C
C This is equivalent to finding the largest eigenvalue and corresponding
C eigenvector of the matrix
C
C   [A2   AUx  AUy  AUz ]
C   [AUx  Ux2  UxUy UzUx]
C   [AUy  UxUy Uy2  UyUz]
C   [AUz  UzUx UyUz Uz2 ]
C
C where
C
C   A2   = Xx Yx + Xy Yy + Xz Yz
C   Ux2  = Xx Yx - Xy Yy - Xz Yz
C   Uy2  = Xy Yy - Xz Yz - Xx Yx
C   Uz2  = Xz Yz - Xx Yx - Xy Yy
C   AUx  = Xz Yy - Xy Yz
C   AUy  = Xx Yz - Xz Yx
C   AUz  = Xy Yx - Xx Yy
C   UxUy = Xx Yy + Xy Yx
C   UyUz = Xy Yz + Xz Yy
C   UzUx = Xz Yx + Xx Yz
C
C The left rotation matrix, u, is obtained from q by
C
C   u = qT1q
C
C INPUT
C   n      - number of points
C   x      - test vector
C   y      - reference vector
C   w      - weight vector
C
C OUTPUT
C   q      - the best-fit quaternion
C   u      - the best-fit left rotation matrix
C   nr     - number of jacobi sweeps required
C
      SUBROUTINE qtrfit (n, x, y, w, q, u, nr)
      IMPLICIT CHARACTER (A-Z)
C
      INTEGER n
      DOUBLEPRECISION x (3, n)
      DOUBLEPRECISION y (3, n)
      DOUBLEPRECISION w (n)
      DOUBLEPRECISION q (0 : 3)
      DOUBLEPRECISION u (3, 3)
      INTEGER nr
C
      DOUBLEPRECISION xxyx, xxyy, xxyz
      DOUBLEPRECISION xyyx, xyyy, xyyz
      DOUBLEPRECISION xzyx, xzyy, xzyz
      DOUBLEPRECISION c (0 : 3, 0 : 3), v (0 : 3, 0 : 3)
      DOUBLEPRECISION d (0 : 3)
      INTEGER i
C
C generate the upper triangle of the quadratic form matrix
C
      xxyx = 0.0D0
      xxyy = 0.0D0
      xxyz = 0.0D0
      xyyx = 0.0D0
      xyyy = 0.0D0
      xyyz = 0.0D0
      xzyx = 0.0D0
      xzyy = 0.0D0
      xzyz = 0.0D0
      DO 11000 i = 1, n
        xxyx = xxyx + x (1, i) * y (1, i) * w (i)
        xxyy = xxyy + x (1, i) * y (2, i) * w (i)
        xxyz = xxyz + x (1, i) * y (3, i) * w (i)
        xyyx = xyyx + x (2, i) * y (1, i) * w (i)
        xyyy = xyyy + x (2, i) * y (2, i) * w (i)
        xyyz = xyyz + x (2, i) * y (3, i) * w (i)
        xzyx = xzyx + x (3, i) * y (1, i) * w (i)
        xzyy = xzyy + x (3, i) * y (2, i) * w (i)
        xzyz = xzyz + x (3, i) * y (3, i) * w (i)
11000   CONTINUE
C
      c (0, 0) = xxyx + xyyy + xzyz
C
      c (0, 1) = xzyy - xyyz
      c (1, 1) = xxyx - xyyy - xzyz
C
      c (0, 2) = xxyz - xzyx
      c (1, 2) = xxyy + xyyx
      c (2, 2) = xyyy - xzyz - xxyx
C
      c (0, 3) = xyyx - xxyy
      c (1, 3) = xzyx + xxyz
      c (2, 3) = xyyz + xzyy
      c (3, 3) = xzyz - xxyx - xyyy
C
C diagonalize c
C
      nr = 16
      CALL jacobi (c, 4, 4, d, v, nr)
C
C extract the desired quaternion
C
      q (0) = v (0, 3)
      q (1) = v (1, 3)
      q (2) = v (2, 3)
      q (3) = v (3, 3)
C
C generate the rotation matrix
C
      CALL q2mat (q, u)
C
      RETURN
      END
C
C RANDOM
C random number generator after Knuth
C
      DOUBLEPRECISION FUNCTION random ()
      IMPLICIT CHARACTER (A-Z)
C
      INTEGER ntable
      INTEGER im1, im2, im3
      INTEGER ia1, ia2, ia3
      INTEGER ic1, ic2, ic3
      DOUBLEPRECISION rm1, rm2, rm3
      PARAMETER (ntable = 97)
      PARAMETER (im1 = 714025, im2 = 214326, im3 = 139968)
      PARAMETER (ia1 = 1366,   ia2 = 3613,   ia3 = 3877)
      PARAMETER (ic1 = 150889, ic2 = 45289,  ic3 = 29573)
      PARAMETER (rm1 = 1.0D0 / im1)
      PARAMETER (rm2 = 1.0D0 / im2)
      PARAMETER (rm3 = 1.0D0 / im3)
C
      INTEGER iseed0, iseed1, iseed2, iseed3
      DOUBLEPRECISION r (ntable)
      INTEGER i
C
      COMMON /rtable/ iseed0, iseed1, iseed2, iseed3, r
      SAVE /rtable/
C
      iseed1 = MOD (ia1 * iseed1 + ic1, im1)
      iseed2 = MOD (ia2 * iseed2 + ic2, im2)
      iseed3 = MOD (ia3 * iseed3 + ic3, im3)
C
      i = 1 + (ntable * iseed3) * rm3
      random = r (i)
      r (i) = (iseed1 + iseed2 * rm2) * rm1
C
      RETURN
      END
      INTEGER FUNCTION ranget ()
      IMPLICIT CHARACTER (A-Z)
C
      INTEGER ntable
      INTEGER im1, im2, im3
      INTEGER ia1, ia2, ia3
      INTEGER ic1, ic2, ic3
      DOUBLEPRECISION rm1, rm2, rm3
      PARAMETER (ntable = 97)
      PARAMETER (im1 = 714025, im2 = 214326, im3 = 139968)
      PARAMETER (ia1 = 1366,   ia2 = 3613,   ia3 = 3877)
      PARAMETER (ic1 = 150889, ic2 = 45289,  ic3 = 29573)
      PARAMETER (rm1 = 1.0D0 / im1)
      PARAMETER (rm2 = 1.0D0 / im2)
      PARAMETER (rm3 = 1.0D0 / im3)
C
      INTEGER iseed0, iseed1, iseed2, iseed3
      DOUBLEPRECISION r (ntable)
C
      COMMON /rtable/ iseed0, iseed1, iseed2, iseed3, r
      SAVE /rtable/
C
      ranget = iseed0
C
      RETURN
      END
C
C RANINI
C seed the random number generator
C
*     SUBROUTINE ranini ()
*     IMPLICIT CHARACTER (A-Z)
C
*     INTEGER bintim (2)
C
*     CALL sys$gettim (bintim)
*     CALL ranset (IAND (bintim (1), '7FFFFFFF'X))
C
*     RETURN
*     END 
C
C RANMOL
C generate a random molecule
C
      SUBROUTINE ranmol (n, x)
      IMPLICIT CHARACTER (A-Z)
C
      INTEGER n
      DOUBLEPRECISION x (3, n)
C
      INTEGER i, j
C
      DOUBLEPRECISION random
      EXTERNAL random
C
C use nested loops to get same result for scalar/vector
C
      DO 10000 j = 1, 3
        DO 10010 i = 1, n
          x (j, i) = random ()
10010     CONTINUE
10000   CONTINUE
C
      RETURN
      END
C
C RANROT
C generate a random (almost) rotation matrix
C
      SUBROUTINE ranrot (angle, u)
      IMPLICIT CHARACTER (A-Z)
C
      DOUBLEPRECISION angle
      DOUBLEPRECISION u (3, 3)
C
      DOUBLEPRECISION a (3)
      DOUBLEPRECISION anorm
      DOUBLEPRECISION theta
      DOUBLEPRECISION pi
C
      DOUBLEPRECISION random
      EXTERNAL random
C
      pi = 4.0D0 * ATAN (1.0D0)
C
      a (1) = 0.5D0 - random ()
      a (2) = 0.5D0 - random ()
      a (3) = 0.5D0 - random ()
C
      anorm = SQRT (a (1) ** 2 + a (2) ** 2 + a (3) ** 2)
      a (1) = a (1) / anorm
      a (2) = a (2) / anorm
      a (3) = a (3) / anorM
C
      theta = angle * pi * random ()
C
      CALL genrot (a, COS (theta), SIN (theta), u)
C
      RETURN
      END
      SUBROUTINE ranset (iseed)
      IMPLICIT CHARACTER (A-Z)
C
      INTEGER ntable
      INTEGER im1, im2, im3
      INTEGER ia1, ia2, ia3
      INTEGER ic1, ic2, ic3
      DOUBLEPRECISION rm1, rm2, rm3
      PARAMETER (ntable = 97)
      PARAMETER (im1 = 714025, im2 = 214326, im3 = 139968)
      PARAMETER (ia1 = 1366,   ia2 = 3613,   ia3 = 3877)
      PARAMETER (ic1 = 150889, ic2 = 45289,  ic3 = 29573)
      PARAMETER (rm1 = 1.0D0 / im1)
      PARAMETER (rm2 = 1.0D0 / im2)
      PARAMETER (rm3 = 1.0D0 / im3)
C
      INTEGER iseed
C
      INTEGER iseed0, iseed1, iseed2, iseed3
      DOUBLEPRECISION r (ntable)
      INTEGER i
C
      COMMON /rtable/ iseed0, iseed1, iseed2, iseed3, r
      SAVE /rtable/
C
      iseed0 = iseed
      iseed1 = MOD (iseed0, im1)
      iseed2 = MOD (iseed0, im2)
      iseed3 = MOD (iseed0, im3)
C
      DO 10000 i = 1, ntable
        iseed1 = MOD (ia1 * iseed1 + ic1, im1)
        iseed2 = MOD (ia2 * iseed2 + ic2, im2)
        iseed3 = MOD (ia3 * iseed3 + ic3, im3)
        r (i) = (iseed1 + iseed2 * rm2) * rm1
10000   CONTINUE
C
      RETURN
      END
C
C ROTMOL
C rotate a molecule
C
      SUBROUTINE rotmol (n, x, y, u)
      IMPLICIT CHARACTER (A-Z)
C
      INTEGER n
      DOUBLEPRECISION x (3, n)
      DOUBLEPRECISION y (3, n)
      DOUBLEPRECISION u (3, 3)
C
      DOUBLEPRECISION yx, yy, yz
      INTEGER i
C
      DO 10000 i = 1, n
        yx = u(1, 1) * x(1, i) + u(1, 2) * x(2, i) + u(1, 3) * x(3, i)
        yy = u(2, 1) * x(1, i) + u(2, 2) * x(2, i) + u(2, 3) * x(3, i)
        yz = u(3, 1) * x(1, i) + u(3, 2) * x(2, i) + u(3, 3) * x(3, i)
C
        y (1, i) = yx
        y (2, i) = yy
        y (3, i) = yz
10000   CONTINUE
C
      RETURN
      END