      SUBROUTINE PLANE(COORD,NPOINT,CENTRE,NORMAL,EIGVEC,ORDVEC,IFAIL)  PLAN001
                                                                        PLAN002
C-----------------------------------------------------------------------PLAN003
C                                                                       PLAN004
C  SUBROUTINE PLANE  -  Calculate a best-fit plane through the supplied PLAN005
C                       data points                                     PLAN006
C                                                                       PLAN007
C----------------------------------------------------------------------+PLAN008
                                                                        PLAN009
      INTEGER       NPOINT, NROT, ORDVEC(3)                             PLAN010
      LOGICAL       IFAIL                                               PLAN011
      REAL          CENTRE(3), COORD(3,NPOINT), D, EIGVEC(3,3),         PLAN012
     &              EIGVAL(3), MATRIX(3,3), NORMAL(3)                   PLAN013
                                                                        PLAN014
C---- Calculate centre of mass of points and adjust coords such that    PLAN015
C     CofM is at the origin                                             PLAN016
      CALL GETCOM(COORD,NPOINT,CENTRE)                                  PLAN017
                                                                        PLAN018
C---- Calculate scalar products matrix                                  PLAN019
      CALL SCPMAT(COORD,NPOINT,MATRIX)                                  PLAN020
                                                                        PLAN021
C---- Calculate eigen values and eigenvectors of matrix                 PLAN022
      CALL JACOBI(MATRIX,3,3,EIGVAL,EIGVEC,IFAIL,NROT)                  PLAN023
      IF (IFAIL) GO TO 900                                              PLAN024
                                                                        PLAN025
C---- Sort the eigen vectors into descending order of eigen value       PLAN026
      CALL EIGSRT(EIGVAL,ORDVEC)                                        PLAN027
                                                                        PLAN028
C---- Compute parameters of best-fit plane                              PLAN029
      NORMAL(1) = EIGVEC(1,ORDVEC(3))                                   PLAN030
      NORMAL(2) = EIGVEC(2,ORDVEC(3))                                   PLAN031
      NORMAL(3) = EIGVEC(3,ORDVEC(3))                                   PLAN032
      D = NORMAL(1) * CENTRE(1) + NORMAL(2) * CENTRE(2)                 PLAN033
     &    + NORMAL(3) * CENTRE(3)                                       PLAN034
                                                                        PLAN035
      GO TO 999                                                         PLAN036
                                                                        PLAN037
C---- Error messages                                                    PLAN038
 900  CONTINUE                                                          PLAN039
      PRINT*, '**** Error. Failed to find eigenvalues of matrix'        PLAN040
      GO TO 990                                                         PLAN041
                                                                        PLAN042
 990  CONTINUE                                                          PLAN043
      IFAIL = .TRUE.                                                    PLAN044
                                                                        PLAN045
 999  CONTINUE                                                          PLAN046
      RETURN                                                            PLAN047
      END                                                               PLAN048
