C      ALGORITHM 660, COLLECTED ALGORITHMS FROM ACM.
C      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
C      VOL. 14, NO. 2, P.149.
C
C			   QS2TEST
C
C   THIS PROGRAM TESTS THE SCATTERED DATA INTERPOLATION
C PACKAGE QSHEP2D BY PRINTING THE MAXIMUM ERRORS ASSOCIATED
C WITH INTERPOLATED VALUES AND GRADIENTS ON A 10 BY 10
C UNIFORM GRID IN THE UNIT SQUARE.  THE DATA SET CONSISTS
C OF 36 NODES WITH DATA VALUES TAKEN FROM A QUADRATIC FUNC-
C TION FOR WHICH THE METHOD IS EXACT.  THE RATIO OF MAXIMUM
C INTERPOLATION ERROR RELATIVE TO THE MACHINE PRECISION IS
C ALSO PRINTED.  THIS SHOULD BE O(1).  THE INTERPOLATED
C VALUES FROM QS2VAL AND QS2GRD ARE COMPARED FOR AGREEMENT.
C
      INTEGER LCELL(3,3), LNEXT(36)
      REAL    X(36), Y(36), F(36), RSQ(36), A(5,36), P(10)
C
C QSHEP2 PARAMETERS AND LOGICAL UNIT FOR OUTPUT
C
      DATA N/36/,  NQ/13/,  NW/19/,  NR/3/,  LOUT/6/
C
C QUADRATIC TEST FUNCTION AND PARTIAL DERIVATIVES
C
      FQ(XX,YY) = ((XX + 2.*YY)/3.)**2
      FX(XX,YY) = 2.*(XX + 2.*YY)/9.
      FY(XX,YY) = 4.*(XX + 2.*YY)/9.
C
C GENERATE A 6 BY 6 GRID OF NODES IN THE UNIT SQUARE WITH
C   THE NATURAL ORDERING.
C
      K = 0
      DO 1 J = 1,6
	YK = FLOAT(6-J)/5.
	DO 1 I = 1,6
	  K = K + 1
	  X(K) = FLOAT(I-1)/5.
    1	  Y(K) = YK
C
C COMPUTE THE DATA VALUES.
C
      DO 2 K = 1,N
    2	F(K) = FQ(X(K),Y(K))
C
C COMPUTE PARAMETERS DEFINING THE INTERPOLANT Q.
C
      CALL QSHEP2 (N,X,Y,F,NQ,NW,NR, LCELL,LNEXT,XMIN,YMIN,
     .		   DX,DY,RMAX,RSQ,A,IER)
      IF (IER .NE. 0) GO TO 6
C
C GENERATE A 10 BY 10 UNIFORM GRID OF INTERPOLATION POINTS
C   (P(I),P(J)) IN THE UNIT SQUARE.  THE FOUR CORNERS COIN-
C   CIDE WITH NODES.
C
      DO 3 I = 1,10
    3	P(I) = FLOAT(I-1)/9.
C
C COMPUTE THE MACHINE PRECISION EPS.
C
      EPS = 1.
    4 EPS = EPS/2.
      EP1 = EPS + 1.
      IF (STORE(EP1) .GT. 1.) GO TO 4
      EPS = EPS*2.
C
C COMPUTE INTERPOLATION ERRORS AND TEST FOR AGREEMENT IN THE
C   Q VALUES RETURNED BY QS2VAL AND QS2GRD.
C
      EQ = 0.
      EQX = 0.
      EQY = 0.
      DO 5 J = 1,10
	PY = P(J)
	DO 5 I = 1,10
	  PX = P(I)
	  Q1 = QS2VAL (PX,PY,N,X,Y,F,NR,LCELL,LNEXT,XMIN,
     .		       YMIN,DX,DY,RMAX,RSQ,A)
	  CALL QS2GRD (PX,PY,N,X,Y,F,NR,LCELL,LNEXT,XMIN,
     .		       YMIN,DX,DY,RMAX,RSQ,A, Q,QX,QY,IER)
	  IF (IER .NE. 0) GO TO 7
	  IF (ABS(Q1-Q) .GT. 3.*ABS(Q)*EPS) GO TO 8
	  EQ = AMAX1(EQ,ABS(FQ(PX,PY)-Q))
	  EQX = AMAX1(EQX,ABS(FX(PX,PY)-QX))
    5	  EQY = AMAX1(EQY,ABS(FY(PX,PY)-QY))
C
C PRINT ERRORS AND THE RATIO EQ/EPS.
C
      RQ = EQ/EPS
      WRITE (LOUT,100)
      WRITE (LOUT,110) EQ, RQ
      WRITE (LOUT,120) EQX
      WRITE (LOUT,130) EQY
      STOP
  100 FORMAT (///1H ,31HMAXIMUM ABSOLUTE ERRORS IN THE ,
     .	      25HINTERPOLANT Q AND PARTIAL/
     .	      1H ,31HDERIVATIVES QX AND QY RELATIVE ,
     .	      24HTO MACHINE PRECISION EPS//
     .	      1H ,10X,8HFUNCTION,3X,9HMAX ERROR,3X,
     .	      13HMAX ERROR/EPS/)
  110 FORMAT (1H ,13X,1HQ,7X,E9.3,7X,F4.2)
  120 FORMAT (1H ,13X,2HQX,6X,E9.3)
  130 FORMAT (1H ,13X,2HQY,6X,E9.3)
C
C ERROR IN QSHEP2
C
    6 WRITE (LOUT,200) IER
      STOP
  200 FORMAT (///1H ,28H*** ERROR IN QSHEP2 -- IER =,I2,
     .	      4H ***)
C
C ERROR IN QS2GRD
C
    7 WRITE (LOUT,210) IER
      STOP
  210 FORMAT (///1H ,28H*** ERROR IN QS2GRD -- IER =,I2,
     .	      4H ***)
C
C VALUES RETURNED BY QS2VAL AND QS2GRD DIFFER BY A RELATIVE
C   AMOUNT GREATER THAN 3*EPS.
C
    8 WRITE (LOUT,220) Q1, Q
      STOP
  220 FORMAT (///1H ,33H*** ERROR -- INTERPOLATED VALUES ,
     .	      37HQ1 (QS2VAL) AND Q2 (QS2GRD) DIFFER --//
     .	      1H ,5X,5HQ1 = ,E21.14,5X,5HQ2 = ,E21.14)
      END
      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 6 J = 1,IRM1
	    JP1 = J + 1
	    CALL GIVENS (B(J,J),B(J,IROW),C,S)
    6	    CALL ROTATE (6-J,C,S,B(JP1,J),B(JP1,IROW))
	  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 12 I = 1,3
	  B(I,6) = SF
	  IP1 = I + 1
	  DO 11 J = IP1,6
   11	    B(J,6) = 0.
	  DO 12 J = I,5
	    JP1 = J + 1
	    CALL GIVENS (B(J,J),B(J,6),C,S)
   12	    CALL ROTATE (6-J,C,S,B(JP1,J),B(JP1,6))
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 15 IB = 1,5
	  I = 6-IB
	  T = 0.
	  IF (I .EQ. 5) GO TO 15
	  IP1 = I + 1
	  DO 14 J = IP1,5
   14	    T = T + B(J,I)*A(J,K)
   15	  A(I,K) = (B(6,I)-T)/B(I,I)
C
C SCALE THE COEFFICIENTS TO ADJUST FOR THE COLUMN SCALING.
C
	DO 16 I = 1,3
   16	  A(I,K) = A(I,K)/AVSQ
	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 17 I = 1,LNP
	  NP = NPTS(I)
   17	  LNEXT(NP) = -LNEXT(NP)
   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
      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 3 J = JMIN,JMAX
	DO 3 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
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 3 J = JMIN,JMAX
	DO 3 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
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 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 1 I = 1,N
	XI = X(I)
	YI = Y(I)
	X(I) = C*XI + S*YI
	Y(I) = -S*XI + C*YI
    1	CONTINUE
      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 2 I = 1,6
    2	ROW(I) = 0.
      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 1 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)
    1	IF (Y(K) .GT. YMX) YMX = Y(K)
      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 2 J = 1,NNR
	DO 2 I = 1,NNR
    2	  LCELL(I,J) = 0
C
C LOOP ON NODES, STORING INDICES IN LCELL AND LNEXT.
C
      NP1 = NN + 1
      DO 3 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
    3	LCELL(I,J) = KB
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