* Last processed by NICE on 12-Jun-2000 15:54:00
* Customized for :  IEEE, LINUX, UNIX, MOTIF, F77
      SUBROUTINE FIT_QUADRATIC (X,Y,Z,W,N,F)
!
! where z = f(1)+f(2)*x+f(3)*y+f(4)*x^2+f(5)*y^2
!
      INTEGER I, J, N, NORDER
      PARAMETER (NORDER = 5)
      REAL*8 G(NORDER,NORDER), H(NORDER), DET
      REAL*4 X(N), Y(N), Z(N), W(N), F(NORDER)
!
      DO I=1,NORDER
         H(I)=0
         DO J=1,NORDER
            G(J,I) = 0
         ENDDO
      ENDDO
      DO I=1,N
         G(1,1) = G(1,1) + W(I)
         G(1,2) = G(1,2) + W(I)*X(I)
         G(1,3) = G(1,3) + W(I)*Y(I)
         G(1,4) = G(1,4) + W(I)*X(I)*X(I)
         G(1,5) = G(1,5) + W(I)*Y(I)*Y(I)
         G(2,3) = G(2,3) + W(I)*X(I)*Y(I)
         G(2,4) = G(2,4) + W(I)*X(I)*X(I)*X(I)
         G(2,5) = G(2,5) + W(I)*Y(I)*Y(I)*X(I)
         G(3,3) = G(3,3) + W(I)*Y(I)*Y(I)*Y(I)
         G(3,4) = G(3,4) + W(I)*X(I)*X(I)*Y(I)
         G(3,5) = G(3,5) + W(I)*Y(I)*Y(I)*Y(I)*Y(I)
         G(4,4) = G(4,4) + W(I)*X(I)*X(I)*X(I)*X(I)
         G(4,5) = G(4,5) + W(I)*X(I)*X(I)*Y(I)*Y(I)
         G(5,5) = G(5,5) + W(I)*Y(I)*Y(I)*Y(I)*Y(I)
         H(1)   = H(1)   + W(I)*Z(I)
         H(2)   = H(2)   + W(I)*Z(I)*X(I)
         H(3)   = H(3)   + W(I)*Z(I)*Y(I)
         H(4)   = H(4)   + W(I)*Z(I)*X(I)*X(I)
         H(5)   = H(5)   + W(I)*Z(I)*Y(I)*Y(I)
      ENDDO
      G(2,1) = G(1,2)
      G(2,2) = G(1,4)
      G(3,1) = G(1,3)
      G(3,2) = G(2,3)
      G(4,1) = G(1,4)
      G(4,2) = G(2,4)
      G(4,3) = G(3,4)
      G(5,1) = G(1,5)
      G(5,2) = G(2,5)
      G(5,3) = G(3,5)
      G(5,4) = G(4,5)
!
      CALL MATINV (G,NORDER,DET)
!
      DO I=1,NORDER
         F(I) = 0
         DO J=1,NORDER
            F(I) = G(I,J)*H(J)+F(I)
         ENDDO
      ENDDO
      RETURN
      END
!
!
c------------------------------------------------------------------------------
c	purpose
c		invert a symmetric matrix and calculate its determinant
c
c	usage
c		call matinv (array, norder, det)
c
c	description of parameters
c		array	- input matrix which is replaced by its inverse
c		norder	- degree of matrix (order of determinant)
c		det	- determinant of input matrix
c
c	subroutines and functions required
c		none
c
c	valid for up to "maxpar" parameters.
c
c------------------------------------------------------------------------------
      SUBROUTINE MATINV (ARRAY, NORDER, DET)
!
      INTEGER   NORDER
      REAL*8    DET
      REAL*8    ARRAY (NORDER,NORDER)
!
! local data
      REAL*8    AMAX, SAVE
      INTEGER   I,J,K,L
      INTEGER    MAXPAR
      PARAMETER (MAXPAR = 5)
      INTEGER    IK(MAXPAR), JK(MAXPAR)
!
      DET = 1
      DO K = 1,NORDER
!
c -----	find largest element array(i,j) in rest of matrix.
!
         AMAX  = 0
100      DO I = K,NORDER
            DO J = K,NORDER
               IF (DABS(AMAX).LE.DABS(ARRAY(I,J))) THEN
                  AMAX  = ARRAY(I,J)
                  IK(K) = I
                  JK(K) = J
               ENDIF
            ENDDO
         ENDDO
!
c       -----	interchange rows and columns to put amax in array(k,k)
!
         IF (AMAX.NE.0) THEN
            I = IK(K)
            IF (I.GE.K) THEN
               IF (I.GT.K) THEN
                  DO J = 1,NORDER
                     SAVE = ARRAY(K,J)
                     ARRAY(K,J) = ARRAY(I,J)
                     ARRAY(I,J) = -SAVE
                  ENDDO
               ENDIF
               J = JK(K)
               IF (J.GE.K) THEN
                  IF (J.GT.K) THEN
                     DO I = 1,NORDER
                        SAVE = ARRAY(I,K)
                        ARRAY(I,K) = ARRAY(I,J)
                        ARRAY(I,J) = -SAVE
                     ENDDO
                  ENDIF
!
c       -----	accumulate elements of inverse matrix
!
                  DO I = 1,NORDER
                     IF (I.NE.K) ARRAY(I,K) = -ARRAY(I,K)/AMAX
                  ENDDO
                  DO I = 1,NORDER
                     DO J = 1,NORDER
                        IF (I.NE.K.AND.J.NE.K) ARRAY(I,J) = ARRAY(I,J)+
     $                  ARRAY(I,K)*ARRAY(K,J)
                     ENDDO
                  ENDDO
                  DO J = 1,NORDER
                     IF (J.NE.K) ARRAY(K,J) = ARRAY(K,J)/AMAX
                  ENDDO
                  ARRAY(K,K) = 1/AMAX
                  DET = DET*AMAX
               ELSE
                  GOTO 100
               ENDIF
            ELSE
               GOTO 100
            ENDIF
         ELSE
            DET = 0
            RETURN
         ENDIF
      ENDDO
!
c       -----	restore ordering of matrix
!
      DO L = 1,NORDER
         K = NORDER-L+1
         J = IK(K)
         IF (J.GT.K) THEN
            DO I = 1,NORDER
               SAVE = ARRAY(I,K)
               ARRAY(I,K) = -ARRAY(I,J)
               ARRAY(I,J) = SAVE
            ENDDO
         ENDIF
         I = JK(K)
         IF (I.GT.K) THEN
            DO J = 1,NORDER
               SAVE = ARRAY(K,J)
               ARRAY(K,J) = -ARRAY(I,J)
               ARRAY(I,J) = SAVE
            ENDDO
         ENDIF
      ENDDO
      RETURN
      END
