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