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