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