ac99610df7155bcfb2f1b337781ac726f38b5977
[scilab.git] / scilab / modules / interpolation / src / fortran / cshep2d.f
1       SUBROUTINE CSHEP2 (N,X,Y,F,NC,NW,NR, LCELL,LNEXT,XMIN,
2      .                   YMIN,DX,DY,RMAX,RW,A,IER)
3       INTEGER N, NC, NW, NR, LCELL(NR,NR), LNEXT(N), IER
4       DOUBLE PRECISION  X(N), Y(N), F(N), XMIN, YMIN, DX,
5      .                  DY, RMAX, RW(N), A(9,N)
6 C
7 C***********************************************************
8 C
9 C                                               From CSHEP2D
10 C                                            Robert J. Renka
11 C                                  Dept. of Computer Science
12 C                                       Univ. of North Texas
13 C                                           renka@cs.unt.edu
14 C                                                   02/13/97
15 C
16 C   This subroutine computes a set of parameters defining a
17 C C2 (twice continuously differentiable) bivariate function
18 C C(X,Y) which interpolates data values F at a set of N
19 C arbitrarily distributed points (X,Y) in the plane (nodes).
20 C The interpolant C may be evaluated at an arbitrary point
21 C by function CS2VAL, and its first partial derivatives are
22 C computed by Subroutine CS2GRD.
23 C
24 C   The interpolation scheme is a modified Cubic Shepard
25 C method:
26 C
27 C C = [W(1)*C(1)+W(2)*C(2)+..+W(N)*C(N)]/[W(1)+W(2)+..+W(N)]
28 C
29 C for bivariate functions W(k) and C(k).  The nodal func-
30 C tions are given by
31 C
32 C  C(k)(x,y) = A(1,k)*(x-X(k))**3 +
33 C              A(2,k)*(x-X(k))**2*(y-Y(k)) +
34 C              A(3,k)*(x-X(k))*(y-Y(k))**2 +
35 C              A(4,k)*(y-Y(k))**3 + A(5,k)*(x-X(k))**2 +
36 C              A(6,k)*(x-X(k))*(y-Y(k)) + A(7,k)*(y-Y(k))**2
37 C              + A(8,k)*(x-X(k)) + A(9,k)*(y-Y(k)) + F(k) .
38 C
39 C Thus, C(k) is a cubic function which interpolates the data
40 C value at node k.  Its coefficients A(,k) are obtained by a
41 C weighted least squares fit to the closest NC data points
42 C with weights similar to W(k).  Note that the radius of
43 C influence for the least squares fit is fixed for each k,
44 C but varies with k.
45 C
46 C The weights are taken to be
47 C
48 C   W(k)(x,y) = ( (R(k)-D(k))+ / R(k)*D(k) )**3 ,
49 C
50 C where (R(k)-D(k))+ = 0 if R(k) < D(k), and D(k)(x,y) is
51 C the Euclidean distance between (x,y) and (X(k),Y(k)).  The
52 C radius of influence R(k) varies with k and is chosen so
53 C that NW nodes are within the radius.  Note that W(k) is
54 C not defined at node (X(k),Y(k)), but C(x,y) has limit F(k)
55 C as (x,y) approaches (X(k),Y(k)).
56 C
57 C On input:
58 C
59 C       N = Number of nodes and data values.  N .GE. 10.
60 C
61 C       X,Y = Arrays of length N containing the Cartesian
62 C             coordinates of the nodes.
63 C
64 C       F = Array of length N containing the data values
65 C           in one-to-one correspondence with the nodes.
66 C
67 C       NC = Number of data points to be used in the least
68 C            squares fit for coefficients defining the nodal
69 C            functions C(k).  Values found to be optimal for
70 C            test data sets ranged from 11 to 25.  A recom-
71 C            mended value for general data sets is NC = 17.
72 C            For nodes lying on (or close to) a rectangular
73 C            grid, the recommended value is NC = 11.  In any
74 C            case, NC must be in the range 9 to Min(40,N-1).
75 C
76 C       NW = Number of nodes within (and defining) the radii
77 C            of influence R(k) which enter into the weights
78 C            W(k).  For N sufficiently large, a recommended
79 C            value is NW = 30.  In general, NW should be
80 C            about 1.5*NC.  1 .LE. NW .LE. Min(40,N-1).
81 C
82 C       NR = Number of rows and columns in the cell grid de-
83 C            fined in Subroutine STORE2.  A rectangle con-
84 C            taining the nodes is partitioned into cells in
85 C            order to increase search efficiency.  NR =
86 C            Sqrt(N/3) is recommended.  NR .GE. 1.
87 C
88 C The above parameters are not altered by this routine.
89 C
90 C       LCELL = Array of length .GE. NR**2.
91 C
92 C       LNEXT = Array of length .GE. N.
93 C
94 C       RW = Array of length .GE. N.
95 C
96 C       A = Array of length .GE. 9N.
97 C
98 C On output:
99 C
100 C       LCELL = NR by NR array of nodal indexes associated
101 C               with cells.  Refer to Subroutine STORE2.
102 C
103 C       LNEXT = Array of length N containing next-node
104 C               indexes.  Refer to Subroutine STORE2.
105 C
106 C       XMIN,YMIN,DX,DY = Minimum nodal coordinates and cell
107 C                         dimensions.  Refer to Subroutine
108 C                         STORE2.
109 C
110 C       RMAX = Largest element in RW -- maximum radius R(k).
111 C
112 C       RW = Array containing the radii R(k) which enter
113 C            into the weights W(k).
114 C
115 C       A = 9 by N array containing the coefficients for
116 C           cubic nodal function C(k) in column k.
117 C
118 C   Note that the output parameters described above are not
119 C defined unless IER = 0.
120 C
121 C       IER = Error indicator:
122 C             IER = 0 if no errors were encountered.
123 C             IER = 1 if N, NC, NW, or NR is outside its
124 C                     valid range.
125 C             IER = 2 if duplicate nodes were encountered.
126 C             IER = 3 if all nodes are collinear.
127 C
128 C Modules required by CSHEP2:  GETNP2, GIVENS, ROTATE,
129 C                                SETUP2, STORE2
130 C
131 C Intrinsic functions called by CSHEP2:  ABS, DBLE, MAX,
132 C                                          MIN, SQRT
133 C
134 C***********************************************************
135 C
136       INTEGER LMX
137       PARAMETER (LMX=40)
138       INTEGER I, IERR, IP1, IRM1, IROW, J, JP1, K, LMAX,
139      .        LNP, NEQ, NN, NNC, NNR, NNW, NP, NPTS(LMX),
140      .        NCWMAX
141       DOUBLE PRECISION B(10,10), C, DDX, DDY, DMIN, DTOL,
142      .                 FK, RC, RS, RSMX, RSOLD, RTOL, RWS,
143      .                 S, SF, SFC, SFS, STF, SUM, T, XK,
144      .                 XMN, YK, YMN
145 C
146       DATA    RTOL/1.D-5/, DTOL/.01/
147 C
148 C Local parameters:
149 C
150 C B =          Transpose of the augmented regression matrix
151 C C =          First component of the plane rotation used to
152 C                zero the lower triangle of B**T -- computed
153 C                by Subroutine GIVENS
154 C DDX,DDY =    Local variables for DX and DY
155 C DMIN =       Minimum of the magnitudes of the diagonal
156 C                elements of the regression matrix after
157 C                zeros are introduced below the diagonal
158 C DTOL =       Tolerance for detecting an ill-conditioned
159 C                system.  The system is accepted when
160 C                DMIN*RC .GE. DTOL.
161 C FK =         Data value at mode K -- F(K)
162 C I =          Index for A, B, and NPTS
163 C IERR =       Error flag for the call to Subroutine STORE2
164 C IP1 =        I+1
165 C IRM1 =       IROW-1
166 C IROW =       Row index for B
167 C J =          Index for A and B
168 C JP1 =        J+1
169 C K =          Nodal function index and column index for A
170 C LMAX =       Maximum number of NPTS elements
171 C LMX =        Maximum value of LMAX
172 C LNP =        Current length of NPTS
173 C NEQ =        Number of equations in the least squares fit
174 C NN,NNC,NNR = Local copies of N, NC, and NR
175 C NNW =        Local copy of NW
176 C NP =         NPTS element
177 C NPTS =       Array containing the indexes of a sequence of
178 C                nodes to be used in the least squares fit
179 C                or to compute RW.  The nodes are ordered
180 C                by distance from K, and the last element
181 C                (usually indexed by LNP) is used only to
182 C                determine RC, or RW(K) if NW > NC.
183 C NCWMAX =     Max(NC,NW)
184 C RC =         Radius of influence which enters into the
185 C                weights for C(K) (see Subroutine SETUP2)
186 C RS =         Squared distance between K and NPTS(LNP) --
187 C                used to compute RC and RW(K)
188 C RSMX =       Maximum squared RW element encountered
189 C RSOLD =      Squared distance between K and NPTS(LNP-1) --
190 C                used to compute a relative change in RS
191 C                between succeeding NPTS elements
192 C RTOL =       Tolerance for detecting a sufficiently large
193 C                relative change in RS.  If the change is
194 C                not greater than RTOL, the nodes are
195 C                treated as being the same distance from K
196 C RWS =        Current squared value of RW(K)
197 C S =          Second component of the plane rotation deter-
198 C                mined by subroutine GIVENS
199 C SF =        Scale factor for the linear terms (columns 8
200 C               and 9) in the least squares fit -- inverse
201 C               of the root-mean-square distance between K
202 C               and the nodes (other than K) in the least
203 C               squares fit
204 C SFS =       Scale factor for the quadratic terms (columns
205 C               5, 6, and 7) in the least squares fit --
206 C               SF*SF
207 C SFC =       Scale factor for the cubic terms (first 4
208 C               columns) in the least squares fit -- SF**3
209 C STF =        Marquardt stabilization factor used to damp
210 C                out the first 4 solution components (third
211 C                partials of the cubic) when the system is
212 C                ill-conditioned.  As STF increases, the
213 C                fitting function approaches a quadratic
214 C                polynomial.
215 C SUM =        Sum of squared Euclidean distances between
216 C                node K and the nodes used in the least
217 C                squares fit (unless additional nodes are
218 C                added for stability)
219 C T =          Temporary variable for accumulating a scalar
220 C                product in the back solve
221 C XK,YK =      Coordinates of node K -- X(K), Y(K)
222 C XMN,YMN =    Local variables for XMIN and YMIN
223 C
224       NN = N
225       NNC = NC
226       NNW = NW
227       NNR = NR
228       NCWMAX = MAX(NNC,NNW)
229       LMAX = MIN(LMX,NN-1)
230       IF (NNC .LT. 9  .OR.  NNW .LT. 1  .OR.  NCWMAX .GT.
231      .    LMAX  .OR.  NNR .LT. 1) GO TO 21
232 C
233 C Create the cell data structure, and initialize RSMX.
234 C
235       CALL STORE2 (NN,X,Y,NNR, LCELL,LNEXT,XMN,YMN,DDX,DDY,
236      .             IERR)
237       IF (IERR .NE. 0) GO TO 23
238       RSMX = 0.
239 C
240 C Outer loop on node K:
241 C
242       DO 16 K = 1,NN
243         XK = X(K)
244         YK = Y(K)
245         FK = F(K)
246 C
247 C Mark node K to exclude it from the search for nearest
248 C   neighbors.
249 C
250         LNEXT(K) = -LNEXT(K)
251 C
252 C Initialize for loop on NPTS.
253 C
254         RS = 0.
255         SUM = 0.
256         RWS = 0.
257         RC = 0.
258         LNP = 0
259 C
260 C Compute NPTS, LNP, RWS, NEQ, RC, and SFS.
261 C
262     1   SUM = SUM + RS
263           IF (LNP .EQ. LMAX) GO TO 2
264           LNP = LNP + 1
265           RSOLD = RS
266           CALL GETNP2 (XK,YK,X,Y,NNR,LCELL,LNEXT,XMN,YMN,
267      .                 DDX,DDY, NP,RS)
268           IF (RS .EQ. 0.) GO TO 22
269           NPTS(LNP) = NP
270           IF ( (RS-RSOLD)/RS .LT. RTOL ) GO TO 1
271           IF (RWS .EQ. 0.  .AND.  LNP .GT. NNW) RWS = RS
272           IF (RC .EQ. 0.  .AND.  LNP .GT. NNC) THEN
273 C
274 C   RC = 0 (not yet computed) and LNP > NC.  RC = Sqrt(RS)
275 C     is sufficiently large to (strictly) include NC nodes.
276 C     The least squares fit will include NEQ = LNP - 1
277 C     equations for 9 .LE. NC .LE. NEQ .LT. LMAX .LE. N-1.
278 C
279             NEQ = LNP - 1
280             RC = SQRT(RS)
281             SFS = DBLE(NEQ)/SUM
282           ENDIF
283 C
284 C   Bottom of loop -- test for termination.
285 C
286           IF (LNP .GT. NCWMAX) GO TO 3
287           GO TO 1
288 C
289 C All LMAX nodes are included in NPTS.  RWS and/or RC**2 is
290 C   (arbitrarily) taken to be 10 percent larger than the
291 C   distance RS to the last node included.
292 C
293     2   IF (RWS .EQ. 0.) RWS = 1.1*RS
294         IF (RC .EQ. 0.) THEN
295           NEQ = LMAX
296           RC = SQRT(1.1*RS)
297           SFS = DBLE(NEQ)/SUM
298         ENDIF
299 C
300 C Store RW(K), update RSMX if necessary, and compute SF
301 C   and SFC.
302 C
303     3   RW(K) = SQRT(RWS)
304         IF (RWS .GT. RSMX) RSMX = RWS
305         SF = SQRT(SFS)
306         SFC = SF*SFS
307 C
308 C A Q-R decomposition is used to solve the least squares
309 C   system.  The transpose of the augmented regression
310 C   matrix is stored in B with columns (rows of B) defined
311 C   as follows:  1-4 are the cubic terms, 5-7 are the quad-
312 C   ratic terms, 8 and 9 are the linear terms, and the last
313 C   column is the right hand side.
314 C
315 C Set up the equations and zero out the lower triangle with
316 C   Givens rotations.
317 C
318         I = 0
319     4     I = I + 1
320           NP = NPTS(I)
321           IROW = MIN(I,10)
322           CALL SETUP2 (XK,YK,FK,X(NP),Y(NP),F(NP),SF,SFS,
323      .                 SFC,RC, B(1,IROW))
324           IF (I .EQ. 1) GO TO 4
325           IRM1 = IROW-1
326           DO 5 J = 1,IRM1
327             JP1 = J + 1
328             CALL GIVENS (B(J,J),B(J,IROW),C,S)
329             CALL ROTATE (10-J,C,S,B(JP1,J),B(JP1,IROW))
330     5       CONTINUE
331           IF (I .LT. NEQ) GO TO 4
332 C
333 C Test the system for ill-conditioning.
334 C
335         DMIN = MIN( ABS(B(1,1)),ABS(B(2,2)),ABS(B(3,3)),
336      .              ABS(B(4,4)),ABS(B(5,5)),ABS(B(6,6)),
337      .              ABS(B(7,7)),ABS(B(8,8)),ABS(B(9,9)) )
338         IF (DMIN*RC .GE. DTOL) GO TO 11
339         IF (NEQ .EQ. LMAX) GO TO 7
340 C
341 C Increase RC and add another equation to the system to
342 C   improve the conditioning.  The number of NPTS elements
343 C   is also increased if necessary.
344 C
345     6   RSOLD = RS
346         NEQ = NEQ + 1
347         IF (NEQ .EQ. LMAX) THEN
348           RC = SQRT(1.1*RS)
349           GO TO 4
350         ENDIF
351         IF (NEQ .LT. LNP) THEN
352 C
353 C   NEQ < LNP.
354 C
355           NP = NPTS(NEQ+1)
356           RS = (X(NP)-XK)**2 + (Y(NP)-YK)**2
357           IF ( (RS-RSOLD)/RS .LT. RTOL ) GO TO 6
358           RC = SQRT(RS)
359           GO TO 4
360         ENDIF
361 C
362 C   NEQ = LNP.  Add an element to NPTS.
363 C
364         LNP = LNP + 1
365         CALL GETNP2 (XK,YK,X,Y,NNR,LCELL,LNEXT,XMN,YMN,
366      .               DDX,DDY, NP,RS)
367         IF (NP .EQ. 0) GO TO 22
368         NPTS(LNP) = NP
369         IF ( (RS-RSOLD)/RS .LT. RTOL ) GO TO 6
370         RC = SQRT(RS)
371         GO TO 4
372 C
373 C Stabilize the system by damping third partials -- add
374 C   multiples of the first four unit vectors to the first
375 C   four equations.
376 C
377     7   STF = 1.0/RC
378         DO 10 I = 1,4
379           B(I,10) = STF
380           IP1 = I + 1
381           DO 8 J = IP1,10
382             B(J,10) = 0.
383     8       CONTINUE
384           DO 9 J = I,9
385             JP1 = J + 1
386             CALL GIVENS (B(J,J),B(J,10),C,S)
387             CALL ROTATE (10-J,C,S,B(JP1,J),B(JP1,10))
388     9       CONTINUE
389    10     CONTINUE
390 C
391 C Test the damped system for ill-conditioning.
392 C
393         DMIN = MIN( ABS(B(5,5)),ABS(B(6,6)),ABS(B(7,7)),
394      .              ABS(B(8,8)),ABS(B(9,9)) )
395         IF (DMIN*RC .LT. DTOL) GO TO 23
396 C
397 C Solve the 9 by 9 triangular system for the coefficients.
398 C
399    11   DO 13 I = 9,1,-1
400           T = 0.
401           IF (I .NE. 9) THEN
402             IP1 = I + 1
403             DO 12 J = IP1,9
404               T = T + B(J,I)*A(J,K)
405    12         CONTINUE
406           ENDIF
407           A(I,K) = (B(10,I)-T)/B(I,I)
408    13     CONTINUE
409 C
410 C Scale the coefficients to adjust for the column scaling.
411 C
412         DO 14 I = 1,4
413           A(I,K) = A(I,K)*SFC
414    14     CONTINUE
415         A(5,K) = A(5,K)*SFS
416         A(6,K) = A(6,K)*SFS
417         A(7,K) = A(7,K)*SFS
418         A(8,K) = A(8,K)*SF
419         A(9,K) = A(9,K)*SF
420 C
421 C Unmark K and the elements of NPTS.
422 C
423         LNEXT(K) = -LNEXT(K)
424         DO 15 I = 1,LNP
425           NP = NPTS(I)
426           LNEXT(NP) = -LNEXT(NP)
427    15     CONTINUE
428    16   CONTINUE
429 C
430 C No errors encountered.
431 C
432       XMIN = XMN
433       YMIN = YMN
434       DX = DDX
435       DY = DDY
436       RMAX = SQRT(RSMX)
437       IER = 0
438       RETURN
439 C
440 C N, NC, NW, or NR is outside its valid range.
441 C
442    21 IER = 1
443       RETURN
444 C
445 C Duplicate nodes were encountered by GETNP2.
446 C
447    22 IER = 2
448       RETURN
449 C
450 C No unique solution due to collinear nodes.
451 C
452    23 XMIN = XMN
453       YMIN = YMN
454       DX = DDX
455       DY = DDY
456       IER = 3
457       RETURN
458       END
459
460       DOUBLE PRECISION FUNCTION CS2VAL (PX,PY,N,X,Y,F,NR,
461      .                LCELL,LNEXT,XMIN,YMIN,DX,DY,RMAX,RW,A)
462       INTEGER N, NR, LCELL(NR,NR), LNEXT(N)
463       DOUBLE PRECISION PX, PY, X(N), Y(N), F(N), XMIN, YMIN,
464      .                 DX, DY, RMAX, RW(N), A(9,N)
465 C
466 C***********************************************************
467 C
468 C                                               From CSHEP2D
469 C                                            Robert J. Renka
470 C                                  Dept. of Computer Science
471 C                                       Univ. of North Texas
472 C                                           renka@cs.unt.edu
473 C                                                   02/03/97
474 C
475 C   This function returns the value C(PX,PY), where C is the
476 C weighted sum of cubic nodal functions defined in Subrou-
477 C tine CSHEP2.  CS2GRD may be called to compute a gradient
478 C of C along with the value, and/or to test for errors.
479 C CS2HES may be called to compute a value, first partial
480 C derivatives, and second partial derivatives at a point.
481 C
482 C On input:
483 C
484 C       PX,PY = Cartesian coordinates of the point P at
485 C               which C is to be evaluated.
486 C
487 C       N = Number of nodes and data values defining C.
488 C           N .GE. 10.
489 C
490 C       X,Y,F = Arrays of length N containing the nodes and
491 C               data values interpolated by C.
492 C
493 C       NR = Number of rows and columns in the cell grid.
494 C            Refer to Subroutine STORE2.  NR .GE. 1.
495 C
496 C       LCELL = NR by NR array of nodal indexes associated
497 C               with cells.  Refer to Subroutine STORE2.
498 C
499 C       LNEXT = Array of length N containing next-node
500 C               indexes.  Refer to Subroutine STORE2.
501 C
502 C       XMIN,YMIN,DX,DY = Minimum nodal coordinates and cell
503 C                         dimensions.  DX and DY must be
504 C                         positive.  Refer to Subroutine
505 C                         STORE2.
506 C
507 C       RMAX = Largest element in RW -- maximum radius R(k).
508 C
509 C       RW = Array containing the radii R(k) which enter
510 C            into the weights W(k) defining C.
511 C
512 C       A = 9 by N array containing the coefficients for
513 C           cubic nodal function C(k) in column k.
514 C
515 C   Input parameters are not altered by this function.  The
516 C parameters other than PX and PY should be input unaltered
517 C from their values on output from CSHEP2.  This function
518 C should not be called if a nonzero error flag was returned
519 C by CSHEP2.
520 C
521 C On output:
522 C
523 C       CS2VAL = Function value C(PX,PY) unless N, NR, DX,
524 C                DY, or RMAX is invalid, in which case no
525 C                value is returned.
526 C
527 C Modules required by CS2VAL:  NONE
528 C
529 C Intrinsic functions called by CS2VAL:  INT, SQRT
530 C
531 C***********************************************************
532 C
533       INTEGER I, IMAX, IMIN, J, JMAX, JMIN, K, KP
534       DOUBLE PRECISION D, DELX, DELY, R, SW, SWC, W, XP, YP
535 C
536 C Local parameters:
537 C
538 C D =         Distance between P and node K
539 C DELX =      XP - X(K)
540 C DELY =      YP - Y(K)
541 C I =         Cell row index in the range IMIN to IMAX
542 C IMIN,IMAX = Range of cell row indexes of the cells
543 C               intersected by a disk of radius RMAX
544 C               centered at P
545 C J =         Cell column index in the range JMIN to JMAX
546 C JMIN,JMAX = Range of cell column indexes of the cells
547 C               intersected by a disk of radius RMAX
548 C               centered at P
549 C K =         Index of a node in cell (I,J)
550 C KP =        Previous value of K in the sequence of nodes
551 C               in cell (I,J)
552 C R =         Radius of influence for node K
553 C SW =        Sum of weights W(K)
554 C SWC =       Sum of weighted nodal function values at P
555 C W =         Weight W(K) value at P:  ((R-D)+/(R*D))**3,
556 C               where (R-D)+ = 0 if R < D
557 C XP,YP =     Local copies of PX and PY -- coordinates of P
558 C
559       XP = PX
560       YP = PY
561       IF (N .LT. 10  .OR.  NR .LT. 1  .OR.  DX .LE. 0.  .OR.
562      .    DY .LE. 0.  .OR.  RMAX .LT. 0.) RETURN
563 C
564 C Set IMIN, IMAX, JMIN, and JMAX to cell indexes defining
565 C   the range of the search for nodes whose radii include
566 C   P.  The cells which must be searched are those inter-
567 C   sected by (or contained in) a circle of radius RMAX
568 C   centered at P.
569 C
570       IMIN = INT((XP-XMIN-RMAX)/DX) + 1
571       IMAX = INT((XP-XMIN+RMAX)/DX) + 1
572       IF (IMIN .LT. 1) IMIN = 1
573       IF (IMAX .GT. NR) IMAX = NR
574       JMIN = INT((YP-YMIN-RMAX)/DY) + 1
575       JMAX = INT((YP-YMIN+RMAX)/DY) + 1
576       IF (JMIN .LT. 1) JMIN = 1
577       IF (JMAX .GT. NR) JMAX = NR
578 C
579 C The following is a test for no cells within the circle
580 C   of radius RMAX.
581 C
582       IF (IMIN .GT. IMAX  .OR.  JMIN .GT. JMAX) GO TO 6
583 C
584 C Accumulate weight values in SW and weighted nodal function
585 C   values in SWC.  The weights are W(K) = ((R-D)+/(R*D))**3
586 C   for R = RW(K) and D = distance between P and node K.
587 C
588       SW = 0.
589       SWC = 0.
590 C
591 C Outer loop on cells (I,J).
592 C
593       DO 4 J = JMIN,JMAX
594         DO 3 I = IMIN,IMAX
595           K = LCELL(I,J)
596           IF (K .EQ. 0) GO TO 3
597 C
598 C Inner loop on nodes K.
599 C
600     1     DELX = XP - X(K)
601           DELY = YP - Y(K)
602           D = SQRT(DELX*DELX + DELY*DELY)
603           R = RW(K)
604           IF (D .GE. R) GO TO 2
605           IF (D .EQ. 0.) GO TO 5
606           W = (1.0/D - 1.0/R)**3
607           SW = SW + W
608           SWC = SWC + W*( ( (A(1,K)*DELX+A(2,K)*DELY+
609      .                       A(5,K))*DELX + (A(3,K)*DELY+
610      .                       A(6,K))*DELY + A(8,K) )*DELX +
611      .                    ( (A(4,K)*DELY+A(7,K))*DELY +
612      .                      A(9,K) )*DELY + F(K) )
613 C
614 C Bottom of loop on nodes in cell (I,J).
615 C
616     2     KP = K
617           K = LNEXT(KP)
618           IF (K .NE. KP) GO TO 1
619     3     CONTINUE
620     4   CONTINUE
621 C
622 C SW = 0 iff P is not within the radius R(K) for any node K.
623 C
624       IF (SW .EQ. 0.) GO TO 6
625       CS2VAL = SWC/SW
626       RETURN
627 C
628 C (PX,PY) = (X(K),Y(K)).
629 C
630     5 CS2VAL = F(K)
631       RETURN
632 C
633 C All weights are 0 at P.
634 C
635     6 CS2VAL = 0.
636       RETURN
637       END
638
639       SUBROUTINE CS2GRD (PX,PY,N,X,Y,F,NR,LCELL,LNEXT,XMIN,
640      .                   YMIN,DX,DY,RMAX,RW,A, C,CX,CY,IER)
641       INTEGER N, NR, LCELL(NR,NR), LNEXT(N), IER
642       DOUBLE PRECISION PX, PY, X(N), Y(N), F(N), XMIN, YMIN,
643      .                 DX, DY, RMAX, RW(N), A(9,N), C, CX,
644      .                 CY
645 C
646 C***********************************************************
647 C
648 C                                               From CSHEP2D
649 C                                            Robert J. Renka
650 C                                  Dept. of Computer Science
651 C                                       Univ. of North Texas
652 C                                           renka@cs.unt.edu
653 C                                                   02/03/97
654 C
655 C   This subroutine computes the value and gradient at P =
656 C (PX,PY) of the interpolatory function C defined in Sub-
657 C routine CSHEP2.  C is a weighted sum of cubic nodal
658 C functions.
659 C
660 C On input:
661 C
662 C       PX,PY = Cartesian coordinates of the point P at
663 C               which C and its partial derivatives are
664 C               to be evaluated.
665 C
666 C       N = Number of nodes and data values defining C.
667 C           N .GE. 10.
668 C
669 C       X,Y,F = Arrays of length N containing the nodes and
670 C               data values interpolated by C.
671 C
672 C       NR = Number of rows and columns in the cell grid.
673 C            Refer to Subroutine STORE2.  NR .GE. 1.
674 C
675 C       LCELL = NR by NR array of nodal indexes associated
676 C               with cells.  Refer to Subroutine STORE2.
677 C
678 C       LNEXT = Array of length N containing next-node
679 C               indexes.  Refer to Subroutine STORE2.
680 C
681 C       XMIN,YMIN,DX,DY = Minimum nodal coordinates and cell
682 C                         dimensions.  DX and DY must be
683 C                         positive.  Refer to Subroutine
684 C                         STORE2.
685 C
686 C       RMAX = Largest element in RW -- maximum radius R(k).
687 C
688 C       RW = Array of length N containing the radii R(k)
689 C            which enter into the weights W(k) defining C.
690 C
691 C       A = 9 by N array containing the coefficients for
692 C           cubic nodal function C(k) in column k.
693 C
694 C   Input parameters are not altered by this subroutine.
695 C The parameters other than PX and PY should be input
696 C unaltered from their values on output from CSHEP2.  This
697 C subroutine should not be called if a nonzero error flag
698 C was returned by CSHEP2.
699 C
700 C On output:
701 C
702 C       C = Value of C at (PX,PY) unless IER .EQ. 1, in
703 C           which case no values are returned.
704 C
705 C       CX,CY = First partial derivatives of C at (PX,PY)
706 C               unless IER .EQ. 1.
707 C
708 C       IER = Error indicator:
709 C             IER = 0 if no errors were encountered.
710 C             IER = 1 if N, NR, DX, DY or RMAX is invalid.
711 C             IER = 2 if no errors were encountered but
712 C                     (PX,PY) is not within the radius R(k)
713 C                     for any node k (and thus C=CX=CY=0).
714 C
715 C Modules required by CS2GRD:  None
716 C
717 C Intrinsic functions called by CS2GRD:  INT, SQRT
718 C
719 C***********************************************************
720 C
721       INTEGER I, IMAX, IMIN, J, JMAX, JMIN, K, KP
722       DOUBLE PRECISION CK, CKX, CKY, D, DELX, DELY, R, SW,
723      .                 SWC, SWCX, SWCY, SWS, SWX, SWY, T, W,
724      .                 WX, WY, XP, YP
725 C
726 C Local parameters:
727 C
728 C CK =        Value of cubic nodal function C(K) at P
729 C CKX,CKY =   Partial derivatives of C(K) with respect to X
730 C               and Y, respectively
731 C D =         Distance between P and node K
732 C DELX =      XP - X(K)
733 C DELY =      YP - Y(K)
734 C I =         Cell row index in the range IMIN to IMAX
735 C IMIN,IMAX = Range of cell row indexes of the cells
736 C               intersected by a disk of radius RMAX
737 C               centered at P
738 C J =         Cell column index in the range JMIN to JMAX
739 C JMIN,JMAX = Range of cell column indexes of the cells
740 C               intersected by a disk of radius RMAX
741 C               centered at P
742 C K =         Index of a node in cell (I,J)
743 C KP =        Previous value of K in the sequence of nodes
744 C               in cell (I,J)
745 C R =         Radius of influence for node K
746 C SW =        Sum of weights W(K)
747 C SWC =       Sum of weighted nodal function values at P
748 C SWCX,SWCY = Partial derivatives of SWC with respect to X
749 C               and Y, respectively
750 C SWS =       SW**2
751 C SWX,SWY =   Partial derivatives of SW with respect to X
752 C               and Y, respectively
753 C T =         Temporary variable
754 C W =         Weight W(K) value at P:  ((R-D)+/(R*D))**3,
755 C               where (R-D)+ = 0 if R < D
756 C WX,WY =     Partial derivatives of W with respect to X
757 C               and Y, respectively
758 C XP,YP =     Local copies of PX and PY -- coordinates of P
759 C
760       XP = PX
761       YP = PY
762       IF (N .LT. 10  .OR.  NR .LT. 1  .OR.  DX .LE. 0.  .OR.
763      .    DY .LE. 0.  .OR.  RMAX .LT. 0.) GO TO 6
764 C
765 C Set IMIN, IMAX, JMIN, and JMAX to cell indexes defining
766 C   the range of the search for nodes whose radii include
767 C   P.  The cells which must be searched are those inter-
768 C   sected by (or contained in) a circle of radius RMAX
769 C   centered at P.
770 C
771       IMIN = INT((XP-XMIN-RMAX)/DX) + 1
772       IMAX = INT((XP-XMIN+RMAX)/DX) + 1
773       IF (IMIN .LT. 1) IMIN = 1
774       IF (IMAX .GT. NR) IMAX = NR
775       JMIN = INT((YP-YMIN-RMAX)/DY) + 1
776       JMAX = INT((YP-YMIN+RMAX)/DY) + 1
777       IF (JMIN .LT. 1) JMIN = 1
778       IF (JMAX .GT. NR) JMAX = NR
779 C
780 C The following is a test for no cells within the circle
781 C   of radius RMAX.
782 C
783       IF (IMIN .GT. IMAX  .OR.  JMIN .GT. JMAX) GO TO 7
784 C
785 C C = SWC/SW = Sum(W(K)*C(K))/Sum(W(K)), where the sum is
786 C   from K = 1 to N, C(K) is the cubic nodal function value,
787 C   and W(K) = ((R-D)+/(R*D))**3 for radius R(K) and dist-
788 C   ance D(K).  Thus
789 C
790 C        CX = (SWCX*SW - SWC*SWX)/SW**2  and
791 C        CY = (SWCY*SW - SWC*SWY)/SW**2
792 C
793 C   where SWCX and SWX are partial derivatives with respect
794 C   to X of SWC and SW, respectively.  SWCY and SWY are de-
795 C   fined similarly.
796 C
797       SW = 0.
798       SWX = 0.
799       SWY = 0.
800       SWC = 0.
801       SWCX = 0.
802       SWCY = 0.
803 C
804 C Outer loop on cells (I,J).
805 C
806       DO 4 J = JMIN,JMAX
807         DO 3 I = IMIN,IMAX
808           K = LCELL(I,J)
809           IF (K .EQ. 0) GO TO 3
810 C
811 C Inner loop on nodes K.
812 C
813     1     DELX = XP - X(K)
814           DELY = YP - Y(K)
815           D = SQRT(DELX*DELX + DELY*DELY)
816           R = RW(K)
817           IF (D .GE. R) GO TO 2
818           IF (D .EQ. 0.) GO TO 5
819           T = (1.0/D - 1.0/R)
820           W = T**3
821           T = -3.0*T*T/(D**3)
822           WX = DELX*T
823           WY = DELY*T
824           T = A(2,K)*DELX + A(3,K)*DELY + A(6,K)
825           CKY = ( 3.0*A(4,K)*DELY + A(3,K)*DELX +
826      .            2.0*A(7,K) )*DELY + T*DELX + A(9,K)
827           T = T*DELY + A(8,K)
828           CKX = ( 3.0*A(1,K)*DELX + A(2,K)*DELY +
829      .            2.0*A(5,K) )*DELX + T
830           CK = ( (A(1,K)*DELX+A(5,K))*DELX + T )*DELX +
831      .         ( (A(4,K)*DELY+A(7,K))*DELY + A(9,K) )*DELY +
832      .         F(K)
833           SW = SW + W
834           SWX = SWX + WX
835           SWY = SWY + WY
836           SWC = SWC + W*CK
837           SWCX = SWCX + WX*CK + W*CKX
838           SWCY = SWCY + WY*CK + W*CKY
839 C
840 C Bottom of loop on nodes in cell (I,J).
841 C
842     2     KP = K
843           K = LNEXT(KP)
844           IF (K .NE. KP) GO TO 1
845     3     CONTINUE
846     4   CONTINUE
847 C
848 C SW = 0 iff P is not within the radius R(K) for any node K.
849 C
850       IF (SW .EQ. 0.) GO TO 7
851       C = SWC/SW
852       SWS = SW*SW
853       CX = (SWCX*SW - SWC*SWX)/SWS
854       CY = (SWCY*SW - SWC*SWY)/SWS
855       IER = 0
856       RETURN
857 C
858 C (PX,PY) = (X(K),Y(K)).
859 C
860     5 C = F(K)
861       CX = A(8,K)
862       CY = A(9,K)
863       IER = 0
864       RETURN
865 C
866 C Invalid input parameter.
867 C
868     6 IER = 1
869       RETURN
870 C
871 C No cells contain a point within RMAX of P, or
872 C   SW = 0 and thus D .GE. RW(K) for all K.
873 C
874     7 C = 0.
875       CX = 0.
876       CY = 0.
877       IER = 2
878       RETURN
879       END
880       SUBROUTINE CS2HES (PX,PY,N,X,Y,F,NR,LCELL,LNEXT,XMIN,
881      .                   YMIN,DX,DY,RMAX,RW,A, C,CX,CY,CXX,
882      .                   CXY,CYY,IER)
883       INTEGER N, NR, LCELL(NR,NR), LNEXT(N), IER
884       DOUBLE PRECISION PX, PY, X(N), Y(N), F(N), XMIN, YMIN,
885      .                 DX, DY, RMAX, RW(N), A(9,N), C, CX,
886      .                 CY, CXX, CXY, CYY
887 C
888 C***********************************************************
889 C
890 C                                               From CSHEP2D
891 C                                            Robert J. Renka
892 C                                  Dept. of Computer Science
893 C                                       Univ. of North Texas
894 C                                           renka@cs.unt.edu
895 C                                                   02/03/97
896 C
897 C   This subroutine computes the value, gradient, and
898 C Hessian at P = (PX,PY) of the interpolatory function C
899 C defined in Subroutine CSHEP2.  C is a weighted sum of
900 C cubic nodal functions.
901 C
902 C On input:
903 C
904 C       PX,PY = Cartesian coordinates of the point P at
905 C               which C and its partial derivatives are
906 C               to be evaluated.
907 C
908 C       N = Number of nodes and data values defining C.
909 C           N .GE. 10.
910 C
911 C       X,Y,F = Arrays of length N containing the nodes and
912 C               data values interpolated by C.
913 C
914 C       NR = Number of rows and columns in the cell grid.
915 C            Refer to Subroutine STORE2.  NR .GE. 1.
916 C
917 C       LCELL = NR by NR array of nodal indexes associated
918 C               with cells.  Refer to Subroutine STORE2.
919 C
920 C       LNEXT = Array of length N containing next-node
921 C               indexes.  Refer to Subroutine STORE2.
922 C
923 C       XMIN,YMIN,DX,DY = Minimum nodal coordinates and cell
924 C                         dimensions.  DX and DY must be
925 C                         positive.  Refer to Subroutine
926 C                         STORE2.
927 C
928 C       RMAX = Largest element in RW -- maximum radius R(k).
929 C
930 C       RW = Array of length N containing the radii R(k)
931 C            which enter into the weights W(k) defining C.
932 C
933 C       A = 9 by N array containing the coefficients for
934 C           cubic nodal function C(k) in column k.
935 C
936 C   Input parameters are not altered by this subroutine.
937 C The parameters other than PX and PY should be input
938 C unaltered from their values on output from CSHEP2.  This
939 C subroutine should not be called if a nonzero error flag
940 C was returned by CSHEP2.
941 C
942 C On output:
943 C
944 C       C = Value of C at (PX,PY) unless IER .EQ. 1, in
945 C           which case no values are returned.
946 C
947 C       CX,CY = First partial derivatives of C at (PX,PY)
948 C               unless IER .EQ. 1.
949 C
950 C       CXX,CXY,CYY = Second partial derivatives of C at
951 C                     (PX,PY) unless IER .EQ. 1.
952 C
953 C       IER = Error indicator:
954 C             IER = 0 if no errors were encountered.
955 C             IER = 1 if N, NR, DX, DY or RMAX is invalid.
956 C             IER = 2 if no errors were encountered but
957 C                     (PX,PY) is not within the radius R(k)
958 C                     for any node k (and thus C = 0).
959 C
960 C Modules required by CS2HES:  None
961 C
962 C Intrinsic functions called by CS2HES:  INT, SQRT
963 C
964 C***********************************************************
965 C
966       INTEGER I, IMAX, IMIN, J, JMAX, JMIN, K, KP
967       DOUBLE PRECISION CK, CKX, CKXX, CKXY, CKY, CKYY, D,
968      .                 DELX, DELY, DXSQ, DYSQ, R, SW, SWC,
969      .                 SWCX, SWCXX, SWCXY, SWCY, SWCYY, SWS,
970      .                 SWX, SWXX, SWXY, SWY, SWYY, T1, T2,
971      .                 T3, T4, W, WX, WXX, WXY, WY, WYY, XP,
972      .                 YP
973 C
974 C Local parameters:
975 C
976 C CK =        Value of cubic nodal function C(K) at P
977 C CKX,CKY =   Partial derivatives of C(K) with respect to X
978 C               and Y, respectively
979 C CKXX,CKXY,CKYY = Second partial derivatives of CK
980 C D =         Distance between P and node K
981 C DELX =      XP - X(K)
982 C DELY =      YP - Y(K)
983 C DXSQ,DYSQ = DELX**2, DELY**2
984 C I =         Cell row index in the range IMIN to IMAX
985 C IMIN,IMAX = Range of cell row indexes of the cells
986 C               intersected by a disk of radius RMAX
987 C               centered at P
988 C J =         Cell column index in the range JMIN to JMAX
989 C JMIN,JMAX = Range of cell column indexes of the cells
990 C               intersected by a disk of radius RMAX
991 C               centered at P
992 C K =         Index of a node in cell (I,J)
993 C KP =        Previous value of K in the sequence of nodes
994 C               in cell (I,J)
995 C R =         Radius of influence for node K
996 C SW =        Sum of weights W(K)
997 C SWC =       Sum of weighted nodal function values at P
998 C SWCX,SWCY = Partial derivatives of SWC with respect to X
999 C               and Y, respectively
1000 C SWCXX,SWCXY,SWCYY = Second partial derivatives of SWC
1001 C SWS =       SW**2
1002 C SWX,SWY =   Partial derivatives of SW with respect to X
1003 C               and Y, respectively
1004 C SWXX,SWXY,SWYY = Second partial derivatives of SW
1005 C T1,T2,T3,T4 = Temporary variables
1006 C W =         Weight W(K) value at P:  ((R-D)+/(R*D))**3,
1007 C               where (R-D)+ = 0 if R < D
1008 C WX,WY =     Partial derivatives of W with respect to X
1009 C               and Y, respectively
1010 C WXX,WXY,WYY = Second partial derivatives of W
1011 C XP,YP =     Local copies of PX and PY -- coordinates of P
1012 C
1013       XP = PX
1014       YP = PY
1015       IF (N .LT. 10  .OR.  NR .LT. 1  .OR.  DX .LE. 0.  .OR.
1016      .    DY .LE. 0.  .OR.  RMAX .LT. 0.) GO TO 6
1017 C
1018 C Set IMIN, IMAX, JMIN, and JMAX to cell indexes defining
1019 C   the range of the search for nodes whose radii include
1020 C   P.  The cells which must be searched are those inter-
1021 C   sected by (or contained in) a circle of radius RMAX
1022 C   centered at P.
1023 C
1024       IMIN = INT((XP-XMIN-RMAX)/DX) + 1
1025       IMAX = INT((XP-XMIN+RMAX)/DX) + 1
1026       IF (IMIN .LT. 1) IMIN = 1
1027       IF (IMAX .GT. NR) IMAX = NR
1028       JMIN = INT((YP-YMIN-RMAX)/DY) + 1
1029       JMAX = INT((YP-YMIN+RMAX)/DY) + 1
1030       IF (JMIN .LT. 1) JMIN = 1
1031       IF (JMAX .GT. NR) JMAX = NR
1032 C
1033 C The following is a test for no cells within the circle
1034 C   of radius RMAX.
1035 C
1036       IF (IMIN .GT. IMAX  .OR.  JMIN .GT. JMAX) GO TO 7
1037 C
1038 C C = SWC/SW = Sum(W(K)*C(K))/Sum(W(K)), where the sum is
1039 C   from K = 1 to N, C(K) is the cubic nodal function value,
1040 C   and W(K) = ((R-D)+/(R*D))**3 for radius R(K) and dist-
1041 C   ance D(K).  Thus
1042 C
1043 C        CX = (SWCX*SW - SWC*SWX)/SW**2  and
1044 C        CY = (SWCY*SW - SWC*SWY)/SW**2
1045 C
1046 C   where SWCX and SWX are partial derivatives with respect
1047 C   to x of SWC and SW, respectively.  SWCY and SWY are de-
1048 C   fined similarly.  The second partials are
1049 C
1050 C        CXX = ( SW*(SWCXX -    2*SWX*CX) - SWC*SWXX )/SW**2
1051 C        CXY = ( SW*(SWCXY-SWX*CY-SWY*CX) - SWC*SWXY )/SW**2
1052 C        CYY = ( SW*(SWCYY -    2*SWY*CY) - SWC*SWYY )/SW**2
1053 C
1054 C   where SWCXX and SWXX are second partials with respect
1055 C   to x, SWCXY and SWXY are mixed partials, and SWCYY and
1056 C   SWYY are second partials with respect to y.
1057 C
1058       SW = 0.
1059       SWX = 0.
1060       SWY = 0.
1061       SWXX = 0.
1062       SWXY = 0.
1063       SWYY = 0.
1064       SWC = 0.
1065       SWCX = 0.
1066       SWCY = 0.
1067       SWCXX = 0.
1068       SWCXY = 0.
1069       SWCYY = 0.
1070 C
1071 C Outer loop on cells (I,J).
1072 C
1073       DO 4 J = JMIN,JMAX
1074         DO 3 I = IMIN,IMAX
1075           K = LCELL(I,J)
1076           IF (K .EQ. 0) GO TO 3
1077 C
1078 C Inner loop on nodes K.
1079 C
1080     1     DELX = XP - X(K)
1081           DELY = YP - Y(K)
1082           DXSQ = DELX*DELX
1083           DYSQ = DELY*DELY
1084           D = SQRT(DXSQ + DYSQ)
1085           R = RW(K)
1086           IF (D .GE. R) GO TO 2
1087           IF (D .EQ. 0.) GO TO 5
1088           T1 = (1.0/D - 1.0/R)
1089           W = T1**3
1090           T2 = -3.0*T1*T1/(D**3)
1091           WX = DELX*T2
1092           WY = DELY*T2
1093           T1 = 3.0*T1*(2.0+3.0*D*T1)/(D**6)
1094           WXX = T1*DXSQ + T2
1095           WXY = T1*DELX*DELY
1096           WYY = T1*DYSQ + T2
1097           T1 = A(1,K)*DELX + A(2,K)*DELY + A(5,K)
1098           T2 = T1 + T1 + A(1,K)*DELX
1099           T3 = A(4,K)*DELY + A(3,K)*DELX + A(7,K)
1100           T4 = T3 + T3 + A(4,K)*DELY
1101           CK = (T1*DELX + A(6,K)*DELY + A(8,K))*DELX +
1102      .         (T3*DELY + A(9,K))*DELY + F(K)
1103           CKX = T2*DELX + (A(3,K)*DELY+A(6,K))*DELY + A(8,K)
1104           CKY = T4*DELY + (A(2,K)*DELX+A(6,K))*DELX + A(9,K)
1105           CKXX = T2 + 3.0*A(1,K)*DELX
1106           CKXY = 2.0*(A(2,K)*DELX + A(3,K)*DELY) + A(6,K)
1107           CKYY = T4 + 3.0*A(4,K)*DELY
1108           SW = SW + W
1109           SWX = SWX + WX
1110           SWY = SWY + WY
1111           SWXX = SWXX + WXX
1112           SWXY = SWXY + WXY
1113           SWYY = SWYY + WYY
1114           SWC = SWC + W*CK
1115           SWCX = SWCX + WX*CK + W*CKX
1116           SWCY = SWCY + WY*CK + W*CKY
1117           SWCXX = SWCXX + W*CKXX + 2.0*WX*CKX + CK*WXX
1118           SWCXY = SWCXY + W*CKXY + WX*CKY + WY*CKX + CK*WXY
1119           SWCYY = SWCYY + W*CKYY + 2.0*WY*CKY + CK*WYY
1120 C
1121 C Bottom of loop on nodes in cell (I,J).
1122 C
1123     2     KP = K
1124           K = LNEXT(KP)
1125           IF (K .NE. KP) GO TO 1
1126     3     CONTINUE
1127     4   CONTINUE
1128 C
1129 C SW = 0 iff P is not within the radius R(K) for any node K.
1130 C
1131       IF (SW .EQ. 0.) GO TO 7
1132       C = SWC/SW
1133       SWS = SW*SW
1134       CX = (SWCX*SW - SWC*SWX)/SWS
1135       CY = (SWCY*SW - SWC*SWY)/SWS
1136       CXX = (SW*(SWCXX-2.0*SWX*CX) - SWC*SWXX)/SWS
1137       CXY = (SW*(SWCXY-SWY*CX-SWX*CY) - SWC*SWXY)/SWS
1138       CYY = (SW*(SWCYY-2.0*SWY*CY) - SWC*SWYY)/SWS
1139       IER = 0
1140       RETURN
1141 C
1142 C (PX,PY) = (X(K),Y(K)).
1143 C
1144     5 C = F(K)
1145       CX = A(8,K)
1146       CY = A(9,K)
1147       CXX = 2.0*A(5,K)
1148       CXY = A(6,K)
1149       CYY = 2.0*A(7,K)
1150       IER = 0
1151       RETURN
1152 C
1153 C Invalid input parameter.
1154 C
1155     6 IER = 1
1156       RETURN
1157 C
1158 C No cells contain a point within RMAX of P, or
1159 C   SW = 0 and thus D .GE. RW(K) for all K.
1160 C
1161     7 C = 0.
1162       CX = 0.
1163       CY = 0.
1164       CXX = 0.
1165       CXY = 0.
1166       CYY = 0.
1167       IER = 2
1168       RETURN
1169       END
1170       SUBROUTINE GETNP2 (PX,PY,X,Y,NR,LCELL,LNEXT,XMIN,YMIN,
1171      .                   DX,DY, NP,DSQ)
1172       INTEGER NR, LCELL(NR,NR), LNEXT(*), NP
1173       DOUBLE PRECISION PX, PY, X(*), Y(*), XMIN, YMIN, DX,
1174      .                 DY, DSQ
1175 C
1176 C***********************************************************
1177 C
1178 C                                               From CSHEP2D
1179 C                                            Robert J. Renka
1180 C                                  Dept. of Computer Science
1181 C                                       Univ. of North Texas
1182 C                                           renka@cs.unt.edu
1183 C                                                   02/03/97
1184 C
1185 C   Given a set of N nodes and the data structure defined in
1186 C Subroutine STORE2, this subroutine uses the cell method to
1187 C find the closest unmarked node NP to a specified point P.
1188 C NP is then marked by setting LNEXT(NP) to -LNEXT(NP).  (A
1189 C node is marked if and only if the corresponding LNEXT ele-
1190 C ment is negative.  The absolute values of LNEXT elements,
1191 C however, must be preserved.)  Thus, the closest M nodes to
1192 C P may be determined by a sequence of M calls to this rou-
1193 C tine.  Note that if the nearest neighbor to node K is to
1194 C be determined (PX = X(K) and PY = Y(K)), then K should be
1195 C marked before the call to this routine.
1196 C
1197 C   The search is begun in the cell containing (or closest
1198 C to) P and proceeds outward in rectangular layers until all
1199 C cells which contain points within distance R of P have
1200 C been searched, where R is the distance from P to the first
1201 C unmarked node encountered (infinite if no unmarked nodes
1202 C are present).
1203 C
1204 C   This code is essentially unaltered from the subroutine
1205 C of the same name in QSHEP2D.
1206 C
1207 C On input:
1208 C
1209 C       PX,PY = Cartesian coordinates of the point P whose
1210 C               nearest unmarked neighbor is to be found.
1211 C
1212 C       X,Y = Arrays of length N, for N .GE. 2, containing
1213 C             the Cartesian coordinates of the nodes.
1214 C
1215 C       NR = Number of rows and columns in the cell grid.
1216 C            Refer to Subroutine STORE2.  NR .GE. 1.
1217 C
1218 C       LCELL = NR by NR array of nodal indexes associated
1219 C               with cells.  Refer to Subroutine STORE2.
1220 C
1221 C       LNEXT = Array of length N containing next-node
1222 C               indexes (or their negatives).  Refer to
1223 C               Subroutine STORE2.
1224 C
1225 C       XMIN,YMIN,DX,DY = Minimum nodal coordinates and cell
1226 C                         dimensions.  DX and DY must be
1227 C                         positive.  Refer to Subroutine
1228 C                         STORE2.
1229 C
1230 C   Input parameters other than LNEXT are not altered by
1231 C this routine.  With the exception of (PX,PY) and the signs
1232 C of LNEXT elements, these parameters should be unaltered
1233 C from their values on output from Subroutine STORE2.
1234 C
1235 C On output:
1236 C
1237 C       NP = Index (for X and Y) of the nearest unmarked
1238 C            node to P, or 0 if all nodes are marked or NR
1239 C            .LT. 1 or DX .LE. 0 or DY .LE. 0.  LNEXT(NP)
1240 C            .LT. 0 IF NP .NE. 0.
1241 C
1242 C       DSQ = Squared Euclidean distance between P and node
1243 C             NP, or 0 if NP = 0.
1244 C
1245 C Modules required by GETNP2:  None
1246 C
1247 C Intrinsic functions called by GETNP2:  ABS, INT, SQRT
1248 C
1249 C***********************************************************
1250 C
1251       INTEGER I, I0, I1, I2, IMAX, IMIN, J, J0, J1, J2,
1252      .        JMAX, JMIN, L, LMIN, LN
1253       LOGICAL FIRST
1254       DOUBLE PRECISION DELX, DELY, R, RSMIN, RSQ, XP, YP
1255 C
1256 C Local parameters:
1257 C
1258 C DELX,DELY =   PX-XMIN, PY-YMIN
1259 C FIRST =       Logical variable with value TRUE iff the
1260 C                 first unmarked node has yet to be
1261 C                 encountered
1262 C I,J =         Cell indexes in the range [I1,I2] X [J1,J2]
1263 C I0,J0 =       Indexes of the cell containing or closest
1264 C                 to P
1265 C I1,I2,J1,J2 = Range of cell indexes defining the layer
1266 C                 whose intersection with the range
1267 C                 [IMIN,IMAX] X [JMIN,JMAX] is currently
1268 C                 being searched
1269 C IMIN,IMAX =   Cell row indexes defining the range of the
1270 C                 search
1271 C JMIN,JMAX =   Cell column indexes defining the range of
1272 C                 the search
1273 C L,LN =        Indexes of nodes in cell (I,J)
1274 C LMIN =        Current candidate for NP
1275 C R =           Distance from P to node LMIN
1276 C RSMIN =       Squared distance from P to node LMIN
1277 C RSQ =         Squared distance from P to node L
1278 C XP,YP =       Local copy of PX,PY -- coordinates of P
1279 C
1280       XP = PX
1281       YP = PY
1282 C
1283 C Test for invalid input parameters.
1284 C
1285       IF (NR .LT. 1  .OR.  DX .LE. 0.  .OR.  DY .LE. 0.)
1286      .  GO TO 9
1287 C
1288 C Initialize parameters.
1289 C
1290       FIRST = .TRUE.
1291       IMIN = 1
1292       IMAX = NR
1293       JMIN = 1
1294       JMAX = NR
1295       DELX = XP - XMIN
1296       DELY = YP - YMIN
1297       I0 = INT(DELX/DX) + 1
1298       IF (I0 .LT. 1) I0 = 1
1299       IF (I0 .GT. NR) I0 = NR
1300       J0 = INT(DELY/DY) + 1
1301       IF (J0 .LT. 1) J0 = 1
1302       IF (J0 .GT. NR) J0 = NR
1303       I1 = I0
1304       I2 = I0
1305       J1 = J0
1306       J2 = J0
1307 C
1308 C Outer loop on layers, inner loop on layer cells, excluding
1309 C   those outside the range [IMIN,IMAX] X [JMIN,JMAX].
1310 C
1311     1 DO 6 J = J1,J2
1312         IF (J .GT. JMAX) GO TO 7
1313         IF (J .LT. JMIN) GO TO 6
1314         DO 5 I = I1,I2
1315           IF (I .GT. IMAX) GO TO 6
1316           IF (I .LT. IMIN) GO TO 5
1317           IF (J .NE. J1  .AND.  J .NE. J2  .AND.  I .NE. I1
1318      .        .AND.  I .NE. I2) GO TO 5
1319 C
1320 C Search cell (I,J) for unmarked nodes L.
1321 C
1322           L = LCELL(I,J)
1323           IF (L .EQ. 0) GO TO 5
1324 C
1325 C   Loop on nodes in cell (I,J).
1326 C
1327     2     LN = LNEXT(L)
1328           IF (LN .LT. 0) GO TO 4
1329 C
1330 C   Node L is not marked.
1331 C
1332           RSQ = (X(L)-XP)**2 + (Y(L)-YP)**2
1333           IF (.NOT. FIRST) GO TO 3
1334 C
1335 C   Node L is the first unmarked neighbor of P encountered.
1336 C     Initialize LMIN to the current candidate for NP, and
1337 C     RSMIN to the squared distance from P to LMIN.  IMIN,
1338 C     IMAX, JMIN, and JMAX are updated to define the smal-
1339 C     lest rectangle containing a circle of radius R =
1340 C     Sqrt(RSMIN) centered at P, and contained in [1,NR] X
1341 C     [1,NR] (except that, if P is outside the rectangle
1342 C     defined by the nodes, it is possible that IMIN > NR,
1343 C     IMAX < 1, JMIN > NR, or JMAX < 1).  FIRST is reset to
1344 C     FALSE.
1345 C
1346           LMIN = L
1347           RSMIN = RSQ
1348           R = SQRT(RSMIN)
1349           IMIN = INT((DELX-R)/DX) + 1
1350           IF (IMIN .LT. 1) IMIN = 1
1351           IMAX = INT((DELX+R)/DX) + 1
1352           IF (IMAX .GT. NR) IMAX = NR
1353           JMIN = INT((DELY-R)/DY) + 1
1354           IF (JMIN .LT. 1) JMIN = 1
1355           JMAX = INT((DELY+R)/DY) + 1
1356           IF (JMAX .GT. NR) JMAX = NR
1357           FIRST = .FALSE.
1358           GO TO 4
1359 C
1360 C   Test for node L closer than LMIN to P.
1361 C
1362     3     IF (RSQ .GE. RSMIN) GO TO 4
1363 C
1364 C   Update LMIN and RSMIN.
1365 C
1366           LMIN = L
1367           RSMIN = RSQ
1368 C
1369 C   Test for termination of loop on nodes in cell (I,J).
1370 C
1371     4     IF (ABS(LN) .EQ. L) GO TO 5
1372           L = ABS(LN)
1373           GO TO 2
1374     5     CONTINUE
1375     6   CONTINUE
1376 C
1377 C Test for termination of loop on cell layers.
1378 C
1379     7 IF (I1 .LE. IMIN  .AND.  I2 .GE. IMAX  .AND.
1380      .    J1 .LE. JMIN  .AND.  J2 .GE. JMAX) GO TO 8
1381       I1 = I1 - 1
1382       I2 = I2 + 1
1383       J1 = J1 - 1
1384       J2 = J2 + 1
1385       GO TO 1
1386 C
1387 C Unless no unmarked nodes were encountered, LMIN is the
1388 C   closest unmarked node to P.
1389 C
1390     8 IF (FIRST) GO TO 9
1391       NP = LMIN
1392       DSQ = RSMIN
1393       LNEXT(LMIN) = -LNEXT(LMIN)
1394       RETURN
1395 C
1396 C Error:  NR, DX, or DY is invalid or all nodes are marked.
1397 C
1398     9 NP = 0
1399       DSQ = 0.
1400       RETURN
1401       END
1402       SUBROUTINE GIVENS ( A,B, C,S)
1403       DOUBLE PRECISION A, B, C, S
1404 C
1405 C***********************************************************
1406 C
1407 C                                               From SRFPACK
1408 C                                            Robert J. Renka
1409 C                                  Dept. of Computer Science
1410 C                                       Univ. of North Texas
1411 C                                           renka@cs.unt.edu
1412 C                                                   09/01/88
1413 C
1414 C   This subroutine constructs the Givens plane rotation,
1415 C
1416 C           ( C  S)
1417 C       G = (     ) , where C*C + S*S = 1,
1418 C           (-S  C)
1419 C
1420 C which zeros the second component of the vector (A,B)**T
1421 C (transposed).  Subroutine ROTATE may be called to apply
1422 C the transformation to a 2 by N matrix.
1423 C
1424 C   This routine is identical to subroutine SROTG from the
1425 C LINPACK BLAS (Basic Linear Algebra Subroutines).
1426 C
1427 C On input:
1428 C
1429 C       A,B = Components of the vector defining the rota-
1430 C             tion.  These are overwritten by values R
1431 C             and Z (described below) which define C and S.
1432 C
1433 C On output:
1434 C
1435 C       A = Signed Euclidean norm R of the input vector:
1436 C           R = +/-SQRT(A*A + B*B)
1437 C
1438 C       B = Value Z such that:
1439 C             C = SQRT(1-Z*Z) and S=Z if ABS(Z) .LE. 1, and
1440 C             C = 1/Z and S = SQRT(1-C*C) if ABS(Z) > 1.
1441 C
1442 C       C = +/-(A/R) or 1 if R = 0.
1443 C
1444 C       S = +/-(B/R) or 0 if R = 0.
1445 C
1446 C Modules required by GIVENS:  None
1447 C
1448 C Intrinsic functions called by GIVENS:  ABS, SQRT
1449 C
1450 C***********************************************************
1451 C
1452       DOUBLE PRECISION AA, BB, R, U, V
1453 C
1454 C Local parameters:
1455 C
1456 C AA,BB = Local copies of A and B
1457 C R =     C*A + S*B = +/-SQRT(A*A+B*B)
1458 C U,V =   Variables used to scale A and B for computing R
1459 C
1460       AA = A
1461       BB = B
1462       IF (ABS(AA) .LE. ABS(BB)) GO TO 1
1463 C
1464 C ABS(A) > ABS(B).
1465 C
1466       U = AA + AA
1467       V = BB/U
1468       R = SQRT(.25 + V*V) * U
1469       C = AA/R
1470       S = V * (C + C)
1471 C
1472 C Note that R has the sign of A, C > 0, and S has
1473 C   SIGN(A)*SIGN(B).
1474 C
1475       B = S
1476       A = R
1477       RETURN
1478 C
1479 C ABS(A) .LE. ABS(B).
1480 C
1481     1 IF (BB .EQ. 0.) GO TO 2
1482       U = BB + BB
1483       V = AA/U
1484 C
1485 C Store R in A.
1486 C
1487       A = SQRT(.25 + V*V) * U
1488       S = BB/A
1489       C = V * (S + S)
1490 C
1491 C Note that R has the sign of B, S > 0, and C has
1492 C   SIGN(A)*SIGN(B).
1493 C
1494       B = 1.
1495       IF (C .NE. 0.) B = 1./C
1496       RETURN
1497 C
1498 C A = B = 0.
1499 C
1500     2 C = 1.
1501       S = 0.
1502       RETURN
1503       END
1504       SUBROUTINE ROTATE (N,C,S, X,Y )
1505       INTEGER N
1506       DOUBLE PRECISION C, S, X(N), Y(N)
1507 C
1508 C***********************************************************
1509 C
1510 C                                               From SRFPACK
1511 C                                            Robert J. Renka
1512 C                                  Dept. of Computer Science
1513 C                                       Univ. of North Texas
1514 C                                           renka@cs.unt.edu
1515 C                                                   09/01/88
1516 C
1517 C                                                ( C  S)
1518 C   This subroutine applies the Givens rotation  (     )  to
1519 C                                                (-S  C)
1520 C                    (X(1) ... X(N))
1521 C the 2 by N matrix  (             ) .
1522 C                    (Y(1) ... Y(N))
1523 C
1524 C   This routine is identical to subroutine SROT from the
1525 C LINPACK BLAS (Basic Linear Algebra Subroutines).
1526 C
1527 C On input:
1528 C
1529 C       N = Number of columns to be rotated.
1530 C
1531 C       C,S = Elements of the Givens rotation.  Refer to
1532 C             subroutine GIVENS.
1533 C
1534 C The above parameters are not altered by this routine.
1535 C
1536 C       X,Y = Arrays of length .GE. N containing the compo-
1537 C             nents of the vectors to be rotated.
1538 C
1539 C On output:
1540 C
1541 C       X,Y = Arrays containing the rotated vectors (not
1542 C             altered if N < 1).
1543 C
1544 C Modules required by ROTATE:  None
1545 C
1546 C***********************************************************
1547 C
1548       INTEGER I
1549       DOUBLE PRECISION XI, YI
1550 C
1551       DO 1 I = 1,N
1552         XI = X(I)
1553         YI = Y(I)
1554         X(I) = C*XI + S*YI
1555         Y(I) = -S*XI + C*YI
1556     1   CONTINUE
1557       RETURN
1558       END
1559       SUBROUTINE SETUP2 (XK,YK,ZK,XI,YI,ZI,S1,S2,S3,R, ROW)
1560       DOUBLE PRECISION XK, YK, ZK, XI, YI, ZI, S1, S2, S3,
1561      .                 R, ROW(10)
1562 C
1563 C***********************************************************
1564 C
1565 C                                               From CSHEP2D
1566 C                                            Robert J. Renka
1567 C                                  Dept. of Computer Science
1568 C                                       Univ. of North Texas
1569 C                                           renka@cs.unt.edu
1570 C                                                   02/03/97
1571 C
1572 C   This subroutine sets up the I-th row of an augmented re-
1573 C gression matrix for a weighted least squares fit of a
1574 C cubic function f(x,y) to a set of data values z, where
1575 C f(XK,YK) = ZK.  The first four columns (cubic terms) are
1576 C scaled by S3, the next three columns (quadratic terms)
1577 C are scaled by S2, and the eighth and ninth columns (lin-
1578 C ear terms) are scaled by S1.
1579 C
1580 C On input:
1581 C
1582 C       XK,YK = Coordinates of node K.
1583 C
1584 C       ZK = Data value at node K to be interpolated by f.
1585 C
1586 C       XI,YI,ZI = Coordinates and data value at node I.
1587 C
1588 C       S1,S2,S3 = Scale factors.
1589 C
1590 C       R = Radius of influence about node K defining the
1591 C           weight.
1592 C
1593 C The above parameters are not altered by this routine.
1594 C
1595 C       ROW = Array of length 10.
1596 C
1597 C On output:
1598 C
1599 C       ROW = Array containing a row of the augmented re-
1600 C             gression matrix.
1601 C
1602 C Modules required by SETUP2:  None
1603 C
1604 C Intrinsic function called by SETUP2:  SQRT
1605 C
1606 C***********************************************************
1607 C
1608       INTEGER I
1609       DOUBLE PRECISION D, DX, DXSQ, DY, DYSQ, W, W1, W2, W3
1610 C
1611 C Local parameters:
1612 C
1613 C D =    Distance between nodes K and I
1614 C DX =   XI - XK
1615 C DXSQ = DX*DX
1616 C DY =   YI - YK
1617 C DYSQ = DY*DY
1618 C I =    DO-loop index
1619 C W =    Weight associated with the row:  (R-D)/(R*D)
1620 C          (0 if D = 0 or D > R)
1621 C W1 =   S1*W
1622 C W2 =   S2*W
1623 C W3 =   W3*W
1624 C
1625       DX = XI - XK
1626       DY = YI - YK
1627       DXSQ = DX*DX
1628       DYSQ = DY*DY
1629       D = SQRT(DXSQ + DYSQ)
1630       IF (D .LE. 0.  .OR.  D .GE. R) GO TO 1
1631       W = (R-D)/R/D
1632       W1 = S1*W
1633       W2 = S2*W
1634       W3 = S3*W
1635       ROW(1) = DXSQ*DX*W3
1636       ROW(2) = DXSQ*DY*W3
1637       ROW(3) = DX*DYSQ*W3
1638       ROW(4) = DYSQ*DY*W3
1639       ROW(5) = DXSQ*W2
1640       ROW(6) = DX*DY*W2
1641       ROW(7) = DYSQ*W2
1642       ROW(8) = DX*W1
1643       ROW(9) = DY*W1
1644       ROW(10) = (ZI - ZK)*W
1645       RETURN
1646 C
1647 C Nodes K and I coincide or node I is outside of the radius
1648 C   of influence.  Set ROW to the zero vector.
1649 C
1650     1 DO 2 I = 1,10
1651         ROW(I) = 0.
1652     2   CONTINUE
1653       RETURN
1654       END
1655       SUBROUTINE STORE2 (N,X,Y,NR, LCELL,LNEXT,XMIN,YMIN,DX,
1656      .                   DY,IER)
1657       INTEGER N, NR, LCELL(NR,NR), LNEXT(N), IER
1658       DOUBLE PRECISION X(N), Y(N), XMIN, YMIN, DX, DY
1659 C
1660 C***********************************************************
1661 C
1662 C                                               From CSHEP2D
1663 C                                            Robert J. Renka
1664 C                                  Dept. of Computer Science
1665 C                                       Univ. of North Texas
1666 C                                           renka@cs.unt.edu
1667 C                                                   03/28/97
1668 C
1669 C   Given a set of N arbitrarily distributed nodes in the
1670 C plane, this subroutine creates a data structure for a
1671 C cell-based method of solving closest-point problems.  The
1672 C smallest rectangle containing the nodes is partitioned
1673 C into an NR by NR uniform grid of cells, and nodes are as-
1674 C sociated with cells.  In particular, the data structure
1675 C stores the indexes of the nodes contained in each cell.
1676 C For a uniform random distribution of nodes, the nearest
1677 C node to an arbitrary point can be determined in constant
1678 C expected time.
1679 C
1680 C   This code is essentially unaltered from the subroutine
1681 C of the same name in QSHEP2D.
1682 C
1683 C On input:
1684 C
1685 C       N = Number of nodes.  N .GE. 2.
1686 C
1687 C       X,Y = Arrays of length N containing the Cartesian
1688 C             coordinates of the nodes.
1689 C
1690 C       NR = Number of rows and columns in the grid.  The
1691 C            cell density (average number of nodes per cell)
1692 C            is D = N/(NR**2).  A recommended value, based
1693 C            on empirical evidence, is D = 3 -- NR =
1694 C            Sqrt(N/3).  NR .GE. 1.
1695 C
1696 C The above parameters are not altered by this routine.
1697 C
1698 C       LCELL = Array of length .GE. NR**2.
1699 C
1700 C       LNEXT = Array of length .GE. N.
1701 C
1702 C On output:
1703 C
1704 C       LCELL = NR by NR cell array such that LCELL(I,J)
1705 C               contains the index (for X and Y) of the
1706 C               first node (node with smallest index) in
1707 C               cell (I,J), or LCELL(I,J) = 0 if no nodes
1708 C               are contained in the cell.  The upper right
1709 C               corner of cell (I,J) has coordinates (XMIN+
1710 C               I*DX,YMIN+J*DY).  LCELL is not defined if
1711 C               IER .NE. 0.
1712 C
1713 C       LNEXT = Array of next-node indexes such that
1714 C               LNEXT(K) contains the index of the next node
1715 C               in the cell which contains node K, or
1716 C               LNEXT(K) = K if K is the last node in the
1717 C               cell for K = 1,...,N.  (The nodes contained
1718 C               in a cell are ordered by their indexes.)
1719 C               If, for example, cell (I,J) contains nodes
1720 C               2, 3, and 5 (and no others), then LCELL(I,J)
1721 C               = 2, LNEXT(2) = 3, LNEXT(3) = 5, and
1722 C               LNEXT(5) = 5.  LNEXT is not defined if
1723 C               IER .NE. 0.
1724 C
1725 C       XMIN,YMIN = Cartesian coordinates of the lower left
1726 C                   corner of the rectangle defined by the
1727 C                   nodes (smallest nodal coordinates) un-
1728 C                   less IER = 1.  The upper right corner is
1729 C                   (XMAX,YMAX) for XMAX = XMIN + NR*DX and
1730 C                   YMAX = YMIN + NR*DY.
1731 C
1732 C       DX,DY = Dimensions of the cells unless IER = 1.  DX
1733 C               = (XMAX-XMIN)/NR and DY = (YMAX-YMIN)/NR,
1734 C               where XMIN, XMAX, YMIN, and YMAX are the
1735 C               extrema of X and Y.
1736 C
1737 C       IER = Error indicator:
1738 C             IER = 0 if no errors were encountered.
1739 C             IER = 1 if N < 2 or NR < 1.
1740 C             IER = 2 if DX = 0 or DY = 0.
1741 C
1742 C Modules required by STORE2:  None
1743 C
1744 C Intrinsic functions called by STORE2:  DBLE, INT
1745 C
1746 C***********************************************************
1747 C
1748       INTEGER I, J, K, L, NN, NNR
1749       DOUBLE PRECISION DELX, DELY, XMN, XMX, YMN, YMX
1750 C
1751 C Local parameters:
1752 C
1753 C DELX,DELY = Components of the cell dimensions -- local
1754 C               copies of DX,DY
1755 C I,J =       Cell indexes
1756 C K =         Nodal index
1757 C L =         Index of a node in cell (I,J)
1758 C NN =        Local copy of N
1759 C NNR =       Local copy of NR
1760 C XMN,XMX =   Range of nodal X coordinates
1761 C YMN,YMX =   Range of nodal Y coordinates
1762 C
1763       NN = N
1764       NNR = NR
1765       IF (NN .LT. 2  .OR.  NNR .LT. 1) GO TO 5
1766 C
1767 C Compute the dimensions of the rectangle containing the
1768 C   nodes.
1769 C
1770       XMN = X(1)
1771       XMX = XMN
1772       YMN = Y(1)
1773       YMX = YMN
1774       DO 1 K = 2,NN
1775         IF (X(K) .LT. XMN) XMN = X(K)
1776         IF (X(K) .GT. XMX) XMX = X(K)
1777         IF (Y(K) .LT. YMN) YMN = Y(K)
1778         IF (Y(K) .GT. YMX) YMX = Y(K)
1779     1   CONTINUE
1780       XMIN = XMN
1781       YMIN = YMN
1782 C
1783 C Compute cell dimensions and test for zero area.
1784 C
1785       DELX = (XMX-XMN)/DBLE(NNR)
1786       DELY = (YMX-YMN)/DBLE(NNR)
1787       DX = DELX
1788       DY = DELY
1789       IF (DELX .EQ. 0.  .OR.  DELY .EQ. 0.) GO TO 6
1790 C
1791 C Initialize LCELL.
1792 C
1793       DO 3 J = 1,NNR
1794         DO 2 I = 1,NNR
1795           LCELL(I,J) = 0
1796     2     CONTINUE
1797     3   CONTINUE
1798 C
1799 C Loop on nodes, storing indexes in LCELL and LNEXT.
1800 C
1801       DO 4 K = NN,1,-1
1802         I = INT((X(K)-XMN)/DELX) + 1
1803         IF (I .GT. NNR) I = NNR
1804         J = INT((Y(K)-YMN)/DELY) + 1
1805         IF (J .GT. NNR) J = NNR
1806         L = LCELL(I,J)
1807         LNEXT(K) = L
1808         IF (L .EQ. 0) LNEXT(K) = K
1809         LCELL(I,J) = K
1810     4   CONTINUE
1811 C
1812 C No errors encountered.
1813 C
1814       IER = 0
1815       RETURN
1816 C
1817 C Invalid input parameter.
1818 C
1819     5 IER = 1
1820       RETURN
1821 C
1822 C DX = 0 or DY = 0.
1823 C
1824     6 IER = 2
1825       RETURN
1826       END