      FUNCTION STORE (X)
      REAL X
C
C***********************************************************
C
C                                    ROBERT RENKA
C                             NORTH TEXAS STATE UNIV.
C                                    (817) 565-2767
C
C   THIS FUNCTION FORCES ITS ARGMENT X TO BE STORED IN MAIN
C MEMORY.  THIS IS USEFUL FOR COMPUTING MACHINE DEPENDENT
C PARAMETERS (SUCH AS THE MACHINE PRECISION) WHERE IT IS
C NECESSARY TO AVOID COMPUTATION IN HIGH PRECISION REG-
C ISTERS.
C
C INPUT PARAMETER -
C
C      X - VALUE TO BE STORED.
C
C X IS NOT ALTERED BY THIS FUNCTION.
C
C OUTPUT PARAMETER -
C
C      STORE - VALUE OF X AFTER IT HAS BEEN STORED AND
C            (POSSIBLY) TRUNCATED OR ROUNDED TO THE
C            MACHINE WORD LENGTH.
C
C MODULES REQUIRED BY STORE - NONE
C
C***********************************************************
C
      COMMON/STCOM/Y
      Y = X
      STORE = Y
      RETURN
      END
      SUBROUTINE QSHEP2 (N,X,Y,F,NQ,NW,NR, LCELL,LNEXT,XMIN,
     &  YMIN,DX,DY,RMAX,RSQ,A,IER)
      INTEGER N, NQ, NW, NR, LCELL(NR,NR), LNEXT(N), IER
      REAL    X(N), Y(N), F(N), XMIN, YMIN, DX, DY, RMAX,
     &  RSQ(N), A(5,N)
C
C***********************************************************
C
C                                    ROBERT RENKA
C                              UNIV. OF NORTH TEXAS
C                                    (817) 565-2767
C                                        01/08/90
C
C   THIS SUBROUTINE COMPUTES A SET OF PARAMETERS A AND RSQ
C DEFINING A SMOOTH (ONCE CONTINUOUSLY DIFFERENTIABLE) BI-
C VARIATE FUNCTION Q(X,Y) WHICH INTERPOLATES DATA VALUES F
C AT SCATTERED NODES (X,Y).  THE INTERPOLANT Q MAY BE EVAL-
C UATED AT AN ARBITRARY POINT BY FUNCTION QS2VAL, AND ITS
C FIRST DERIVATIVES ARE COMPUTED BY SUBROUTINE QS2GRD.
C   THE INTERPOLATION SCHEME IS A MODIFIED QUADRATIC SHEPARD
C METHOD --
C
C Q = (W(1)*Q(1)+W(2)*Q(2)+..+W(N)*Q(N))/(W(1)+W(2)+..+W(N))
C
C FOR BIVARIATE FUNCTIONS W(K) AND Q(K).  THE NODAL FUNC-
C TIONS ARE GIVEN BY
C
C Q(K)(X,Y) = A(1,K)*(X-X(K))**2 + A(2,K)*(X-X(K))*(Y-Y(K))
C          + A(3,K)*(Y-Y(K))**2 + A(4,K)*(X-X(K))
C          + A(5,K)*(Y-Y(K))       + F(K) .
C
C THUS, Q(K) IS A QUADRATIC FUNCTION WHICH INTERPOLATES THE
C DATA VALUE AT NODE K.  ITS COEFFICIENTS A(,K) ARE OBTAINED
C BY A WEIGHTED LEAST SQUARES FIT TO THE CLOSEST NQ DATA
C POINTS WITH WEIGHTS SIMILAR TO W(K).      NOTE THAT THE RADIUS
C OF INFLUENCE FOR THE LEAST SQUARES FIT IS FIXED FOR EACH
C K, BUT VARIES WITH K.
C   THE WEIGHTS ARE TAKEN TO BE
C
C W(K)(X,Y) = ( (R(K)-D(K))+ / R(K)*D(K) )**2
C
C WHERE (R(K)-D(K))+ = 0 IF R(K) .LE. D(K) AND D(K)(X,Y) IS
C THE EUCLIDEAN DISTANCE BETWEEN (X,Y) AND (X(K),Y(K)).  THE
C RADIUS OF INFLUENCE R(K) VARIES WITH K AND IS CHOSEN SO
C THAT NW NODES ARE WITHIN THE RADIUS.      NOTE THAT W(K) IS
C NOT DEFINED AT NODE (X(K),Y(K)), BUT Q(X,Y) HAS LIMIT F(K)
C AS (X,Y) APPROACHES (X(K),Y(K)).
C
C ON INPUT --
C
C      N = NUMBER OF NODES AND ASSOCIATED DATA VALUES.
C          N .GE. 6.
C
C      X,Y = ARRAYS OF LENGTH N CONTAINING THE CARTESIAN
C            COORDINATES OF THE NODES.
C
C      F = ARRAY OF LENGTH N CONTAINING THE DATA VALUES
C          IN ONE-TO-ONE CORRESPONDENCE WITH THE NODES.
C
C      NQ = NUMBER OF DATA POINTS TO BE USED IN THE LEAST
C           SQUARES FIT FOR COEFFICIENTS DEFINING THE NODAL
C           FUNCTIONS Q(K).  A HIGHLY RECOMMENDED VALUE IS
C           NQ = 13.  5 .LE. NQ .LE. MIN(40,N-1).
C
C      NW = NUMBER OF NODES WITHIN (AND DEFINING) THE RADII
C           OF INFLUENCE R(K) WHICH ENTER INTO THE WEIGHTS
C           W(K).  FOR N SUFFICIENTLY LARGE, A RECOMMENDED
C           VALUE IS NW = 19.      1 .LE. NW .LE. MIN(40,N-1).
C
C      NR = NUMBER OF ROWS AND COLUMNS IN THE CELL GRID DE-
C           FINED IN SUBROUTINE STORE2.  A RECTANGLE CON-
C           TAINING THE NODES IS PARTITIONED INTO CELLS IN
C           ORDER TO INCREASE SEARCH EFFICIENCY.  NR =
C           SQRT(N/3) IS RECOMMENDED.      NR .GE. 1.
C
C THE ABOVE PARAMETERS ARE NOT ALTERED BY THIS ROUTINE.
C
C      LCELL = ARRAY OF LENGTH .GE. NR**2.
C
C      LNEXT = ARRAY OF LENGTH .GE. N.
C
C      RSQ = ARRAY OF LENGTH .GE. N.
C
C      A = ARRAY OF LENGTH .GE. 5N.
C
C ON OUTPUT --
C
C      LCELL = NR BY NR ARRAY OF NODAL INDICES ASSOCIATED
C            WITH CELLS.  REFER TO STORE2.
C
C      LNEXT = ARRAY OF LENGTH N CONTAINING NEXT-NODE INDI-
C            CES.  REFER TO STORE2.
C
C      XMIN,YMIN,DX,DY = MINIMUM NODAL COORDINATES AND CELL
C                    DIMENSIONS.  REFER TO STORE2.
C
C      RMAX = SQUARE ROOT OF THE LARGEST ELEMENT IN RSQ --
C             MAXIMUM RADIUS R(K).
C
C      RSQ = ARRAY CONTAINING THE SQUARES OF THE RADII R(K)
C            WHICH ENTER INTO THE WEIGHTS W(K).
C
C      A = 5 BY N ARRAY CONTAINING THE COEFFICIENTS FOR
C          QUADRATIC NODAL FUNCTION Q(K) IN COLUMN K.
C
C   NOTE THAT THE ABOVE OUTPUT PARAMETERS ARE NOT DEFINED
C UNLESS IER = 0.
C
C      IER = ERROR INDICATOR --
C            IER = 0 IF NO ERRORS WERE ENCOUNTERED.
C            IER = 1 IF N, NQ, NW, OR NR IS OUT OF RANGE.
C            IER = 2 IF DUPLICATE NODES WERE ENCOUNTERED.
C            IER = 3 IF ALL NODES ARE COLLINEAR.
C
C MODULES REQUIRED BY QSHEP2 -- GETNP2, GIVENS, ROTATE,
C                          SETUP2, STORE2
C
C INTRINSIC FUNCTIONS CALLED BY QSHEP2 -- ABS, AMIN1, FLOAT,
C                                  MAX0, MIN0, SQRT
C
C***********************************************************
C
      INTEGER I, IB, IERR, IP1, IRM1, IROW, J, JP1, K, LMAX,
     &  LNP, NEQ, NN, NNQ, NNR, NNW, NP, NPTS(40),
     &  NQWMAX
      REAL    AV, AVSQ, B(6,6), C, DDX, DDY, DMIN, DTOL, FK,
     &  RQ, RS, RSMX, RSOLD, RTOL, RWS, S, SF, SUM, T,
     &  XK, XMN, YK, YMN
C
      DATA    RTOL/1.E-5/, DTOL/.01/, SF/1./
C
C LOCAL PARAMETERS --
C
C AV =             ROOT-MEAN-SQUARE DISTANCE BETWEEN K AND THE
C             NODES IN THE LEAST SQUARES SYSTEM (UNLESS
C             ADDITIONAL NODES ARE INTRODUCED FOR STABIL-
C             ITY).      THE FIRST 3 COLUMNS OF THE MATRIX
C             ARE SCALED BY 1/AVSQ, THE LAST 2 BY 1/AV
C AVSQ =       AV*AV
C B =             TRANSPOSE OF THE AUGMENTED REGRESSION MATRIX
C C =             FIRST COMPONENT OF THE PLANE ROTATION USED TO
C             ZERO THE LOWER TRIANGLE OF B**T -- COMPUTED
C             BY SUBROUTINE GIVENS
C DDX,DDY =    LOCAL VARIABLES FOR DX AND DY
C DMIN =       MINIMUM OF THE MAGNITUDES OF THE DIAGONAL
C             ELEMENTS OF THE REGRESSION MATRIX AFTER
C             ZEROS ARE INTRODUCED BELOW THE DIAGONAL
C DTOL =       TOLERANCE FOR DETECTING AN ILL-CONDITIONED
C             SYSTEM.  THE SYSTEM IS ACCEPTED WHEN DMIN
C             .GE. DTOL
C FK =             DATA VALUE AT NODE K -- F(K)
C I =             INDEX FOR A, B, AND NPTS
C IB =             DO-LOOP INDEX FOR BACK SOLVE
C IERR =       ERROR FLAG FOR THE CALL TO STORE2
C IP1 =        I+1
C IRM1 =       IROW-1
C IROW =       ROW INDEX FOR B
C J =             INDEX FOR A AND B
C JP1 =        J+1
C K =             NODAL FUNCTION INDEX AND COLUMN INDEX FOR A
C LMAX =       MAXIMUM NUMBER OF NPTS ELEMENTS (MUST BE CON-
C             SISTENT WITH THE DIMENSION STATEMENT ABOVE)
C LNP =        CURRENT LENGTH OF NPTS
C NEQ =        NUMBER OF EQUATIONS IN THE LEAST SQUARES FIT
C NN,NNQ,NNR = LOCAL COPIES OF N, NQ, AND NR
C NNW =        LOCAL COPY OF NW
C NP =             NPTS ELEMENT
C NPTS =       ARRAY CONTAINING THE INDICES OF A SEQUENCE OF
C             NODES TO BE USED IN THE LEAST SQUARES FIT
C             OR TO COMPUTE RSQ.  THE NODES ARE ORDERED
C             BY DISTANCE FROM K AND THE LAST ELEMENT
C             (USUALLY INDEXED BY LNP) IS USED ONLY TO
C             DETERMINE RQ, OR RSQ(K) IF NW .GT. NQ
C NQWMAX =     MAX(NQ,NW)
C RQ =             RADIUS OF INFLUENCE WHICH ENTERS INTO THE
C             WEIGHTS FOR Q(K) (SEE SUBROUTINE SETUP2)
C RS =             SQUARED DISTANCE BETWEEN K AND NPTS(LNP) --
C             USED TO COMPUTE RQ AND RSQ(K)
C RSMX =       MAXIMUM RSQ ELEMENT ENCOUNTERED
C RSOLD =      SQUARED DISTANCE BETWEEN K AND NPTS(LNP-1) --
C             USED TO COMPUTE A RELATIVE CHANGE IN RS
C             BETWEEN SUCCEEDING NPTS ELEMENTS
C RTOL =       TOLERANCE FOR DETECTING A SUFFICIENTLY LARGE
C             RELATIVE CHANGE IN RS.  IF THE CHANGE IS
C             NOT GREATER THAN RTOL, THE NODES ARE
C             TREATED AS BEING THE SAME DISTANCE FROM K
C RWS =        CURRENT VALUE OF RSQ(K)
C S =             SECOND COMPONENT OF THE PLANE GIVENS ROTATION
C SF =             MARQUARDT STABILIZATION FACTOR USED TO DAMP
C             OUT THE FIRST 3 SOLUTION COMPONENTS (SECOND
C             PARTIALS OF THE QUADRATIC) WHEN THE SYSTEM
C             IS ILL-CONDITIONED.  AS SF INCREASES, THE
C             FITTING FUNCTION APPROACHES A LINEAR
C SUM =        SUM OF SQUARED EUCLIDEAN DISTANCES BETWEEN
C             NODE K AND THE NODES USED IN THE LEAST
C             SQUARES FIT (UNLESS ADDITIONAL NODES ARE
C             ADDED FOR STABILITY)
C T =             TEMPORARY VARIABLE FOR ACCUMULATING A SCALAR
C             PRODUCT IN THE BACK SOLVE
C XK,YK =      COORDINATES OF NODE K -- X(K), Y(K)
C XMN,YMN =    LOCAL VARIABLES FOR XMIN AND YMIN
C
      NN = N
      NNQ = NQ
      NNW = NW
      NNR = NR
      NQWMAX = MAX0(NNQ,NNW)
      LMAX = MIN0(40,NN-1)
      IF (5 .GT. NNQ  .OR. 1 .GT. NNW .OR. NQWMAX .GT.
     &  LMAX .OR. NNR .LT. 1) GO TO 20
C
C CREATE THE CELL DATA STRUCTURE, AND INITIALIZE RSMX.
C
      CALL STORE2 (NN,X,Y,NNR, LCELL,LNEXT,XMN,YMN,DDX,DDY,
     &  IERR)
      IF (IERR .NE. 0) GO TO 22
      RSMX = 0.
C
C OUTER LOOP ON NODE K
C
      DO 18 K = 1,NN
      XK = X(K)
      YK = Y(K)
      FK = F(K)
C
C MARK NODE K TO EXCLUDE IT FROM THE SEARCH FOR NEAREST
C   NEIGHBORS.
C
      LNEXT(K) = -LNEXT(K)
C
C INITIALIZE FOR LOOP ON NPTS.
C
      RS = 0.
      SUM = 0.
      RWS = 0.
      RQ = 0.
      LNP = 0
C
C COMPUTE NPTS, LNP, RWS, NEQ, RQ, AND AVSQ.
C
    1      SUM = SUM + RS
        IF (LNP .EQ. LMAX) GO TO 3
        LNP = LNP + 1
        RSOLD = RS
        CALL GETNP2 (XK,YK,X,Y,NNR,LCELL,LNEXT,XMN,YMN,
     &      DDX,DDY, NP,RS)
        IF (RS .EQ. 0.) GO TO 21
        NPTS(LNP) = NP
        IF ( (RS-RSOLD)/RS .LT. RTOL ) GO TO 1
        IF (RWS .EQ. 0.  .AND.  LNP .GT. NNW) RWS = RS
        IF (RQ .NE. 0. .OR. LNP .LE. NNQ) GO TO 2
C
C   RQ = 0 (NOT YET COMPUTED) AND LNP .GT. NQ.      RQ =
C     SQRT(RS) IS SUFFICIENTLY LARGE TO (STRICTLY) INCLUDE
C     NQ NODES.  THE LEAST SQUARES FIT WILL INCLUDE NEQ =
C     LNP - 1 EQUATIONS FOR 5 .LE. NQ .LE. NEQ .LT. LMAX
C     .LE. N-1.
C
        NEQ = LNP - 1
        RQ = SQRT(RS)
        AVSQ = SUM/FLOAT(NEQ)
C
C   BOTTOM OF LOOP -- TEST FOR TERMINATION.
C
    2        IF (LNP .GT. NQWMAX) GO TO 4
        GO TO 1
C
C ALL LMAX NODES ARE INCLUDED IN NPTS.      RWS AND/OR RQ**2 IS
C   (ARBITRARILY) TAKEN TO BE 10 PERCENT LARGER THAN THE
C   DISTANCE RS TO THE LAST NODE INCLUDED.
C
    3      IF (RWS .EQ. 0.) RWS = 1.1*RS
      IF (RQ .NE. 0.) GO TO 4
        NEQ = LMAX
        RQ = SQRT(1.1*RS)
        AVSQ = SUM/FLOAT(NEQ)
C
C STORE RSQ(K), UPDATE RSMX IF NECESSARY, AND COMPUTE AV.
C
    4      RSQ(K) = RWS
      IF (RWS .GT. RSMX) RSMX = RWS
      AV = SQRT(AVSQ)
C
C SET UP THE AUGMENTED REGRESSION MATRIX (TRANSPOSED) AS THE
C   COLUMNS OF B, AND ZERO OUT THE LOWER TRIANGLE (UPPER
C   TRIANGLE OF B) WITH GIVENS ROTATIONS -- QR DECOMPOSITION
C   WITH ORTHOGONAL MATRIX Q NOT STORED.
C
      I = 0
    5        I = I + 1
        NP = NPTS(I)
        IROW = MIN0(I,6)
        CALL SETUP2 (XK,YK,FK,X(NP),Y(NP),F(NP),AV,AVSQ,
     &      RQ, B(1,IROW))
        IF (I .EQ. 1) GO TO 5
        IRM1 = IROW-1
        DO J = 1,IRM1
          JP1 = J + 1
          CALL GIVENS (B(J,J),B(J,IROW),C,S)
               CALL ROTATE (6-J,C,S,B(JP1,J),B(JP1,IROW))
          end do
        IF (I .LT. NEQ) GO TO 5
C
C TEST THE SYSTEM FOR ILL-CONDITIONING.
C
      DMIN = AMIN1( ABS(B(1,1)),ABS(B(2,2)),ABS(B(3,3)),
     &    ABS(B(4,4)),ABS(B(5,5)) )
      IF (DMIN*RQ .GE. DTOL) GO TO 13
      IF (NEQ .EQ. LMAX) GO TO 10
C
C INCREASE RQ AND ADD ANOTHER EQUATION TO THE SYSTEM TO
C   IMPROVE THE CONDITIONING.  THE NUMBER OF NPTS ELEMENTS
C   IS ALSO INCREASED IF NECESSARY.
C
    7      RSOLD = RS
      NEQ = NEQ + 1
      IF (NEQ .EQ. LMAX) GO TO 9
      IF (NEQ .EQ. LNP) GO TO 8
C
C   NEQ .LT. LNP
C
      NP = NPTS(NEQ+1)
      RS = (X(NP)-XK)**2 + (Y(NP)-YK)**2
      IF ( (RS-RSOLD)/RS .LT. RTOL ) GO TO 7
      RQ = SQRT(RS)
      GO TO 5
C
C   ADD AN ELEMENT TO NPTS.
C
    8      LNP = LNP + 1
      CALL GETNP2 (XK,YK,X,Y,NNR,LCELL,LNEXT,XMN,YMN,
     &    DDX,DDY, NP,RS)
      IF (NP .EQ. 0) GO TO 21
      NPTS(LNP) = NP
      IF ( (RS-RSOLD)/RS .LT. RTOL ) GO TO 7
      RQ = SQRT(RS)
      GO TO 5
C
    9      RQ = SQRT(1.1*RS)
      GO TO 5
C
C STABILIZE THE SYSTEM BY DAMPING SECOND PARTIALS -- ADD
C   MULTIPLES OF THE FIRST THREE UNIT VECTORS TO THE FIRST
C   THREE EQUATIONS.
C
   10      DO I = 1,3
        B(I,6) = SF
        IP1 = I + 1
        DO J = IP1,6
               B(J,6) = 0.
          end do
        DO J = I,5
          JP1 = J + 1
          CALL GIVENS (B(J,J),B(J,6),C,S)
               CALL ROTATE (6-J,C,S,B(JP1,J),B(JP1,6))
          end do
        end do
C
C TEST THE STABILIZED SYSTEM FOR ILL-CONDITIONING.
C
      DMIN = AMIN1( ABS(B(1,1)),ABS(B(2,2)),ABS(B(3,3)),
     &    ABS(B(4,4)),ABS(B(5,5)) )
      IF (DMIN*RQ .LT. DTOL) GO TO 22
C
C SOLVE THE 5 BY 5 TRIANGULAR SYSTEM FOR THE COEFFICIENTS
C
   13      DO IB = 1,5
        I = 6-IB
        T = 0.
        IF (I .EQ. 5) GO TO 15
        IP1 = I + 1
        DO J = IP1,5
               T = T + B(J,I)*A(J,K)
          end do
   15        A(I,K) = (B(6,I)-T)/B(I,I)
        end do
C
C SCALE THE COEFFICIENTS TO ADJUST FOR THE COLUMN SCALING.
C
      DO I = 1,3
             A(I,K) = A(I,K)/AVSQ
        end do
      A(4,K) = A(4,K)/AV
      A(5,K) = A(5,K)/AV
C
C UNMARK K AND THE ELEMENTS OF NPTS.
C
      LNEXT(K) = -LNEXT(K)
      DO I = 1,LNP
        NP = NPTS(I)
             LNEXT(NP) = -LNEXT(NP)
        end do

   18      CONTINUE
C
C NO ERRORS ENCOUNTERED.
C
      XMIN = XMN
      YMIN = YMN
      DX = DDX
      DY = DDY
      RMAX = SQRT(RSMX)
      IER = 0
      RETURN
C
C N, NQ, NW, OR NR IS OUT OF RANGE.
C
   20 IER = 1
      RETURN
C
C DUPLICATE NODES WERE ENCOUNTERED BY GETNP2.
C
   21 IER = 2
      RETURN
C
C NO UNIQUE SOLUTION DUE TO COLLINEAR NODES.
C
   22 XMIN = XMN
      YMIN = YMN
      DX = DDX
      DY = DDY
      IER = 3
      RETURN
      END
      FUNCTION QS2VAL ( PX, PY, N, X, Y, F, NR, LCELL, 
     &  LNEXT, XMIN, YMIN, DX, DY, RMAX, RSQ, A )

      INTEGER N, NR, LCELL(NR,NR), LNEXT(N)
      REAL PX, PY, X(N), Y(N), F(N), XMIN, YMIN, DX, DY,
     &  RMAX, RSQ(N), A(5,N)
C
C***********************************************************
C
C                                    ROBERT RENKA
C                              UNIV. OF NORTH TEXAS
C                                    (817) 565-2767
C                                        10/28/87
C
C   THIS FUNCTION RETURNS THE VALUE Q(PX,PY) WHERE Q IS THE
C WEIGHTED SUM OF QUADRATIC NODAL FUNCTIONS DEFINED IN SUB-
C ROUTINE QSHEP2.  QS2GRD MAY BE CALLED TO COMPUTE A GRADI-
C ENT OF Q ALONG WITH THE VALUE, AND/OR TO TEST FOR ERRORS.
C
C ON INPUT --
C
C      PX,PY = CARTESIAN COORDINATES OF THE POINT P AT
C            WHICH Q IS TO BE EVALUATED.
C
C      N = NUMBER OF NODES AND DATA VALUES DEFINING Q.
C          N .GE. 6.
C
C      X,Y,F = ARRAYS OF LENGTH N CONTAINING THE NODES AND
C            DATA VALUES INTERPOLATED BY Q.
C
C      NR = NUMBER OF ROWS AND COLUMNS IN THE CELL GRID.
C           REFER TO STORE2.  NR .GE. 1.
C
C      LCELL = NR BY NR ARRAY OF NODAL INDICES ASSOCIATED
C            WITH CELLS.  REFER TO STORE2.
C
C      LNEXT = ARRAY OF LENGTH N CONTAINING NEXT-NODE INDI-
C            CES.  REFER TO STORE2.
C
C      XMIN,YMIN,DX,DY = MINIMUM NODAL COORDINATES AND CELL
C                    DIMENSIONS.  DX AND DY MUST BE
C                    POSITIVE.  REFER TO STORE2.
C
C      RMAX = SQUARE ROOT OF THE LARGEST ELEMENT IN RSQ --
C             MAXIMUM RADIUS.
C
C      RSQ = ARRAY OF LENGTH N CONTAINING THE SQUARED RADII
C            WHICH ENTER INTO THE WEIGHTS DEFINING Q.
C
C      A = 5 BY N ARRAY CONTAINING THE COEFFICIENTS FOR THE
C          NODAL FUNCTIONS DEFINING Q.
C
C   INPUT PARAMETERS ARE NOT ALTERED BY THIS FUNCTION.      THE
C PARAMETERS OTHER THAN PX AND PY SHOULD BE INPUT UNALTERED
C FROM THEIR VALUES ON OUTPUT FROM QSHEP2.  THIS FUNCTION
C SHOULD NOT BE CALLED IF A NONZERO ERROR FLAG WAS RETURNED
C BY QSHEP2.
C
C ON OUTPUT --
C
C      QS2VAL = FUNCTION VALUE Q(PX,PY) UNLESS N, NR, DX,
C             DY, OR RMAX IS INVALID, IN WHICH CASE NO
C             VALUE IS RETURNED.
C
C MODULES REQUIRED BY QS2VAL -- NONE
C
C INTRINSIC FUNCTIONS CALLED BY QS2VAL -- IFIX, SQRT
C
C***********************************************************
C
      qs2val = 0.0

      XP = PX
      YP = PY
      IF (N .LT. 6 .OR. NR .LT. 1 .OR. DX .LE. 0. .OR.
     &  DY .LE. 0. .OR. RMAX .LT. 0.) RETURN
C
C SET IMIN, IMAX, JMIN, AND JMAX TO CELL INDICES DEFINING
C   THE RANGE OF THE SEARCH FOR NODES WHOSE RADII INCLUDE
C   P.      THE CELLS WHICH MUST BE SEARCHED ARE THOSE INTER-
C   SECTED BY (OR CONTAINED IN) A CIRCLE OF RADIUS RMAX
C   CENTERED AT P.
C
      IMIN = IFIX((XP-XMIN-RMAX)/DX) + 1
      IMAX = IFIX((XP-XMIN+RMAX)/DX) + 1
      IF (IMIN .LT. 1) IMIN = 1
      IF (IMAX .GT. NR) IMAX = NR
      JMIN = IFIX((YP-YMIN-RMAX)/DY) + 1
      JMAX = IFIX((YP-YMIN+RMAX)/DY) + 1
      IF (JMIN .LT. 1) JMIN = 1
      IF (JMAX .GT. NR) JMAX = NR
C
C THE FOLLOWING IS A TEST FOR NO CELLS WITHIN THE CIRCLE
C   OF RADIUS RMAX.
C
      IF (IMIN .GT. IMAX .OR. JMIN .GT. JMAX) GO TO 5
C
C ACCUMULATE WEIGHT VALUES IN SW AND WEIGHTED NODAL FUNCTION
C   VALUES IN SWQ.  THE WEIGHTS ARE W(K) = ((R-D)+/(R*D))**2
C   FOR R**2 = RSQ(K) AND D = DISTANCE BETWEEN P AND NODE K.
C
      SW = 0.
      SWQ = 0.
C
C OUTER LOOP ON CELLS (I,J).
C
      DO J = JMIN,JMAX
      DO I = IMIN,IMAX
        K = LCELL(I,J)
        IF (K .EQ. 0) GO TO 3
C
C INNER LOOP ON NODES K.
C
    1        DELX = XP - X(K)
        DELY = YP - Y(K)
        DXSQ = DELX*DELX
        DYSQ = DELY*DELY
        DS = DXSQ + DYSQ
        RS = RSQ(K)
        IF (DS .GE. RS) GO TO 2
        IF (DS .EQ. 0.) GO TO 4
        RDS = RS*DS
        RD = SQRT(RDS)
        W = (RS+DS-RD-RD)/RDS
        SW = SW + W
        SWQ = SWQ + W*(A(1,K)*DXSQ + A(2,K)*DELX*DELY +
     &      A(3,K)*DYSQ + A(4,K)*DELX +
     &      A(5,K)*DELY + F(K))
C
C BOTTOM OF LOOP ON NODES IN CELL (I,J).
C
    2        KP = K
        K = LNEXT(KP)
        IF (K .NE. KP) GO TO 1
    3        CONTINUE
        end do
      end do
C
C SW = 0 IFF P IS NOT WITHIN THE RADIUS R(K) FOR ANY NODE K.
C
      IF (SW .EQ. 0.) GO TO 5
      QS2VAL = SWQ/SW
      RETURN
C
C (PX,PY) = (X(K),Y(K))
C
    4 QS2VAL = F(K)
      RETURN
C
C ALL WEIGHTS ARE 0 AT P.
C
    5 QS2VAL = 0.
      RETURN
      END
      SUBROUTINE QS2GRD (PX,PY,N,X,Y,F,NR,LCELL,LNEXT,XMIN,
     &  YMIN,DX,DY,RMAX,RSQ,A, Q,QX,QY,IER)
      INTEGER N, NR, LCELL(NR,NR), LNEXT(N), IER
      REAL PX, PY, X(N), Y(N), F(N), XMIN, YMIN, DX, DY,
     &  RMAX, RSQ(N), A(5,N), Q, QX, QY
C
C***********************************************************
C
C                                    ROBERT RENKA
C                              UNIV. OF NORTH TEXAS
C                                    (817) 565-2767
C                                        10/28/87
C
C   THIS SUBROUTINE COMPUTES THE VALUE AND GRADIENT AT
C (PX,PY) OF THE INTERPOLATORY FUNCTION Q DEFINED IN SUB-
C ROUTINE QSHEP2.  Q(X,Y) IS A WEIGHTED SUM OF QUADRATIC
C NODAL FUNCTIONS.
C
C ON INPUT --
C
C      PX,PY = CARTESIAN COORDINATES OF THE POINT AT WHICH
C            Q AND ITS PARTIALS ARE TO BE EVALUATED.
C
C      N = NUMBER OF NODES AND DATA VALUES DEFINING Q.
C          N .GE. 6.
C
C      X,Y,F = ARRAYS OF LENGTH N CONTAINING THE NODES AND
C            DATA VALUES INTERPOLATED BY Q.
C
C      NR = NUMBER OF ROWS AND COLUMNS IN THE CELL GRID.
C           REFER TO STORE2.  NR .GE. 1.
C
C      LCELL = NR BY NR ARRAY OF NODAL INDICES ASSOCIATED
C            WITH CELLS.  REFER TO STORE2.
C
C      LNEXT = ARRAY OF LENGTH N CONTAINING NEXT-NODE INDI-
C            CES.  REFER TO STORE2.
C
C      XMIN,YMIN,DX,DY = MINIMUM NODAL COORDINATES AND CELL
C                    DIMENSIONS.  DX AND DY MUST BE
C                    POSITIVE.  REFER TO STORE2.
C
C      RMAX = SQUARE ROOT OF THE LARGEST ELEMENT IN RSQ --
C             MAXIMUM RADIUS.
C
C      RSQ = ARRAY OF LENGTH N CONTAINING THE SQUARED RADII
C            WHICH ENTER INTO THE WEIGHTS DEFINING Q.
C
C      A = 5 BY N ARRAY CONTAINING THE COEFFICIENTS FOR THE
C          NODAL FUNCTIONS DEFINING Q.
C
C   INPUT PARAMETERS ARE NOT ALTERED BY THIS SUBROUTINE.
C THE PARAMETERS OTHER THAN PX AND PY SHOULD BE INPUT UNAL-
C TERED FROM THEIR VALUES ON OUTPUT FROM QSHEP2.  THIS SUB-
C ROUTINE SHOULD NOT BE CALLED IF A NONZERO ERROR FLAG WAS
C RETURNED BY QSHEP2.
C
C ON OUTPUT --
C
C      Q = VALUE OF Q AT (PX,PY) UNLESS IER .EQ. 1, IN
C          WHICH CASE NO VALUES ARE RETURNED.
C
C      QX,QY = FIRST PARTIAL DERIVATIVES OF Q AT (PX,PY)
C            UNLESS IER .EQ. 1.
C
C      IER = ERROR INDICATOR
C            IER = 0 IF NO ERRORS WERE ENCOUNTERED.
C            IER = 1 IF N, NR, DX, DY OR RMAX IS INVALID.
C            IER = 2 IF NO ERRORS WERE ENCOUNTERED BUT
C                  (PX,PY) IS NOT WITHIN THE RADIUS R(K)
C                  FOR ANY NODE K (AND THUS Q=QX=QY=0).
C
C MODULES REQUIRED BY QS2GRD -- NONE
C
C INTRINSIC FUNCTIONS CALLED BY QS2GRD -- IFIX, SQRT
C
C***********************************************************
C
      XP = PX
      YP = PY
      IF (N .LT. 6  .OR. NR .LT. 1 .OR. DX .LE. 0. .OR.
     &  DY .LE. 0. .OR. RMAX .LT. 0.) GO TO 5
C
C SET IMIN, IMAX, JMIN, AND JMAX TO CELL INDICES DEFINING
C   THE RANGE OF THE SEARCH FOR NODES WHOSE RADII INCLUDE
C   P.      THE CELLS WHICH MUST BE SEARCHED ARE THOSE INTER-
C   SECTED BY (OR CONTAINED IN) A CIRCLE OF RADIUS RMAX
C   CENTERED AT P.
C
      IMIN = IFIX((XP-XMIN-RMAX)/DX) + 1
      IMAX = IFIX((XP-XMIN+RMAX)/DX) + 1
      IF (IMIN .LT. 1) IMIN = 1
      IF (IMAX .GT. NR) IMAX = NR
      JMIN = IFIX((YP-YMIN-RMAX)/DY) + 1
      JMAX = IFIX((YP-YMIN+RMAX)/DY) + 1
      IF (JMIN .LT. 1) JMIN = 1
      IF (JMAX .GT. NR) JMAX = NR
C
C THE FOLLOWING IS A TEST FOR NO CELLS WITHIN THE CIRCLE
C   OF RADIUS RMAX.
C
      IF (IMIN .GT. IMAX .OR. JMIN .GT. JMAX) GO TO 6
C
C Q = SWQ/SW = SUM(W(K)*Q(K))/SUM(W(K)) WHERE THE SUM IS
C   FROM K = 1 TO N, Q(K) IS THE QUADRATIC NODAL FUNCTION,
C   AND W(K) = ((R-D)+/(R*D))**2 FOR RADIUS R(K) AND DIST-
C   ANCE D(K).      THUS
C
C       QX = (SWQX*SW - SWQ*SWX)/SW**2  AND
C       QY = (SWQY*SW - SWQ*SWY)/SW**2
C
C   WHERE SWQX AND SWX ARE PARTIAL DERIVATIVES WITH RESPECT
C   TO X OF SWQ AND SW, RESPECTIVELY.  SWQY AND SWY ARE DE-
C   FINED SIMILARLY.
C
      SW = 0.
      SWX = 0.
      SWY = 0.
      SWQ = 0.
      SWQX = 0.
      SWQY = 0.
C
C OUTER LOOP ON CELLS (I,J).
C
      DO J = JMIN,JMAX
      DO I = IMIN,IMAX
        K = LCELL(I,J)
        IF (K .EQ. 0) GO TO 3
C
C INNER LOOP ON NODES K.
C
    1        DELX = XP - X(K)
        DELY = YP - Y(K)
        DXSQ = DELX*DELX
        DYSQ = DELY*DELY
        DS = DXSQ + DYSQ
        RS = RSQ(K)
        IF (DS .GE. RS) GO TO 2
        IF (DS .EQ. 0.) GO TO 4
        RDS = RS*DS
        RD = SQRT(RDS)
        W = (RS+DS-RD-RD)/RDS
        T = 2.*(RD-RS)/(DS*RDS)
        WX = DELX*T
        WY = DELY*T
        QKX = 2.*A(1,K)*DELX + A(2,K)*DELY
        QKY = A(2,K)*DELX + 2.*A(3,K)*DELY
        QK = (QKX*DELX + QKY*DELY)/2.
        QKX = QKX + A(4,K)
        QKY = QKY + A(5,K)
        QK = QK + A(4,K)*DELX + A(5,K)*DELY + F(K)
        SW = SW + W
        SWX = SWX + WX
        SWY = SWY + WY
        SWQ = SWQ + W*QK
        SWQX = SWQX + WX*QK + W*QKX
        SWQY = SWQY + WY*QK + W*QKY
C
C BOTTOM OF LOOP ON NODES IN CELL (I,J).
C
    2        KP = K
        K = LNEXT(KP)
        IF (K .NE. KP) GO TO 1
    3        CONTINUE
        end do
      end do
C
C SW = 0 IFF P IS NOT WITHIN THE RADIUS R(K) FOR ANY NODE K.
C
      IF (SW .EQ. 0.) GO TO 6
      Q = SWQ/SW
      SWS = SW*SW
      QX = (SWQX*SW - SWQ*SWX)/SWS
      QY = (SWQY*SW - SWQ*SWY)/SWS
      IER = 0
      RETURN
C
C (PX,PY) = (X(K),Y(K))
C
    4 Q = F(K)
      QX = A(4,K)
      QY = A(5,K)
      IER = 0
      RETURN
C
C INVALID INPUT PARAMETER.
C
    5 IER = 1
      RETURN
C
C NO CELLS CONTAIN A POINT WITHIN RMAX OF P, OR
C   SW = 0 AND THUS DS .GE. RSQ(K) FOR ALL K.
C
    6 Q = 0.
      QX = 0.
      QY = 0.
      IER = 2
      RETURN
      END
      SUBROUTINE GETNP2 (PX,PY,X,Y,NR,LCELL,LNEXT,XMIN,YMIN,
     &  DX,DY, NP,DSQ)
      INTEGER NR, LCELL(NR,NR), LNEXT(1), NP
      REAL    PX, PY, X(1), Y(1), XMIN, YMIN, DX, DY, DSQ
C
C***********************************************************
C
C                                    ROBERT RENKA
C                              UNIV. OF NORTH TEXAS
C                                    (817) 565-2767
C
C   GIVEN A SET OF N NODES AND THE DATA STRUCTURE DEFINED IN
C SUBROUTINE STORE2, THIS SUBROUTINE USES THE CELL METHOD TO
C FIND THE CLOSEST UNMARKED NODE NP TO A SPECIFIED POINT P.
C NP IS THEN MARKED BY SETTING LNEXT(NP) TO -LNEXT(NP).  (A
C NODE IS MARKED IF AND ONLY IF THE CORRESPONDING LNEXT ELE-
C MENT IS NEGATIVE.  THE ABSOLUTE VALUES OF LNEXT ELEMENTS,
C HOWEVER, MUST BE PRESERVED.)      THUS, THE CLOSEST M NODES TO
C P MAY BE DETERMINED BY A SEQUENCE OF M CALLS TO THIS ROU-
C TINE.  NOTE THAT IF THE NEAREST NEIGHBOR TO NODE K IS TO
C BE DETERMINED (PX = X(K) AND PY = Y(K)), THEN K SHOULD BE
C MARKED BEFORE THE CALL TO THIS ROUTINE.
C   THE SEARCH IS BEGUN IN THE CELL CONTAINING (OR CLOSEST
C TO) P AND PROCEEDS OUTWARD IN RECTANGULAR LAYERS UNTIL ALL
C CELLS WHICH CONTAIN POINTS WITHIN DISTANCE R OF P HAVE
C BEEN SEARCHED, WHERE R IS THE DISTANCE FROM P TO THE FIRST
C UNMARKED NODE ENCOUNTERED (INFINITE IF NO UNMARKED NODES
C ARE PRESENT).
C
C ON INPUT --
C
C      PX,PY = CARTESIAN COORDINATES OF THE POINT P WHOSE
C            NEAREST UNMARKED NEIGHBOR IS TO BE FOUND.
C
C      X,Y = ARRAYS OF LENGTH N, FOR N .GE. 2, CONTAINING
C            THE CARTESIAN COORDINATES OF THE NODES.
C
C      NR = NUMBER OF ROWS AND COLUMNS IN THE CELL GRID.
C           NR .GE. 1.
C
C      LCELL = NR BY NR ARRAY OF NODAL INDICES ASSOCIATED
C            WITH CELLS.
C
C      LNEXT = ARRAY OF LENGTH N CONTAINING NEXT-NODE INDI-
C            CES (OR THEIR NEGATIVES).
C
C      XMIN,YMIN,DX,DY = MINIMUM NODAL COORDINATES AND CELL
C                    DIMENSIONS.  DX AND DY MUST BE
C                    POSITIVE.
C
C   INPUT PARAMETERS OTHER THAN LNEXT ARE NOT ALTERED BY
C THIS ROUTINE.  WITH THE EXCEPTION OF (PX,PY) AND THE SIGNS
C OF LNEXT ELEMENTS, THESE PARAMETERS SHOULD BE UNALTERED
C FROM THEIR VALUES ON OUTPUT FROM SUBROUTINE STORE2.
C
C ON OUTPUT --
C
C      NP = INDEX (FOR X AND Y) OF THE NEAREST UNMARKED
C           NODE TO P, OR 0 IF ALL NODES ARE MARKED OR NR
C           .LT. 1 OR DX .LE. 0 OR DY .LE. 0.      LNEXT(NP)
C           .LT. 0 IF NP .NE. 0.
C
C      DSQ = SQUARED EUCLIDEAN DISTANCE BETWEEN P AND NODE
C            NP, OR 0 IF NP = 0.
C
C MODULES REQUIRED BY GETNP2 -- NONE
C
C INTRINSIC FUNCTIONS CALLED BY GETNP2 -- IABS, IFIX, SQRT
C
C***********************************************************
C
      LOGICAL FIRST
      XP = PX
      YP = PY
C
C TEST FOR INVALID INPUT PARAMETERS.
C
      IF (NR .LT. 1 .OR. DX .LE. 0. .OR. DY .LE. 0.)
     .  GO TO 9
C
C INITIALIZE PARAMETERS --
C
C   FIRST = TRUE IFF THE FIRST UNMARKED NODE HAS YET TO BE
C          ENCOUNTERED,
C   IMIN,IMAX,JMIN,JMAX = CELL INDICES DEFINING THE RANGE OF
C                    THE SEARCH,
C   DELX,DELY = PX-XMIN AND PY-YMIN,
C   I0,J0 = CELL CONTAINING OR CLOSEST TO P,
C   I1,I2,J1,J2 = CELL INDICES OF THE LAYER WHOSE INTERSEC-
C              TION WITH THE RANGE DEFINED BY IMIN,...,
C              JMAX IS CURRENTLY BEING SEARCHED.
C
      FIRST = .TRUE.
      IMIN = 1
      IMAX = NR
      JMIN = 1
      JMAX = NR
      DELX = XP - XMIN
      DELY = YP - YMIN
      I0 = IFIX(DELX/DX) + 1
      IF (I0 .LT. 1) I0 = 1
      IF (I0 .GT. NR) I0 = NR
      J0 = IFIX(DELY/DY) + 1
      IF (J0 .LT. 1) J0 = 1
      IF (J0 .GT. NR) J0 = NR
      I1 = I0
      I2 = I0
      J1 = J0
      J2 = J0
C
C OUTER LOOP ON LAYERS, INNER LOOP ON LAYER CELLS, EXCLUDING
C   THOSE OUTSIDE THE RANGE (IMIN,IMAX) X (JMIN,JMAX).
C
    1 continue

      DO 6 J = J1,J2
      IF (J .GT. JMAX) GO TO 7
      IF (J .LT. JMIN) GO TO 6
      DO 5 I = I1,I2
        IF (I .GT. IMAX) GO TO 6
        IF (I .LT. IMIN) GO TO 5
        IF (J .NE. J1 .AND. J .NE. J2  .AND.  I .NE. I1
     &      .AND. I .NE. I2) GO TO 5
C
C SEARCH CELL (I,J) FOR UNMARKED NODES L.
C
        L = LCELL(I,J)
        IF (L .EQ. 0) GO TO 5
C
C   LOOP ON NODES IN CELL (I,J).
C
    2        LN = LNEXT(L)
        IF (LN .LT. 0) GO TO 4
C
C   NODE L IS NOT MARKED.
C
        RSQ = (X(L)-XP)**2 + (Y(L)-YP)**2
        IF (.NOT. FIRST) GO TO 3
C
C   NODE L IS THE FIRST UNMARKED NEIGHBOR OF P ENCOUNTERED.
C     INITIALIZE LMIN TO THE CURRENT CANDIDATE FOR NP, AND
C     RSMIN TO THE SQUARED DISTANCE FROM P TO LMIN.  IMIN,
C     IMAX, JMIN, AND JMAX ARE UPDATED TO DEFINE THE SMAL-
C     LEST RECTANGLE CONTAINING A CIRCLE OF RADIUS R =
C     SQRT(RSMIN) CENTERED AT P, AND CONTAINED IN (1,NR) X
C     (1,NR) (EXCEPT THAT, IF P IS OUTSIDE THE RECTANGLE
C     DEFINED BY THE NODES, IT IS POSSIBLE THAT IMIN .GT.
C     NR, IMAX .LT. 1, JMIN .GT. NR, OR JMAX .LT. 1).  FIRST
C     IS RESET TO FALSE.
C
        LMIN = L
        RSMIN = RSQ
        R = SQRT(RSMIN)
        IMIN = IFIX((DELX-R)/DX) + 1
        IF (IMIN .LT. 1) IMIN = 1
        IMAX = IFIX((DELX+R)/DX) + 1
        IF (IMAX .GT. NR) IMAX = NR
        JMIN = IFIX((DELY-R)/DY) + 1
        IF (JMIN .LT. 1) JMIN = 1
        JMAX = IFIX((DELY+R)/DY) + 1
        IF (JMAX .GT. NR) JMAX = NR
        FIRST = .FALSE.
        GO TO 4
C
C   TEST FOR NODE L CLOSER THAN LMIN TO P.
C
    3        IF (RSQ .GE. RSMIN) GO TO 4
C
C   UPDATE LMIN AND RSMIN.
C
        LMIN = L
        RSMIN = RSQ
C
C   TEST FOR TERMINATION OF LOOP ON NODES IN CELL (I,J).
C
    4        IF (IABS(LN) .EQ. L) GO TO 5
        L = IABS(LN)
        GO TO 2
    5        CONTINUE
    6      CONTINUE
C
C TEST FOR TERMINATION OF LOOP ON CELL LAYERS.
C
    7 IF (I1 .LE. IMIN .AND. I2 .GE. IMAX  .AND.
     .  J1 .LE. JMIN .AND. J2 .GE. JMAX) GO TO 8
      I1 = I1 - 1
      I2 = I2 + 1
      J1 = J1 - 1
      J2 = J2 + 1
      GO TO 1
C
C UNLESS NO UNMARKED NODES WERE ENCOUNTERED, LMIN IS THE
C   CLOSEST UNMARKED NODE TO P.
C
    8 IF (FIRST) GO TO 9
      NP = LMIN
      DSQ = RSMIN
      LNEXT(LMIN) = -LNEXT(LMIN)
      RETURN
C
C ERROR -- NR, DX, OR DY IS INVALID OR ALL NODES ARE MARKED.
C
    9 NP = 0
      DSQ = 0.
      RETURN
      END
      SUBROUTINE GIVENS ( A,B, C,S)
      REAL A, B, C, S
C
C***********************************************************
C
C                                    ROBERT RENKA
C                              UNIV. OF NORTH TEXAS
C                                    (817) 565-2767
C
C   THIS ROUTINE CONSTRUCTS THE GIVENS PLANE ROTATION --
C     ( C  S)
C G = (     ) WHERE C*C + S*S = 1 -- WHICH ZEROS THE SECOND
C     (-S  C)
C ENTRY OF THE 2-VECTOR (A B)-TRANSPOSE.  A CALL TO GIVENS
C IS NORMALLY FOLLOWED BY A CALL TO ROTATE WHICH APPLIES
C THE TRANSFORMATION TO A 2 BY N MATRIX.  THIS ROUTINE WAS
C TAKEN FROM LINPACK.
C
C ON INPUT --
C
C      A,B = COMPONENTS OF THE 2-VECTOR TO BE ROTATED.
C
C ON OUTPUT --
C
C      A = VALUE OVERWRITTEN BY R = +/-SQRT(A*A + B*B)
C
C      B = VALUE OVERWRITTEN BY A VALUE Z WHICH ALLOWS C
C          AND S TO BE RECOVERED AS FOLLOWS --
C            C = SQRT(1-Z*Z), S=Z     IF ABS(Z) .LE. 1.
C            C = 1/Z, S = SQRT(1-C*C) IF ABS(Z) .GT. 1.
C
C      C = +/-(A/R)
C
C      S = +/-(B/R)
C
C MODULES REQUIRED BY GIVENS -- NONE
C
C INTRINSIC FUNCTIONS CALLED BY GIVENS - ABS, SQRT
C
C***********************************************************
C
      REAL AA, BB, R, U, V
C
C LOCAL PARAMETERS --
C
C AA,BB = LOCAL COPIES OF A AND B
C R =        C*A + S*B = +/-SQRT(A*A+B*B)
C U,V =   VARIABLES USED TO SCALE A AND B FOR COMPUTING R
C
      AA = A
      BB = B
      IF (ABS(AA) .LE. ABS(BB)) GO TO 1
C
C ABS(A) .GT. ABS(B)
C
      U = AA + AA
      V = BB/U
      R = SQRT(.25 + V*V) * U
      C = AA/R
      S = V * (C + C)
C
C NOTE THAT R HAS THE SIGN OF A, C .GT. 0, AND S HAS
C   SIGN(A)*SIGN(B).
C
      B = S
      A = R
      RETURN
C
C ABS(A) .LE. ABS(B)
C
    1 IF (BB .EQ. 0.) GO TO 2
      U = BB + BB
      V = AA/U
C
C STORE R IN A.
C
      A = SQRT(.25 + V*V) * U
      S = BB/A
      C = V * (S + S)
C
C NOTE THAT R HAS THE SIGN OF B, S .GT. 0, AND C HAS
C   SIGN(A)*SIGN(B).
C
      B = 1.
      IF (C .NE. 0.) B = 1./C
      RETURN
C
C A = B = 0.
C
    2 C = 1.
      S = 0.
      RETURN
      END
      SUBROUTINE ROTATE (N,C,S, X,Y )
      INTEGER N
      REAL    C, S, X(N), Y(N)
C
C***********************************************************
C
C                                    ROBERT RENKA
C                              UNIV. OF NORTH TEXAS
C                                    (817) 565-2767
C
C                                   ( C  S)
C   THIS ROUTINE APPLIES THE GIVENS ROTATION (         ) TO THE
C                                   (-S  C)
C            (X(1) ... X(N))
C 2 BY N MATRIX (            ).
C            (Y(1) ... Y(N))
C
C ON INPUT --
C
C      N = NUMBER OF COLUMNS TO BE ROTATED.
C
C      C,S = ELEMENTS OF THE GIVENS ROTATION.      THESE MAY BE
C            DETERMINED BY SUBROUTINE GIVENS.
C
C      X,Y = ARRAYS OF LENGTH .GE. N CONTAINING THE VECTORS
C            TO BE ROTATED.
C
C PARAMETERS N, C, AND S ARE NOT ALTERED BY THIS ROUTINE.
C
C ON OUTPUT --
C
C      X,Y = ROTATED VECTORS.
C
C MODULES REQUIRED BY ROTATE -- NONE
C
C***********************************************************
C
      INTEGER I
      REAL    XI, YI
C
C LOCAL PARAMETERS --
C
C I =        DO-LOOP INDEX
C XI,YI = X(I), Y(I)
C
      IF (N .LE. 0 .OR. (C .EQ. 1. .AND. S .EQ. 0.)) RETURN
      DO I = 1,N
      XI = X(I)
      YI = Y(I)
      X(I) = C*XI + S*YI
      Y(I) = -S*XI + C*YI
      end do

      RETURN
      END
      SUBROUTINE SETUP2 (XK,YK,FK,XI,YI,FI,S1,S2,R, ROW)
      REAL XK, YK, FK, XI, YI, FI, S1, S2, R, ROW(6)
C
C***********************************************************
C
C                                    ROBERT RENKA
C                              UNIV. OF NORTH TEXAS
C                                    (817) 565-2767
C
C   THIS ROUTINE SETS UP THE I-TH ROW OF AN AUGMENTED RE-
C GRESSION MATRIX FOR A WEIGHTED LEAST-SQUARES FIT OF A
C QUADRATIC FUNCTION Q(X,Y) TO A SET OF DATA VALUES F, WHERE
C Q(XK,YK) = FK.  THE FIRST 3 COLUMNS (QUADRATIC TERMS) ARE
C SCALED BY 1/S2 AND THE FOURTH AND FIFTH COLUMNS (LINEAR
C TERMS) ARE SCALED BY 1/S1.  THE WEIGHT IS (R-D)/(R*D) IF
C R .GT. D AND 0 IF R .LE. D, WHERE D IS THE DISTANCE
C BETWEEN NODES I AND K.
C
C ON INPUT --
C
C      XK,YK,FK = COORDINATES AND DATA VALUE AT NODE K --
C               INTERPOLATED BY Q.
C
C      XI,YI,FI = COORDINATES AND DATA VALUE AT NODE I.
C
C      S1,S2 = RECIPROCALS OF THE SCALE FACTORS.
C
C      R = RADIUS OF INFLUENCE ABOUT NODE K DEFINING THE
C          WEIGHT.
C
C      ROW = ARRAY OF LENGTH 6.
C
C INPUT PARAMETERS ARE NOT ALTERED BY THIS ROUTINE.
C
C ON OUTPUT --
C
C      ROW = VECTOR CONTAINING A ROW OF THE AUGMENTED
C            REGRESSION MATRIX.
C
C MODULES REQUIRED BY SETUP2 -- NONE
C
C INTRINSIC FUNCTION CALLED BY SETUP2 -- SQRT
C
C***********************************************************
C
      INTEGER I
      REAL    DX, DY, DXSQ, DYSQ, D, W, W1, W2
C
C LOCAL PARAMETERS -
C
C I =       DO-LOOP INDEX
C DX =       XI - XK
C DY =       YI - YK
C DXSQ = DX*DX
C DYSQ = DY*DY
C D =       DISTANCE BETWEEN NODES K AND I
C W =       WEIGHT ASSOCIATED WITH THE ROW
C W1 =       W/S1
C W2 =       W/S2
C
      DX = XI - XK
      DY = YI - YK
      DXSQ = DX*DX
      DYSQ = DY*DY
      D = SQRT(DXSQ + DYSQ)
      IF (D .LE. 0.  .OR.  D .GE. R) GO TO 1
      W = (R-D)/R/D
      W1 = W/S1
      W2 = W/S2
      ROW(1) = DXSQ*W2
      ROW(2) = DX*DY*W2
      ROW(3) = DYSQ*W2
      ROW(4) = DX*W1
      ROW(5) = DY*W1
      ROW(6) = (FI - FK)*W
      RETURN
C
C NODES K AND I COINCIDE OR NODE I IS OUTSIDE OF THE RADIUS
C   OF INFLUENCE.  SET ROW TO THE ZERO VECTOR.
C
    1 DO I = 1,6
           ROW(I) = 0.0
      end do

      RETURN
      END
      SUBROUTINE STORE2 (N,X,Y,NR, LCELL,LNEXT,XMIN,YMIN,DX,
     &  DY,IER)
      INTEGER N, NR, LCELL(NR,NR), LNEXT(N), IER
      REAL    X(N), Y(N), XMIN, YMIN, DX, DY
C
C***********************************************************
C
C                                    ROBERT RENKA
C                              UNIV. OF NORTH TEXAS
C                                    (817) 565-2767
C
C   GIVEN A SET OF N ARBITRARILY DISTRIBUTED NODES IN THE
C PLANE, THIS SUBROUTINE CREATES A DATA STRUCTURE FOR A
C CELL-BASED METHOD OF SOLVING CLOSEST-POINT PROBLEMS.      THE
C SMALLEST RECTANGLE CONTAINING THE NODES IS PARTITIONED
C INTO AN NR BY NR UNIFORM GRID OF CELLS, AND NODES ARE AS-
C SOCIATED WITH CELLS.      IN PARTICULAR, THE DATA STRUCTURE
C STORES THE INDICES OF THE NODES CONTAINED IN EACH CELL.
C FOR A UNIFORM RANDOM DISTRIBUTION OF NODES, THE NEAREST
C NODE TO AN ARBITRARY POINT CAN BE DETERMINED IN CONSTANT
C EXPECTED TIME.
C
C ON INPUT --
C
C      N = NUMBER OF NODES.  N .GE. 2.
C
C      X,Y = ARRAYS OF LENGTH N CONTAINING THE CARTESIAN
C            COORDINATES OF THE NODES.
C
C      NR = NUMBER OF ROWS AND COLUMNS IN THE GRID.  THE
C           CELL DENSITY (AVERAGE NUMBER OF NODES PER CELL)
C           IS D = N/(NR**2).      A RECOMMENDED VALUE, BASED
C           ON EMPIRICAL EVIDENCE, IS D = 3 -- NR =
C           SQRT(N/3).  NR .GE. 1.
C
C THE ABOVE PARAMETERS ARE NOT ALTERED BY THIS ROUTINE.
C
C      LCELL = ARRAY OF LENGTH .GE. NR**2.
C
C      LNEXT = ARRAY OF LENGTH .GE. N.
C
C ON OUTPUT --
C
C      LCELL = NR BY NR CELL ARRAY SUCH THAT LCELL(I,J)
C            CONTAINS THE INDEX (FOR X AND Y) OF THE
C            FIRST NODE (NODE WITH SMALLEST INDEX) IN
C            CELL (I,J), OR LCELL(I,J) = 0 IF NO NODES
C            ARE CONTAINED IN THE CELL.  THE UPPER RIGHT
C            CORNER OF CELL (I,J) HAS COORDINATES (XMIN+
C            I*DX,YMIN+J*DY).  LCELL IS NOT DEFINED IF
C            IER .NE. 0.
C
C      LNEXT = ARRAY OF NEXT-NODE INDICES SUCH THAT
C            LNEXT(K) CONTAINS THE INDEX OF THE NEXT NODE
C            IN THE CELL WHICH CONTAINS NODE K, OR
C            LNEXT(K) = K IF K IS THE LAST NODE IN THE
C            CELL FOR K = 1,...,N.  (THE NODES CONTAINED
C            IN A CELL ARE ORDERED BY THEIR INDICES.)
C            IF, FOR EXAMPLE, CELL (I,J) CONTAINS NODES
C            2, 3, AND 5 (AND NO OTHERS), THEN LCELL(I,J)
C            = 2, LNEXT(2) = 3, LNEXT(3) = 5, AND
C            LNEXT(5) = 5.  LNEXT IS NOT DEFINED IF
C            IER .NE. 0.
C
C      XMIN,YMIN = CARTESIAN COORDINATES OF THE LOWER LEFT
c       CORNER OF THE RECTANGLE DEFINED BY THE
C       NODES (SMALLEST NODAL COORDINATES) UN-
C       LESS IER = 1.  THE UPPER RIGHT CORNER IS
C       (XMAX,YMAX) FOR XMAX = XMIN + NR*DX AND
C       YMAX = YMIN + NR*DY.
C
C      DX,DY = DIMENSIONS OF THE CELLS UNLESS IER = 1.  DX
C       = (XMAX-XMIN)/NR AND DY = (YMAX-YMIN)/NR
C       WHERE XMIN, XMAX, YMIN, AND YMAX ARE THE
C       EXTREMA OF X AND Y.
C
C      IER = ERROR INDICATOR --
C            IER = 0 IF NO ERRORS WERE ENCOUNTERED.
C            IER = 1 IF N .LT. 2 OR NR .LT. 1.
C            IER = 2 IF DX = 0 OR DY = 0.
C
C MODULES REQUIRED BY STORE2 -- NONE
C
C INTRINSIC FUNCTIONS CALLED BY STORE2 -- FLOAT, IFIX
C
C***********************************************************
C
      NN = N
      NNR = NR
      IF (NN .LT. 2  .OR.  NNR .LT. 1) GO TO 4
C
C COMPUTE THE DIMENSIONS OF THE RECTANGLE CONTAINING THE
C   NODES.
C
      XMN = X(1)
      XMX = XMN
      YMN = Y(1)
      YMX = YMN
      DO K = 2,NN
      IF (X(K) .LT. XMN) XMN = X(K)
      IF (X(K) .GT. XMX) XMX = X(K)
      IF (Y(K) .LT. YMN) YMN = Y(K)
           IF (Y(K) .GT. YMX) YMX = Y(K)
      end do
      XMIN = XMN
      YMIN = YMN
C
C COMPUTE CELL DIMENSIONS AND TEST FOR ZERO AREA.
C
      DELX = (XMX-XMN)/FLOAT(NNR)
      DELY = (YMX-YMN)/FLOAT(NNR)
      DX = DELX
      DY = DELY
      IF (DELX .EQ. 0.  .OR. DELY .EQ. 0.) GO TO 5
C
C INITIALIZE LCELL.
C
      DO J = 1,NNR
      DO I = 1,NNR
             LCELL(I,J) = 0
        end do
      end do
C
C LOOP ON NODES, STORING INDICES IN LCELL AND LNEXT.
C
      NP1 = NN + 1
      DO K = 1,NN
      KB = NP1 - K
      I = IFIX((X(KB)-XMN)/DELX) + 1
      IF (I .GT. NNR) I = NNR
      J = IFIX((Y(KB)-YMN)/DELY) + 1
      IF (J .GT. NNR) J = NNR
      L = LCELL(I,J)
      LNEXT(KB) = L
      IF (L .EQ. 0) LNEXT(KB) = KB
           LCELL(I,J) = KB
      end do
C
C NO ERRORS ENCOUNTERED
C
      IER = 0
      RETURN
C
C INVALID INPUT PARAMETER
C
    4 IER = 1
      RETURN
C
C DX = 0 OR DY = 0
C
    5 IER = 2
      RETURN
      END
