Bug 13268 fixed: There was a missing return in C++ code (ScilabObjet::getArgumentId)
[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 C
562       CS2VAL = 0.
563 C
564       IF (N .LT. 10  .OR.  NR .LT. 1  .OR.  DX .LE. 0.  .OR.
565      .    DY .LE. 0.  .OR.  RMAX .LT. 0.) RETURN
566 C
567 C Set IMIN, IMAX, JMIN, and JMAX to cell indexes defining
568 C   the range of the search for nodes whose radii include
569 C   P.  The cells which must be searched are those inter-
570 C   sected by (or contained in) a circle of radius RMAX
571 C   centered at P.
572 C
573       IMIN = INT((XP-XMIN-RMAX)/DX) + 1
574       IMAX = INT((XP-XMIN+RMAX)/DX) + 1
575       IF (IMIN .LT. 1) IMIN = 1
576       IF (IMAX .GT. NR) IMAX = NR
577       JMIN = INT((YP-YMIN-RMAX)/DY) + 1
578       JMAX = INT((YP-YMIN+RMAX)/DY) + 1
579       IF (JMIN .LT. 1) JMIN = 1
580       IF (JMAX .GT. NR) JMAX = NR
581 C
582 C The following is a test for no cells within the circle
583 C   of radius RMAX.
584 C
585       IF (IMIN .GT. IMAX  .OR.  JMIN .GT. JMAX) GO TO 6
586 C
587 C Accumulate weight values in SW and weighted nodal function
588 C   values in SWC.  The weights are W(K) = ((R-D)+/(R*D))**3
589 C   for R = RW(K) and D = distance between P and node K.
590 C
591       SW = 0.
592       SWC = 0.
593 C
594 C Outer loop on cells (I,J).
595 C
596       DO 4 J = JMIN,JMAX
597         DO 3 I = IMIN,IMAX
598           K = LCELL(I,J)
599           IF (K .EQ. 0) GO TO 3
600 C
601 C Inner loop on nodes K.
602 C
603     1     DELX = XP - X(K)
604           DELY = YP - Y(K)
605           D = SQRT(DELX*DELX + DELY*DELY)
606           R = RW(K)
607           IF (D .GE. R) GO TO 2
608           IF (D .EQ. 0.) GO TO 5
609           W = (1.0/D - 1.0/R)**3
610           SW = SW + W
611           SWC = SWC + W*( ( (A(1,K)*DELX+A(2,K)*DELY+
612      .                       A(5,K))*DELX + (A(3,K)*DELY+
613      .                       A(6,K))*DELY + A(8,K) )*DELX +
614      .                    ( (A(4,K)*DELY+A(7,K))*DELY +
615      .                      A(9,K) )*DELY + F(K) )
616 C
617 C Bottom of loop on nodes in cell (I,J).
618 C
619     2     KP = K
620           K = LNEXT(KP)
621           IF (K .NE. KP) GO TO 1
622     3     CONTINUE
623     4   CONTINUE
624 C
625 C SW = 0 iff P is not within the radius R(K) for any node K.
626 C
627       IF (SW .EQ. 0.) GO TO 6
628       CS2VAL = SWC/SW
629       RETURN
630 C
631 C (PX,PY) = (X(K),Y(K)).
632 C
633     5 CS2VAL = F(K)
634       RETURN
635 C
636 C All weights are 0 at P.
637 C
638     6 CS2VAL = 0.
639       RETURN
640       END
641
642       SUBROUTINE CS2GRD (PX,PY,N,X,Y,F,NR,LCELL,LNEXT,XMIN,
643      .                   YMIN,DX,DY,RMAX,RW,A, C,CX,CY,IER)
644       INTEGER N, NR, LCELL(NR,NR), LNEXT(N), IER
645       DOUBLE PRECISION PX, PY, X(N), Y(N), F(N), XMIN, YMIN,
646      .                 DX, DY, RMAX, RW(N), A(9,N), C, CX,
647      .                 CY
648 C
649 C***********************************************************
650 C
651 C                                               From CSHEP2D
652 C                                            Robert J. Renka
653 C                                  Dept. of Computer Science
654 C                                       Univ. of North Texas
655 C                                           renka@cs.unt.edu
656 C                                                   02/03/97
657 C
658 C   This subroutine computes the value and gradient at P =
659 C (PX,PY) of the interpolatory function C defined in Sub-
660 C routine CSHEP2.  C is a weighted sum of cubic nodal
661 C functions.
662 C
663 C On input:
664 C
665 C       PX,PY = Cartesian coordinates of the point P at
666 C               which C and its partial derivatives are
667 C               to be evaluated.
668 C
669 C       N = Number of nodes and data values defining C.
670 C           N .GE. 10.
671 C
672 C       X,Y,F = Arrays of length N containing the nodes and
673 C               data values interpolated by C.
674 C
675 C       NR = Number of rows and columns in the cell grid.
676 C            Refer to Subroutine STORE2.  NR .GE. 1.
677 C
678 C       LCELL = NR by NR array of nodal indexes associated
679 C               with cells.  Refer to Subroutine STORE2.
680 C
681 C       LNEXT = Array of length N containing next-node
682 C               indexes.  Refer to Subroutine STORE2.
683 C
684 C       XMIN,YMIN,DX,DY = Minimum nodal coordinates and cell
685 C                         dimensions.  DX and DY must be
686 C                         positive.  Refer to Subroutine
687 C                         STORE2.
688 C
689 C       RMAX = Largest element in RW -- maximum radius R(k).
690 C
691 C       RW = Array of length N containing the radii R(k)
692 C            which enter into the weights W(k) defining C.
693 C
694 C       A = 9 by N array containing the coefficients for
695 C           cubic nodal function C(k) in column k.
696 C
697 C   Input parameters are not altered by this subroutine.
698 C The parameters other than PX and PY should be input
699 C unaltered from their values on output from CSHEP2.  This
700 C subroutine should not be called if a nonzero error flag
701 C was returned by CSHEP2.
702 C
703 C On output:
704 C
705 C       C = Value of C at (PX,PY) unless IER .EQ. 1, in
706 C           which case no values are returned.
707 C
708 C       CX,CY = First partial derivatives of C at (PX,PY)
709 C               unless IER .EQ. 1.
710 C
711 C       IER = Error indicator:
712 C             IER = 0 if no errors were encountered.
713 C             IER = 1 if N, NR, DX, DY or RMAX is invalid.
714 C             IER = 2 if no errors were encountered but
715 C                     (PX,PY) is not within the radius R(k)
716 C                     for any node k (and thus C=CX=CY=0).
717 C
718 C Modules required by CS2GRD:  None
719 C
720 C Intrinsic functions called by CS2GRD:  INT, SQRT
721 C
722 C***********************************************************
723 C
724       INTEGER I, IMAX, IMIN, J, JMAX, JMIN, K, KP
725       DOUBLE PRECISION CK, CKX, CKY, D, DELX, DELY, R, SW,
726      .                 SWC, SWCX, SWCY, SWS, SWX, SWY, T, W,
727      .                 WX, WY, XP, YP
728 C
729 C Local parameters:
730 C
731 C CK =        Value of cubic nodal function C(K) at P
732 C CKX,CKY =   Partial derivatives of C(K) with respect to X
733 C               and Y, respectively
734 C D =         Distance between P and node K
735 C DELX =      XP - X(K)
736 C DELY =      YP - Y(K)
737 C I =         Cell row index in the range IMIN to IMAX
738 C IMIN,IMAX = Range of cell row indexes of the cells
739 C               intersected by a disk of radius RMAX
740 C               centered at P
741 C J =         Cell column index in the range JMIN to JMAX
742 C JMIN,JMAX = Range of cell column indexes of the cells
743 C               intersected by a disk of radius RMAX
744 C               centered at P
745 C K =         Index of a node in cell (I,J)
746 C KP =        Previous value of K in the sequence of nodes
747 C               in cell (I,J)
748 C R =         Radius of influence for node K
749 C SW =        Sum of weights W(K)
750 C SWC =       Sum of weighted nodal function values at P
751 C SWCX,SWCY = Partial derivatives of SWC with respect to X
752 C               and Y, respectively
753 C SWS =       SW**2
754 C SWX,SWY =   Partial derivatives of SW with respect to X
755 C               and Y, respectively
756 C T =         Temporary variable
757 C W =         Weight W(K) value at P:  ((R-D)+/(R*D))**3,
758 C               where (R-D)+ = 0 if R < D
759 C WX,WY =     Partial derivatives of W with respect to X
760 C               and Y, respectively
761 C XP,YP =     Local copies of PX and PY -- coordinates of P
762 C
763       XP = PX
764       YP = PY
765       IF (N .LT. 10  .OR.  NR .LT. 1  .OR.  DX .LE. 0.  .OR.
766      .    DY .LE. 0.  .OR.  RMAX .LT. 0.) GO TO 6
767 C
768 C Set IMIN, IMAX, JMIN, and JMAX to cell indexes defining
769 C   the range of the search for nodes whose radii include
770 C   P.  The cells which must be searched are those inter-
771 C   sected by (or contained in) a circle of radius RMAX
772 C   centered at P.
773 C
774       IMIN = INT((XP-XMIN-RMAX)/DX) + 1
775       IMAX = INT((XP-XMIN+RMAX)/DX) + 1
776       IF (IMIN .LT. 1) IMIN = 1
777       IF (IMAX .GT. NR) IMAX = NR
778       JMIN = INT((YP-YMIN-RMAX)/DY) + 1
779       JMAX = INT((YP-YMIN+RMAX)/DY) + 1
780       IF (JMIN .LT. 1) JMIN = 1
781       IF (JMAX .GT. NR) JMAX = NR
782 C
783 C The following is a test for no cells within the circle
784 C   of radius RMAX.
785 C
786       IF (IMIN .GT. IMAX  .OR.  JMIN .GT. JMAX) GO TO 7
787 C
788 C C = SWC/SW = Sum(W(K)*C(K))/Sum(W(K)), where the sum is
789 C   from K = 1 to N, C(K) is the cubic nodal function value,
790 C   and W(K) = ((R-D)+/(R*D))**3 for radius R(K) and dist-
791 C   ance D(K).  Thus
792 C
793 C        CX = (SWCX*SW - SWC*SWX)/SW**2  and
794 C        CY = (SWCY*SW - SWC*SWY)/SW**2
795 C
796 C   where SWCX and SWX are partial derivatives with respect
797 C   to X of SWC and SW, respectively.  SWCY and SWY are de-
798 C   fined similarly.
799 C
800       SW = 0.
801       SWX = 0.
802       SWY = 0.
803       SWC = 0.
804       SWCX = 0.
805       SWCY = 0.
806 C
807 C Outer loop on cells (I,J).
808 C
809       DO 4 J = JMIN,JMAX
810         DO 3 I = IMIN,IMAX
811           K = LCELL(I,J)
812           IF (K .EQ. 0) GO TO 3
813 C
814 C Inner loop on nodes K.
815 C
816     1     DELX = XP - X(K)
817           DELY = YP - Y(K)
818           D = SQRT(DELX*DELX + DELY*DELY)
819           R = RW(K)
820           IF (D .GE. R) GO TO 2
821           IF (D .EQ. 0.) GO TO 5
822           T = (1.0/D - 1.0/R)
823           W = T**3
824           T = -3.0*T*T/(D**3)
825           WX = DELX*T
826           WY = DELY*T
827           T = A(2,K)*DELX + A(3,K)*DELY + A(6,K)
828           CKY = ( 3.0*A(4,K)*DELY + A(3,K)*DELX +
829      .            2.0*A(7,K) )*DELY + T*DELX + A(9,K)
830           T = T*DELY + A(8,K)
831           CKX = ( 3.0*A(1,K)*DELX + A(2,K)*DELY +
832      .            2.0*A(5,K) )*DELX + T
833           CK = ( (A(1,K)*DELX+A(5,K))*DELX + T )*DELX +
834      .         ( (A(4,K)*DELY+A(7,K))*DELY + A(9,K) )*DELY +
835      .         F(K)
836           SW = SW + W
837           SWX = SWX + WX
838           SWY = SWY + WY
839           SWC = SWC + W*CK
840           SWCX = SWCX + WX*CK + W*CKX
841           SWCY = SWCY + WY*CK + W*CKY
842 C
843 C Bottom of loop on nodes in cell (I,J).
844 C
845     2     KP = K
846           K = LNEXT(KP)
847           IF (K .NE. KP) GO TO 1
848     3     CONTINUE
849     4   CONTINUE
850 C
851 C SW = 0 iff P is not within the radius R(K) for any node K.
852 C
853       IF (SW .EQ. 0.) GO TO 7
854       C = SWC/SW
855       SWS = SW*SW
856       CX = (SWCX*SW - SWC*SWX)/SWS
857       CY = (SWCY*SW - SWC*SWY)/SWS
858       IER = 0
859       RETURN
860 C
861 C (PX,PY) = (X(K),Y(K)).
862 C
863     5 C = F(K)
864       CX = A(8,K)
865       CY = A(9,K)
866       IER = 0
867       RETURN
868 C
869 C Invalid input parameter.
870 C
871     6 IER = 1
872       RETURN
873 C
874 C No cells contain a point within RMAX of P, or
875 C   SW = 0 and thus D .GE. RW(K) for all K.
876 C
877     7 C = 0.
878       CX = 0.
879       CY = 0.
880       IER = 2
881       RETURN
882       END
883       SUBROUTINE CS2HES (PX,PY,N,X,Y,F,NR,LCELL,LNEXT,XMIN,
884      .                   YMIN,DX,DY,RMAX,RW,A, C,CX,CY,CXX,
885      .                   CXY,CYY,IER)
886       INTEGER N, NR, LCELL(NR,NR), LNEXT(N), IER
887       DOUBLE PRECISION PX, PY, X(N), Y(N), F(N), XMIN, YMIN,
888      .                 DX, DY, RMAX, RW(N), A(9,N), C, CX,
889      .                 CY, CXX, CXY, CYY
890 C
891 C***********************************************************
892 C
893 C                                               From CSHEP2D
894 C                                            Robert J. Renka
895 C                                  Dept. of Computer Science
896 C                                       Univ. of North Texas
897 C                                           renka@cs.unt.edu
898 C                                                   02/03/97
899 C
900 C   This subroutine computes the value, gradient, and
901 C Hessian at P = (PX,PY) of the interpolatory function C
902 C defined in Subroutine CSHEP2.  C is a weighted sum of
903 C cubic nodal functions.
904 C
905 C On input:
906 C
907 C       PX,PY = Cartesian coordinates of the point P at
908 C               which C and its partial derivatives are
909 C               to be evaluated.
910 C
911 C       N = Number of nodes and data values defining C.
912 C           N .GE. 10.
913 C
914 C       X,Y,F = Arrays of length N containing the nodes and
915 C               data values interpolated by C.
916 C
917 C       NR = Number of rows and columns in the cell grid.
918 C            Refer to Subroutine STORE2.  NR .GE. 1.
919 C
920 C       LCELL = NR by NR array of nodal indexes associated
921 C               with cells.  Refer to Subroutine STORE2.
922 C
923 C       LNEXT = Array of length N containing next-node
924 C               indexes.  Refer to Subroutine STORE2.
925 C
926 C       XMIN,YMIN,DX,DY = Minimum nodal coordinates and cell
927 C                         dimensions.  DX and DY must be
928 C                         positive.  Refer to Subroutine
929 C                         STORE2.
930 C
931 C       RMAX = Largest element in RW -- maximum radius R(k).
932 C
933 C       RW = Array of length N containing the radii R(k)
934 C            which enter into the weights W(k) defining C.
935 C
936 C       A = 9 by N array containing the coefficients for
937 C           cubic nodal function C(k) in column k.
938 C
939 C   Input parameters are not altered by this subroutine.
940 C The parameters other than PX and PY should be input
941 C unaltered from their values on output from CSHEP2.  This
942 C subroutine should not be called if a nonzero error flag
943 C was returned by CSHEP2.
944 C
945 C On output:
946 C
947 C       C = Value of C at (PX,PY) unless IER .EQ. 1, in
948 C           which case no values are returned.
949 C
950 C       CX,CY = First partial derivatives of C at (PX,PY)
951 C               unless IER .EQ. 1.
952 C
953 C       CXX,CXY,CYY = Second partial derivatives of C at
954 C                     (PX,PY) unless IER .EQ. 1.
955 C
956 C       IER = Error indicator:
957 C             IER = 0 if no errors were encountered.
958 C             IER = 1 if N, NR, DX, DY or RMAX is invalid.
959 C             IER = 2 if no errors were encountered but
960 C                     (PX,PY) is not within the radius R(k)
961 C                     for any node k (and thus C = 0).
962 C
963 C Modules required by CS2HES:  None
964 C
965 C Intrinsic functions called by CS2HES:  INT, SQRT
966 C
967 C***********************************************************
968 C
969       INTEGER I, IMAX, IMIN, J, JMAX, JMIN, K, KP
970       DOUBLE PRECISION CK, CKX, CKXX, CKXY, CKY, CKYY, D,
971      .                 DELX, DELY, DXSQ, DYSQ, R, SW, SWC,
972      .                 SWCX, SWCXX, SWCXY, SWCY, SWCYY, SWS,
973      .                 SWX, SWXX, SWXY, SWY, SWYY, T1, T2,
974      .                 T3, T4, W, WX, WXX, WXY, WY, WYY, XP,
975      .                 YP
976 C
977 C Local parameters:
978 C
979 C CK =        Value of cubic nodal function C(K) at P
980 C CKX,CKY =   Partial derivatives of C(K) with respect to X
981 C               and Y, respectively
982 C CKXX,CKXY,CKYY = Second partial derivatives of CK
983 C D =         Distance between P and node K
984 C DELX =      XP - X(K)
985 C DELY =      YP - Y(K)
986 C DXSQ,DYSQ = DELX**2, DELY**2
987 C I =         Cell row index in the range IMIN to IMAX
988 C IMIN,IMAX = Range of cell row indexes of the cells
989 C               intersected by a disk of radius RMAX
990 C               centered at P
991 C J =         Cell column index in the range JMIN to JMAX
992 C JMIN,JMAX = Range of cell column indexes of the cells
993 C               intersected by a disk of radius RMAX
994 C               centered at P
995 C K =         Index of a node in cell (I,J)
996 C KP =        Previous value of K in the sequence of nodes
997 C               in cell (I,J)
998 C R =         Radius of influence for node K
999 C SW =        Sum of weights W(K)
1000 C SWC =       Sum of weighted nodal function values at P
1001 C SWCX,SWCY = Partial derivatives of SWC with respect to X
1002 C               and Y, respectively
1003 C SWCXX,SWCXY,SWCYY = Second partial derivatives of SWC
1004 C SWS =       SW**2
1005 C SWX,SWY =   Partial derivatives of SW with respect to X
1006 C               and Y, respectively
1007 C SWXX,SWXY,SWYY = Second partial derivatives of SW
1008 C T1,T2,T3,T4 = Temporary variables
1009 C W =         Weight W(K) value at P:  ((R-D)+/(R*D))**3,
1010 C               where (R-D)+ = 0 if R < D
1011 C WX,WY =     Partial derivatives of W with respect to X
1012 C               and Y, respectively
1013 C WXX,WXY,WYY = Second partial derivatives of W
1014 C XP,YP =     Local copies of PX and PY -- coordinates of P
1015 C
1016       XP = PX
1017       YP = PY
1018       IF (N .LT. 10  .OR.  NR .LT. 1  .OR.  DX .LE. 0.  .OR.
1019      .    DY .LE. 0.  .OR.  RMAX .LT. 0.) GO TO 6
1020 C
1021 C Set IMIN, IMAX, JMIN, and JMAX to cell indexes defining
1022 C   the range of the search for nodes whose radii include
1023 C   P.  The cells which must be searched are those inter-
1024 C   sected by (or contained in) a circle of radius RMAX
1025 C   centered at P.
1026 C
1027       IMIN = INT((XP-XMIN-RMAX)/DX) + 1
1028       IMAX = INT((XP-XMIN+RMAX)/DX) + 1
1029       IF (IMIN .LT. 1) IMIN = 1
1030       IF (IMAX .GT. NR) IMAX = NR
1031       JMIN = INT((YP-YMIN-RMAX)/DY) + 1
1032       JMAX = INT((YP-YMIN+RMAX)/DY) + 1
1033       IF (JMIN .LT. 1) JMIN = 1
1034       IF (JMAX .GT. NR) JMAX = NR
1035 C
1036 C The following is a test for no cells within the circle
1037 C   of radius RMAX.
1038 C
1039       IF (IMIN .GT. IMAX  .OR.  JMIN .GT. JMAX) GO TO 7
1040 C
1041 C C = SWC/SW = Sum(W(K)*C(K))/Sum(W(K)), where the sum is
1042 C   from K = 1 to N, C(K) is the cubic nodal function value,
1043 C   and W(K) = ((R-D)+/(R*D))**3 for radius R(K) and dist-
1044 C   ance D(K).  Thus
1045 C
1046 C        CX = (SWCX*SW - SWC*SWX)/SW**2  and
1047 C        CY = (SWCY*SW - SWC*SWY)/SW**2
1048 C
1049 C   where SWCX and SWX are partial derivatives with respect
1050 C   to x of SWC and SW, respectively.  SWCY and SWY are de-
1051 C   fined similarly.  The second partials are
1052 C
1053 C        CXX = ( SW*(SWCXX -    2*SWX*CX) - SWC*SWXX )/SW**2
1054 C        CXY = ( SW*(SWCXY-SWX*CY-SWY*CX) - SWC*SWXY )/SW**2
1055 C        CYY = ( SW*(SWCYY -    2*SWY*CY) - SWC*SWYY )/SW**2
1056 C
1057 C   where SWCXX and SWXX are second partials with respect
1058 C   to x, SWCXY and SWXY are mixed partials, and SWCYY and
1059 C   SWYY are second partials with respect to y.
1060 C
1061       SW = 0.
1062       SWX = 0.
1063       SWY = 0.
1064       SWXX = 0.
1065       SWXY = 0.
1066       SWYY = 0.
1067       SWC = 0.
1068       SWCX = 0.
1069       SWCY = 0.
1070       SWCXX = 0.
1071       SWCXY = 0.
1072       SWCYY = 0.
1073 C
1074 C Outer loop on cells (I,J).
1075 C
1076       DO 4 J = JMIN,JMAX
1077         DO 3 I = IMIN,IMAX
1078           K = LCELL(I,J)
1079           IF (K .EQ. 0) GO TO 3
1080 C
1081 C Inner loop on nodes K.
1082 C
1083     1     DELX = XP - X(K)
1084           DELY = YP - Y(K)
1085           DXSQ = DELX*DELX
1086           DYSQ = DELY*DELY
1087           D = SQRT(DXSQ + DYSQ)
1088           R = RW(K)
1089           IF (D .GE. R) GO TO 2
1090           IF (D .EQ. 0.) GO TO 5
1091           T1 = (1.0/D - 1.0/R)
1092           W = T1**3
1093           T2 = -3.0*T1*T1/(D**3)
1094           WX = DELX*T2
1095           WY = DELY*T2
1096           T1 = 3.0*T1*(2.0+3.0*D*T1)/(D**6)
1097           WXX = T1*DXSQ + T2
1098           WXY = T1*DELX*DELY
1099           WYY = T1*DYSQ + T2
1100           T1 = A(1,K)*DELX + A(2,K)*DELY + A(5,K)
1101           T2 = T1 + T1 + A(1,K)*DELX
1102           T3 = A(4,K)*DELY + A(3,K)*DELX + A(7,K)
1103           T4 = T3 + T3 + A(4,K)*DELY
1104           CK = (T1*DELX + A(6,K)*DELY + A(8,K))*DELX +
1105      .         (T3*DELY + A(9,K))*DELY + F(K)
1106           CKX = T2*DELX + (A(3,K)*DELY+A(6,K))*DELY + A(8,K)
1107           CKY = T4*DELY + (A(2,K)*DELX+A(6,K))*DELX + A(9,K)
1108           CKXX = T2 + 3.0*A(1,K)*DELX
1109           CKXY = 2.0*(A(2,K)*DELX + A(3,K)*DELY) + A(6,K)
1110           CKYY = T4 + 3.0*A(4,K)*DELY
1111           SW = SW + W
1112           SWX = SWX + WX
1113           SWY = SWY + WY
1114           SWXX = SWXX + WXX
1115           SWXY = SWXY + WXY
1116           SWYY = SWYY + WYY
1117           SWC = SWC + W*CK
1118           SWCX = SWCX + WX*CK + W*CKX
1119           SWCY = SWCY + WY*CK + W*CKY
1120           SWCXX = SWCXX + W*CKXX + 2.0*WX*CKX + CK*WXX
1121           SWCXY = SWCXY + W*CKXY + WX*CKY + WY*CKX + CK*WXY
1122           SWCYY = SWCYY + W*CKYY + 2.0*WY*CKY + CK*WYY
1123 C
1124 C Bottom of loop on nodes in cell (I,J).
1125 C
1126     2     KP = K
1127           K = LNEXT(KP)
1128           IF (K .NE. KP) GO TO 1
1129     3     CONTINUE
1130     4   CONTINUE
1131 C
1132 C SW = 0 iff P is not within the radius R(K) for any node K.
1133 C
1134       IF (SW .EQ. 0.) GO TO 7
1135       C = SWC/SW
1136       SWS = SW*SW
1137       CX = (SWCX*SW - SWC*SWX)/SWS
1138       CY = (SWCY*SW - SWC*SWY)/SWS
1139       CXX = (SW*(SWCXX-2.0*SWX*CX) - SWC*SWXX)/SWS
1140       CXY = (SW*(SWCXY-SWY*CX-SWX*CY) - SWC*SWXY)/SWS
1141       CYY = (SW*(SWCYY-2.0*SWY*CY) - SWC*SWYY)/SWS
1142       IER = 0
1143       RETURN
1144 C
1145 C (PX,PY) = (X(K),Y(K)).
1146 C
1147     5 C = F(K)
1148       CX = A(8,K)
1149       CY = A(9,K)
1150       CXX = 2.0*A(5,K)
1151       CXY = A(6,K)
1152       CYY = 2.0*A(7,K)
1153       IER = 0
1154       RETURN
1155 C
1156 C Invalid input parameter.
1157 C
1158     6 IER = 1
1159       RETURN
1160 C
1161 C No cells contain a point within RMAX of P, or
1162 C   SW = 0 and thus D .GE. RW(K) for all K.
1163 C
1164     7 C = 0.
1165       CX = 0.
1166       CY = 0.
1167       CXX = 0.
1168       CXY = 0.
1169       CYY = 0.
1170       IER = 2
1171       RETURN
1172       END
1173       SUBROUTINE GETNP2 (PX,PY,X,Y,NR,LCELL,LNEXT,XMIN,YMIN,
1174      .                   DX,DY, NP,DSQ)
1175       INTEGER NR, LCELL(NR,NR), LNEXT(*), NP
1176       DOUBLE PRECISION PX, PY, X(*), Y(*), XMIN, YMIN, DX,
1177      .                 DY, DSQ
1178 C
1179 C***********************************************************
1180 C
1181 C                                               From CSHEP2D
1182 C                                            Robert J. Renka
1183 C                                  Dept. of Computer Science
1184 C                                       Univ. of North Texas
1185 C                                           renka@cs.unt.edu
1186 C                                                   02/03/97
1187 C
1188 C   Given a set of N nodes and the data structure defined in
1189 C Subroutine STORE2, this subroutine uses the cell method to
1190 C find the closest unmarked node NP to a specified point P.
1191 C NP is then marked by setting LNEXT(NP) to -LNEXT(NP).  (A
1192 C node is marked if and only if the corresponding LNEXT ele-
1193 C ment is negative.  The absolute values of LNEXT elements,
1194 C however, must be preserved.)  Thus, the closest M nodes to
1195 C P may be determined by a sequence of M calls to this rou-
1196 C tine.  Note that if the nearest neighbor to node K is to
1197 C be determined (PX = X(K) and PY = Y(K)), then K should be
1198 C marked before the call to this routine.
1199 C
1200 C   The search is begun in the cell containing (or closest
1201 C to) P and proceeds outward in rectangular layers until all
1202 C cells which contain points within distance R of P have
1203 C been searched, where R is the distance from P to the first
1204 C unmarked node encountered (infinite if no unmarked nodes
1205 C are present).
1206 C
1207 C   This code is essentially unaltered from the subroutine
1208 C of the same name in QSHEP2D.
1209 C
1210 C On input:
1211 C
1212 C       PX,PY = Cartesian coordinates of the point P whose
1213 C               nearest unmarked neighbor is to be found.
1214 C
1215 C       X,Y = Arrays of length N, for N .GE. 2, containing
1216 C             the Cartesian coordinates of the nodes.
1217 C
1218 C       NR = Number of rows and columns in the cell grid.
1219 C            Refer to Subroutine STORE2.  NR .GE. 1.
1220 C
1221 C       LCELL = NR by NR array of nodal indexes associated
1222 C               with cells.  Refer to Subroutine STORE2.
1223 C
1224 C       LNEXT = Array of length N containing next-node
1225 C               indexes (or their negatives).  Refer to
1226 C               Subroutine STORE2.
1227 C
1228 C       XMIN,YMIN,DX,DY = Minimum nodal coordinates and cell
1229 C                         dimensions.  DX and DY must be
1230 C                         positive.  Refer to Subroutine
1231 C                         STORE2.
1232 C
1233 C   Input parameters other than LNEXT are not altered by
1234 C this routine.  With the exception of (PX,PY) and the signs
1235 C of LNEXT elements, these parameters should be unaltered
1236 C from their values on output from Subroutine STORE2.
1237 C
1238 C On output:
1239 C
1240 C       NP = Index (for X and Y) of the nearest unmarked
1241 C            node to P, or 0 if all nodes are marked or NR
1242 C            .LT. 1 or DX .LE. 0 or DY .LE. 0.  LNEXT(NP)
1243 C            .LT. 0 IF NP .NE. 0.
1244 C
1245 C       DSQ = Squared Euclidean distance between P and node
1246 C             NP, or 0 if NP = 0.
1247 C
1248 C Modules required by GETNP2:  None
1249 C
1250 C Intrinsic functions called by GETNP2:  ABS, INT, SQRT
1251 C
1252 C***********************************************************
1253 C
1254       INTEGER I, I0, I1, I2, IMAX, IMIN, J, J0, J1, J2,
1255      .        JMAX, JMIN, L, LMIN, LN
1256       LOGICAL FIRST
1257       DOUBLE PRECISION DELX, DELY, R, RSMIN, RSQ, XP, YP
1258 C
1259 C Local parameters:
1260 C
1261 C DELX,DELY =   PX-XMIN, PY-YMIN
1262 C FIRST =       Logical variable with value TRUE iff the
1263 C                 first unmarked node has yet to be
1264 C                 encountered
1265 C I,J =         Cell indexes in the range [I1,I2] X [J1,J2]
1266 C I0,J0 =       Indexes of the cell containing or closest
1267 C                 to P
1268 C I1,I2,J1,J2 = Range of cell indexes defining the layer
1269 C                 whose intersection with the range
1270 C                 [IMIN,IMAX] X [JMIN,JMAX] is currently
1271 C                 being searched
1272 C IMIN,IMAX =   Cell row indexes defining the range of the
1273 C                 search
1274 C JMIN,JMAX =   Cell column indexes defining the range of
1275 C                 the search
1276 C L,LN =        Indexes of nodes in cell (I,J)
1277 C LMIN =        Current candidate for NP
1278 C R =           Distance from P to node LMIN
1279 C RSMIN =       Squared distance from P to node LMIN
1280 C RSQ =         Squared distance from P to node L
1281 C XP,YP =       Local copy of PX,PY -- coordinates of P
1282 C
1283       XP = PX
1284       YP = PY
1285 C
1286 C Test for invalid input parameters.
1287 C
1288       IF (NR .LT. 1  .OR.  DX .LE. 0.  .OR.  DY .LE. 0.)
1289      .  GO TO 9
1290 C
1291 C Initialize parameters.
1292 C
1293       FIRST = .TRUE.
1294       IMIN = 1
1295       IMAX = NR
1296       JMIN = 1
1297       JMAX = NR
1298       DELX = XP - XMIN
1299       DELY = YP - YMIN
1300       I0 = INT(DELX/DX) + 1
1301       IF (I0 .LT. 1) I0 = 1
1302       IF (I0 .GT. NR) I0 = NR
1303       J0 = INT(DELY/DY) + 1
1304       IF (J0 .LT. 1) J0 = 1
1305       IF (J0 .GT. NR) J0 = NR
1306       I1 = I0
1307       I2 = I0
1308       J1 = J0
1309       J2 = J0
1310 C
1311 C Outer loop on layers, inner loop on layer cells, excluding
1312 C   those outside the range [IMIN,IMAX] X [JMIN,JMAX].
1313 C
1314     1 DO 6 J = J1,J2
1315         IF (J .GT. JMAX) GO TO 7
1316         IF (J .LT. JMIN) GO TO 6
1317         DO 5 I = I1,I2
1318           IF (I .GT. IMAX) GO TO 6
1319           IF (I .LT. IMIN) GO TO 5
1320           IF (J .NE. J1  .AND.  J .NE. J2  .AND.  I .NE. I1
1321      .        .AND.  I .NE. I2) GO TO 5
1322 C
1323 C Search cell (I,J) for unmarked nodes L.
1324 C
1325           L = LCELL(I,J)
1326           IF (L .EQ. 0) GO TO 5
1327 C
1328 C   Loop on nodes in cell (I,J).
1329 C
1330     2     LN = LNEXT(L)
1331           IF (LN .LT. 0) GO TO 4
1332 C
1333 C   Node L is not marked.
1334 C
1335           RSQ = (X(L)-XP)**2 + (Y(L)-YP)**2
1336           IF (.NOT. FIRST) GO TO 3
1337 C
1338 C   Node L is the first unmarked neighbor of P encountered.
1339 C     Initialize LMIN to the current candidate for NP, and
1340 C     RSMIN to the squared distance from P to LMIN.  IMIN,
1341 C     IMAX, JMIN, and JMAX are updated to define the smal-
1342 C     lest rectangle containing a circle of radius R =
1343 C     Sqrt(RSMIN) centered at P, and contained in [1,NR] X
1344 C     [1,NR] (except that, if P is outside the rectangle
1345 C     defined by the nodes, it is possible that IMIN > NR,
1346 C     IMAX < 1, JMIN > NR, or JMAX < 1).  FIRST is reset to
1347 C     FALSE.
1348 C
1349           LMIN = L
1350           RSMIN = RSQ
1351           R = SQRT(RSMIN)
1352           IMIN = INT((DELX-R)/DX) + 1
1353           IF (IMIN .LT. 1) IMIN = 1
1354           IMAX = INT((DELX+R)/DX) + 1
1355           IF (IMAX .GT. NR) IMAX = NR
1356           JMIN = INT((DELY-R)/DY) + 1
1357           IF (JMIN .LT. 1) JMIN = 1
1358           JMAX = INT((DELY+R)/DY) + 1
1359           IF (JMAX .GT. NR) JMAX = NR
1360           FIRST = .FALSE.
1361           GO TO 4
1362 C
1363 C   Test for node L closer than LMIN to P.
1364 C
1365     3     IF (RSQ .GE. RSMIN) GO TO 4
1366 C
1367 C   Update LMIN and RSMIN.
1368 C
1369           LMIN = L
1370           RSMIN = RSQ
1371 C
1372 C   Test for termination of loop on nodes in cell (I,J).
1373 C
1374     4     IF (ABS(LN) .EQ. L) GO TO 5
1375           L = ABS(LN)
1376           GO TO 2
1377     5     CONTINUE
1378     6   CONTINUE
1379 C
1380 C Test for termination of loop on cell layers.
1381 C
1382     7 IF (I1 .LE. IMIN  .AND.  I2 .GE. IMAX  .AND.
1383      .    J1 .LE. JMIN  .AND.  J2 .GE. JMAX) GO TO 8
1384       I1 = I1 - 1
1385       I2 = I2 + 1
1386       J1 = J1 - 1
1387       J2 = J2 + 1
1388       GO TO 1
1389 C
1390 C Unless no unmarked nodes were encountered, LMIN is the
1391 C   closest unmarked node to P.
1392 C
1393     8 IF (FIRST) GO TO 9
1394       NP = LMIN
1395       DSQ = RSMIN
1396       LNEXT(LMIN) = -LNEXT(LMIN)
1397       RETURN
1398 C
1399 C Error:  NR, DX, or DY is invalid or all nodes are marked.
1400 C
1401     9 NP = 0
1402       DSQ = 0.
1403       RETURN
1404       END
1405       SUBROUTINE GIVENS ( A,B, C,S)
1406       DOUBLE PRECISION A, B, C, S
1407 C
1408 C***********************************************************
1409 C
1410 C                                               From SRFPACK
1411 C                                            Robert J. Renka
1412 C                                  Dept. of Computer Science
1413 C                                       Univ. of North Texas
1414 C                                           renka@cs.unt.edu
1415 C                                                   09/01/88
1416 C
1417 C   This subroutine constructs the Givens plane rotation,
1418 C
1419 C           ( C  S)
1420 C       G = (     ) , where C*C + S*S = 1,
1421 C           (-S  C)
1422 C
1423 C which zeros the second component of the vector (A,B)**T
1424 C (transposed).  Subroutine ROTATE may be called to apply
1425 C the transformation to a 2 by N matrix.
1426 C
1427 C   This routine is identical to subroutine SROTG from the
1428 C LINPACK BLAS (Basic Linear Algebra Subroutines).
1429 C
1430 C On input:
1431 C
1432 C       A,B = Components of the vector defining the rota-
1433 C             tion.  These are overwritten by values R
1434 C             and Z (described below) which define C and S.
1435 C
1436 C On output:
1437 C
1438 C       A = Signed Euclidean norm R of the input vector:
1439 C           R = +/-SQRT(A*A + B*B)
1440 C
1441 C       B = Value Z such that:
1442 C             C = SQRT(1-Z*Z) and S=Z if ABS(Z) .LE. 1, and
1443 C             C = 1/Z and S = SQRT(1-C*C) if ABS(Z) > 1.
1444 C
1445 C       C = +/-(A/R) or 1 if R = 0.
1446 C
1447 C       S = +/-(B/R) or 0 if R = 0.
1448 C
1449 C Modules required by GIVENS:  None
1450 C
1451 C Intrinsic functions called by GIVENS:  ABS, SQRT
1452 C
1453 C***********************************************************
1454 C
1455       DOUBLE PRECISION AA, BB, R, U, V
1456 C
1457 C Local parameters:
1458 C
1459 C AA,BB = Local copies of A and B
1460 C R =     C*A + S*B = +/-SQRT(A*A+B*B)
1461 C U,V =   Variables used to scale A and B for computing R
1462 C
1463       AA = A
1464       BB = B
1465       IF (ABS(AA) .LE. ABS(BB)) GO TO 1
1466 C
1467 C ABS(A) > ABS(B).
1468 C
1469       U = AA + AA
1470       V = BB/U
1471       R = SQRT(.25 + V*V) * U
1472       C = AA/R
1473       S = V * (C + C)
1474 C
1475 C Note that R has the sign of A, C > 0, and S has
1476 C   SIGN(A)*SIGN(B).
1477 C
1478       B = S
1479       A = R
1480       RETURN
1481 C
1482 C ABS(A) .LE. ABS(B).
1483 C
1484     1 IF (BB .EQ. 0.) GO TO 2
1485       U = BB + BB
1486       V = AA/U
1487 C
1488 C Store R in A.
1489 C
1490       A = SQRT(.25 + V*V) * U
1491       S = BB/A
1492       C = V * (S + S)
1493 C
1494 C Note that R has the sign of B, S > 0, and C has
1495 C   SIGN(A)*SIGN(B).
1496 C
1497       B = 1.
1498       IF (C .NE. 0.) B = 1./C
1499       RETURN
1500 C
1501 C A = B = 0.
1502 C
1503     2 C = 1.
1504       S = 0.
1505       RETURN
1506       END
1507       SUBROUTINE ROTATE (N,C,S, X,Y )
1508       INTEGER N
1509       DOUBLE PRECISION C, S, X(N), Y(N)
1510 C
1511 C***********************************************************
1512 C
1513 C                                               From SRFPACK
1514 C                                            Robert J. Renka
1515 C                                  Dept. of Computer Science
1516 C                                       Univ. of North Texas
1517 C                                           renka@cs.unt.edu
1518 C                                                   09/01/88
1519 C
1520 C                                                ( C  S)
1521 C   This subroutine applies the Givens rotation  (     )  to
1522 C                                                (-S  C)
1523 C                    (X(1) ... X(N))
1524 C the 2 by N matrix  (             ) .
1525 C                    (Y(1) ... Y(N))
1526 C
1527 C   This routine is identical to subroutine SROT from the
1528 C LINPACK BLAS (Basic Linear Algebra Subroutines).
1529 C
1530 C On input:
1531 C
1532 C       N = Number of columns to be rotated.
1533 C
1534 C       C,S = Elements of the Givens rotation.  Refer to
1535 C             subroutine GIVENS.
1536 C
1537 C The above parameters are not altered by this routine.
1538 C
1539 C       X,Y = Arrays of length .GE. N containing the compo-
1540 C             nents of the vectors to be rotated.
1541 C
1542 C On output:
1543 C
1544 C       X,Y = Arrays containing the rotated vectors (not
1545 C             altered if N < 1).
1546 C
1547 C Modules required by ROTATE:  None
1548 C
1549 C***********************************************************
1550 C
1551       INTEGER I
1552       DOUBLE PRECISION XI, YI
1553 C
1554       DO 1 I = 1,N
1555         XI = X(I)
1556         YI = Y(I)
1557         X(I) = C*XI + S*YI
1558         Y(I) = -S*XI + C*YI
1559     1   CONTINUE
1560       RETURN
1561       END
1562       SUBROUTINE SETUP2 (XK,YK,ZK,XI,YI,ZI,S1,S2,S3,R, ROW)
1563       DOUBLE PRECISION XK, YK, ZK, XI, YI, ZI, S1, S2, S3,
1564      .                 R, ROW(10)
1565 C
1566 C***********************************************************
1567 C
1568 C                                               From CSHEP2D
1569 C                                            Robert J. Renka
1570 C                                  Dept. of Computer Science
1571 C                                       Univ. of North Texas
1572 C                                           renka@cs.unt.edu
1573 C                                                   02/03/97
1574 C
1575 C   This subroutine sets up the I-th row of an augmented re-
1576 C gression matrix for a weighted least squares fit of a
1577 C cubic function f(x,y) to a set of data values z, where
1578 C f(XK,YK) = ZK.  The first four columns (cubic terms) are
1579 C scaled by S3, the next three columns (quadratic terms)
1580 C are scaled by S2, and the eighth and ninth columns (lin-
1581 C ear terms) are scaled by S1.
1582 C
1583 C On input:
1584 C
1585 C       XK,YK = Coordinates of node K.
1586 C
1587 C       ZK = Data value at node K to be interpolated by f.
1588 C
1589 C       XI,YI,ZI = Coordinates and data value at node I.
1590 C
1591 C       S1,S2,S3 = Scale factors.
1592 C
1593 C       R = Radius of influence about node K defining the
1594 C           weight.
1595 C
1596 C The above parameters are not altered by this routine.
1597 C
1598 C       ROW = Array of length 10.
1599 C
1600 C On output:
1601 C
1602 C       ROW = Array containing a row of the augmented re-
1603 C             gression matrix.
1604 C
1605 C Modules required by SETUP2:  None
1606 C
1607 C Intrinsic function called by SETUP2:  SQRT
1608 C
1609 C***********************************************************
1610 C
1611       INTEGER I
1612       DOUBLE PRECISION D, DX, DXSQ, DY, DYSQ, W, W1, W2, W3
1613 C
1614 C Local parameters:
1615 C
1616 C D =    Distance between nodes K and I
1617 C DX =   XI - XK
1618 C DXSQ = DX*DX
1619 C DY =   YI - YK
1620 C DYSQ = DY*DY
1621 C I =    DO-loop index
1622 C W =    Weight associated with the row:  (R-D)/(R*D)
1623 C          (0 if D = 0 or D > R)
1624 C W1 =   S1*W
1625 C W2 =   S2*W
1626 C W3 =   W3*W
1627 C
1628       DX = XI - XK
1629       DY = YI - YK
1630       DXSQ = DX*DX
1631       DYSQ = DY*DY
1632       D = SQRT(DXSQ + DYSQ)
1633       IF (D .LE. 0.  .OR.  D .GE. R) GO TO 1
1634       W = (R-D)/R/D
1635       W1 = S1*W
1636       W2 = S2*W
1637       W3 = S3*W
1638       ROW(1) = DXSQ*DX*W3
1639       ROW(2) = DXSQ*DY*W3
1640       ROW(3) = DX*DYSQ*W3
1641       ROW(4) = DYSQ*DY*W3
1642       ROW(5) = DXSQ*W2
1643       ROW(6) = DX*DY*W2
1644       ROW(7) = DYSQ*W2
1645       ROW(8) = DX*W1
1646       ROW(9) = DY*W1
1647       ROW(10) = (ZI - ZK)*W
1648       RETURN
1649 C
1650 C Nodes K and I coincide or node I is outside of the radius
1651 C   of influence.  Set ROW to the zero vector.
1652 C
1653     1 DO 2 I = 1,10
1654         ROW(I) = 0.
1655     2   CONTINUE
1656       RETURN
1657       END
1658       SUBROUTINE STORE2 (N,X,Y,NR, LCELL,LNEXT,XMIN,YMIN,DX,
1659      .                   DY,IER)
1660       INTEGER N, NR, LCELL(NR,NR), LNEXT(N), IER
1661       DOUBLE PRECISION X(N), Y(N), XMIN, YMIN, DX, DY
1662 C
1663 C***********************************************************
1664 C
1665 C                                               From CSHEP2D
1666 C                                            Robert J. Renka
1667 C                                  Dept. of Computer Science
1668 C                                       Univ. of North Texas
1669 C                                           renka@cs.unt.edu
1670 C                                                   03/28/97
1671 C
1672 C   Given a set of N arbitrarily distributed nodes in the
1673 C plane, this subroutine creates a data structure for a
1674 C cell-based method of solving closest-point problems.  The
1675 C smallest rectangle containing the nodes is partitioned
1676 C into an NR by NR uniform grid of cells, and nodes are as-
1677 C sociated with cells.  In particular, the data structure
1678 C stores the indexes of the nodes contained in each cell.
1679 C For a uniform random distribution of nodes, the nearest
1680 C node to an arbitrary point can be determined in constant
1681 C expected time.
1682 C
1683 C   This code is essentially unaltered from the subroutine
1684 C of the same name in QSHEP2D.
1685 C
1686 C On input:
1687 C
1688 C       N = Number of nodes.  N .GE. 2.
1689 C
1690 C       X,Y = Arrays of length N containing the Cartesian
1691 C             coordinates of the nodes.
1692 C
1693 C       NR = Number of rows and columns in the grid.  The
1694 C            cell density (average number of nodes per cell)
1695 C            is D = N/(NR**2).  A recommended value, based
1696 C            on empirical evidence, is D = 3 -- NR =
1697 C            Sqrt(N/3).  NR .GE. 1.
1698 C
1699 C The above parameters are not altered by this routine.
1700 C
1701 C       LCELL = Array of length .GE. NR**2.
1702 C
1703 C       LNEXT = Array of length .GE. N.
1704 C
1705 C On output:
1706 C
1707 C       LCELL = NR by NR cell array such that LCELL(I,J)
1708 C               contains the index (for X and Y) of the
1709 C               first node (node with smallest index) in
1710 C               cell (I,J), or LCELL(I,J) = 0 if no nodes
1711 C               are contained in the cell.  The upper right
1712 C               corner of cell (I,J) has coordinates (XMIN+
1713 C               I*DX,YMIN+J*DY).  LCELL is not defined if
1714 C               IER .NE. 0.
1715 C
1716 C       LNEXT = Array of next-node indexes such that
1717 C               LNEXT(K) contains the index of the next node
1718 C               in the cell which contains node K, or
1719 C               LNEXT(K) = K if K is the last node in the
1720 C               cell for K = 1,...,N.  (The nodes contained
1721 C               in a cell are ordered by their indexes.)
1722 C               If, for example, cell (I,J) contains nodes
1723 C               2, 3, and 5 (and no others), then LCELL(I,J)
1724 C               = 2, LNEXT(2) = 3, LNEXT(3) = 5, and
1725 C               LNEXT(5) = 5.  LNEXT is not defined if
1726 C               IER .NE. 0.
1727 C
1728 C       XMIN,YMIN = Cartesian coordinates of the lower left
1729 C                   corner of the rectangle defined by the
1730 C                   nodes (smallest nodal coordinates) un-
1731 C                   less IER = 1.  The upper right corner is
1732 C                   (XMAX,YMAX) for XMAX = XMIN + NR*DX and
1733 C                   YMAX = YMIN + NR*DY.
1734 C
1735 C       DX,DY = Dimensions of the cells unless IER = 1.  DX
1736 C               = (XMAX-XMIN)/NR and DY = (YMAX-YMIN)/NR,
1737 C               where XMIN, XMAX, YMIN, and YMAX are the
1738 C               extrema of X and Y.
1739 C
1740 C       IER = Error indicator:
1741 C             IER = 0 if no errors were encountered.
1742 C             IER = 1 if N < 2 or NR < 1.
1743 C             IER = 2 if DX = 0 or DY = 0.
1744 C
1745 C Modules required by STORE2:  None
1746 C
1747 C Intrinsic functions called by STORE2:  DBLE, INT
1748 C
1749 C***********************************************************
1750 C
1751       INTEGER I, J, K, L, NN, NNR
1752       DOUBLE PRECISION DELX, DELY, XMN, XMX, YMN, YMX
1753 C
1754 C Local parameters:
1755 C
1756 C DELX,DELY = Components of the cell dimensions -- local
1757 C               copies of DX,DY
1758 C I,J =       Cell indexes
1759 C K =         Nodal index
1760 C L =         Index of a node in cell (I,J)
1761 C NN =        Local copy of N
1762 C NNR =       Local copy of NR
1763 C XMN,XMX =   Range of nodal X coordinates
1764 C YMN,YMX =   Range of nodal Y coordinates
1765 C
1766       NN = N
1767       NNR = NR
1768       IF (NN .LT. 2  .OR.  NNR .LT. 1) GO TO 5
1769 C
1770 C Compute the dimensions of the rectangle containing the
1771 C   nodes.
1772 C
1773       XMN = X(1)
1774       XMX = XMN
1775       YMN = Y(1)
1776       YMX = YMN
1777       DO 1 K = 2,NN
1778         IF (X(K) .LT. XMN) XMN = X(K)
1779         IF (X(K) .GT. XMX) XMX = X(K)
1780         IF (Y(K) .LT. YMN) YMN = Y(K)
1781         IF (Y(K) .GT. YMX) YMX = Y(K)
1782     1   CONTINUE
1783       XMIN = XMN
1784       YMIN = YMN
1785 C
1786 C Compute cell dimensions and test for zero area.
1787 C
1788       DELX = (XMX-XMN)/DBLE(NNR)
1789       DELY = (YMX-YMN)/DBLE(NNR)
1790       DX = DELX
1791       DY = DELY
1792       IF (DELX .EQ. 0.  .OR.  DELY .EQ. 0.) GO TO 6
1793 C
1794 C Initialize LCELL.
1795 C
1796       DO 3 J = 1,NNR
1797         DO 2 I = 1,NNR
1798           LCELL(I,J) = 0
1799     2     CONTINUE
1800     3   CONTINUE
1801 C
1802 C Loop on nodes, storing indexes in LCELL and LNEXT.
1803 C
1804       DO 4 K = NN,1,-1
1805         I = INT((X(K)-XMN)/DELX) + 1
1806         IF (I .GT. NNR) I = NNR
1807         J = INT((Y(K)-YMN)/DELY) + 1
1808         IF (J .GT. NNR) J = NNR
1809         L = LCELL(I,J)
1810         LNEXT(K) = L
1811         IF (L .EQ. 0) LNEXT(K) = K
1812         LCELL(I,J) = K
1813     4   CONTINUE
1814 C
1815 C No errors encountered.
1816 C
1817       IER = 0
1818       RETURN
1819 C
1820 C Invalid input parameter.
1821 C
1822     5 IER = 1
1823       RETURN
1824 C
1825 C DX = 0 or DY = 0.
1826 C
1827     6 IER = 2
1828       RETURN
1829       END