C ALGORITHM 752, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 22, NO. 1, March, 1996, P. 9--17. C ######## including remark from renka (to appear) 4/dec/1998 #! /bin/sh # This is a shell archive, meaning: # 1. Remove everything above the #! /bin/sh line. # 2. Save the resulting text in a file. # 3. Execute the file with /bin/sh (not csh) to create the files: # Drivers/ # Drivers/Fortran77/ # Drivers/Fortran77/Sp/ # Drivers/Fortran77/Sp/RES # Drivers/Fortran77/Sp/RES.eps # Drivers/Fortran77/Sp/driver.f # Src/ # Src/Fortran77/ # Src/Fortran77/Sp/ # Src/Fortran77/Sp/srfpack.f # This archive created: Fri Dec 4 14:10:59 1998 export PATH; PATH=/bin:$PATH if test ! -d 'Drivers' then mkdir 'Drivers' fi cd 'Drivers' if test ! -d 'Fortran77' then mkdir 'Fortran77' fi cd 'Fortran77' if test ! -d 'Sp' then mkdir 'Sp' fi cd 'Sp' if test -f 'RES' then echo shar: will not over-write existing file "'RES'" else cat << SHAR_EOF > 'RES' SRFPACK Test: N = 12 Output from BNODES, AREA, and VOLUME BNODES: 4 Boundary Nodes, 29 Arcs, 18 Triangles AREAP: Area of convex hull = 0.10000000E+01 VOLUME: Area of convex hull = 0.10000001E+01 VOLUME: Area excluding constraint regions = 0.91000009E+00 GRADG Test: Maximum error = 0.00000000E+00 INTRC0 Test: Maximum error = 0.11920929E-06 GRADL Test: Maximum error = 0.35762787E-06 UNIF (INTRC1) Test: Maximum error = 0.11920929E-06 57 of 169 points required extrapolation ZGRADG/GETSIG Test: Maximum interpolation error = 0.51644452E-01 Tension Factors, N = 12 Nodes N1 N2 Tension 1 2 0.00000000 1 3 21.22009659 1 9 0.00000000 1 4 26.16072083 1 7 0.00000000 2 8 0.00000000 2 5 3.92235970 2 10 1.63794458 2 3 3.85027361 3 10 0.00000000 3 9 3.16620874 4 7 3.68543696 4 9 2.47128916 4 12 9.13391972 5 8 29.93745804 5 11 4.42266846 5 10 2.86108732 6 8 19.67565346 6 7 3.61773205 6 12 3.87078404 6 11 2.74047351 7 12 1.91911578 7 8 0.00000000 8 11 0.00000000 9 10 3.71387148 9 11 4.88137627 9 12 3.72744823 10 11 3.82832456 11 12 3.94437337 NA = 29 Arcs A contour plot file was successfully created. SHAR_EOF fi # end of overwriting check if test -f 'RES.eps' then echo shar: will not over-write existing file "'RES.eps'" else cat << SHAR_EOF > 'RES.eps' %!PS-Adobe-3.0 EPSF-3.0 %%BoundingBox: 36 126 576 666 %%Title: Contour Plot %%Creator: SRFPACK %%EndComments 36 126 moveto 36 666 lineto 576 666 lineto 576 126 lineto closepath stroke 81.000000 171.000000 translate 449.999969 449.999969 scale 0.004444 setlinewidth 0.386320 -0.100000 moveto 0.400000 -0.088527 lineto 0.500000 -0.006402 lineto 0.508153 0.000000 lineto 0.500000 0.006436 lineto 0.451171 0.100000 lineto 0.500000 0.153824 lineto 0.540608 0.200000 lineto 0.600000 0.254558 lineto 0.653935 0.300000 lineto 0.700000 0.346552 lineto 0.744069 0.400000 lineto 0.800000 0.459934 lineto 0.848164 0.500000 lineto 0.900000 0.546322 lineto 0.997466 0.500000 lineto 1.000000 0.497128 lineto 1.002238 0.500000 lineto 1.085377 0.600000 lineto 1.100000 0.617369 lineto stroke 0.984278 -0.100000 moveto 1.000000 -0.080803 lineto 1.083979 0.000000 lineto 1.100000 0.012349 lineto stroke 0.190768 -0.100000 moveto 0.200000 -0.091678 lineto 0.300000 -0.002483 lineto 0.302849 0.000000 lineto 0.336259 0.100000 lineto 0.400000 0.166105 lineto 0.431635 0.200000 lineto 0.500000 0.265461 lineto 0.536856 0.300000 lineto 0.600000 0.364824 lineto 0.634915 0.400000 lineto 0.700000 0.463103 lineto 0.734974 0.500000 lineto 0.800000 0.566544 lineto 0.836123 0.600000 lineto 0.900000 0.662533 lineto 1.000000 0.698944 lineto 1.000923 0.700000 lineto 1.090733 0.800000 lineto 1.100000 0.810226 lineto stroke 0.056183 -0.100000 moveto 0.100000 -0.059593 lineto 0.165080 0.000000 lineto 0.200000 0.040645 lineto 0.239653 0.100000 lineto 0.300000 0.160861 lineto 0.337880 0.200000 lineto 0.400000 0.261578 lineto 0.440059 0.300000 lineto 0.500000 0.360289 lineto 0.538760 0.400000 lineto 0.600000 0.461127 lineto 0.640009 0.500000 lineto 0.700000 0.558856 lineto 0.740232 0.600000 lineto 0.800000 0.660052 lineto 0.840721 0.700000 lineto 0.900000 0.759437 lineto 0.959318 0.800000 lineto 1.000000 0.835541 lineto 1.059281 0.900000 lineto 1.100000 0.943935 lineto stroke -0.051327 -0.100000 moveto 0.000000 -0.048781 lineto 0.052340 0.000000 lineto 0.100000 0.047556 lineto 0.146797 0.100000 lineto 0.200000 0.152988 lineto 0.246627 0.200000 lineto 0.300000 0.252960 lineto 0.347452 0.300000 lineto 0.400000 0.352002 lineto 0.447735 0.400000 lineto 0.500000 0.452087 lineto 0.547998 0.500000 lineto 0.600000 0.551935 lineto 0.648950 0.600000 lineto 0.700000 0.651098 lineto 0.748673 0.700000 lineto 0.800000 0.751868 lineto 0.848008 0.800000 lineto 0.900000 0.852658 lineto 0.952456 0.900000 lineto 1.000000 0.947816 lineto 1.048841 1.000000 lineto 1.100000 1.051011 lineto stroke 1.052475 1.100000 moveto 1.000000 1.047677 lineto 0.948833 1.000000 lineto 0.900000 0.950885 lineto 0.855033 0.900000 lineto 0.800000 0.844259 lineto 0.756972 0.800000 lineto 0.700000 0.742429 lineto 0.658040 0.700000 lineto 0.600000 0.641906 lineto 0.558829 0.600000 lineto 0.500000 0.541246 lineto 0.458762 0.500000 lineto 0.400000 0.441437 lineto 0.357481 0.400000 lineto 0.300000 0.343114 lineto 0.256143 0.300000 lineto 0.200000 0.244294 lineto 0.154804 0.200000 lineto 0.100000 0.145419 lineto 0.049027 0.100000 lineto 0.000000 0.051080 lineto -0.047789 0.000000 lineto -0.100000 -0.052101 lineto stroke 0.945271 1.100000 moveto 0.900000 1.058263 lineto 0.836391 1.000000 lineto 0.800000 0.957070 lineto 0.762168 0.900000 lineto 0.700000 0.835979 lineto 0.666163 0.800000 lineto 0.600000 0.733092 lineto 0.566403 0.700000 lineto 0.500000 0.632348 lineto 0.468857 0.600000 lineto 0.400000 0.531129 lineto 0.367659 0.500000 lineto 0.300000 0.434147 lineto 0.265898 0.400000 lineto 0.200000 0.335234 lineto 0.162738 0.300000 lineto 0.100000 0.238565 lineto 0.042833 0.200000 lineto 0.000000 0.163246 lineto -0.058178 0.100000 lineto -0.100000 0.054852 lineto stroke 0.811110 1.100000 moveto 0.800000 1.089981 lineto 0.700000 1.000358 lineto 0.699589 1.000000 lineto 0.665001 0.900000 lineto 0.600000 0.830319 lineto 0.572904 0.800000 lineto 0.500000 0.728250 lineto 0.471098 0.700000 lineto 0.400000 0.625919 lineto 0.374280 0.600000 lineto 0.300000 0.528705 lineto 0.271704 0.500000 lineto 0.200000 0.428198 lineto 0.168123 0.400000 lineto 0.100000 0.336071 lineto 0.001089 0.300000 lineto 0.000000 0.299449 lineto -0.089558 0.200000 lineto -0.100000 0.188480 lineto stroke 0.617194 1.100000 moveto 0.600000 1.085558 lineto 0.500000 1.002279 lineto 0.497100 1.000000 lineto 0.500000 0.997623 lineto 0.549325 0.900000 lineto 0.500000 0.842091 lineto 0.465213 0.800000 lineto 0.400000 0.736749 lineto 0.358553 0.700000 lineto 0.300000 0.640803 lineto 0.264301 0.600000 lineto 0.200000 0.534458 lineto 0.156367 0.500000 lineto 0.100000 0.452935 lineto 0.000000 0.499687 lineto -0.083727 0.400000 lineto -0.100000 0.380718 lineto stroke 0.012284 1.100000 moveto 0.000000 1.084869 lineto -0.087055 1.000000 lineto -0.100000 0.989896 lineto stroke showpage %%EOF  SHAR_EOF fi # end of overwriting check if test -f 'driver.f' then echo shar: will not over-write existing file "'driver.f'" else cat << SHAR_EOF > 'driver.f' C C C SRFTEST: Portable Test Driver for SRFPACK C 07/09/98 C C C This driver tests software package SRFPACK for construc- C ting an interpolatory surface to data values at scattered C points in the plane. The software modules are tested on C data for which the methods are exact, and maximum errors C are printed. These should be close to the machine pre- C cision. The final test using nonzero tension factors is C included primarily to demonstrate how the modules are C called. Given the small number of nodes, large interpo- C lation errors are expected in this case. C C By default, tests are performed on a simple data set C consisting of 12 nodes whose convex hull is the unit C square with a .3 by .3 square constraint region in the C center. However, by enabling the READ statements below C (C# in the first two columns), testing may be performed on C an arbitrary set of up to NMAX nodes with up to NCMAX C constraint curves. (Refer to the PARAMETER statements be- C low.) Note that the choice of interpolation/extrapolation C points assumes the nodes to approximately cover the unit C square. A data set consists of the following sequence of C records: C C N = Number of nodes (format I4): 3 to NMAX. C NCC = Number of constraint curves (format I4): 0 to C NCMAX. C (LCC(I), I = 1,NCC) = Indexes of the first node in each C constraint curve (format I4). 4 .LE. LCC(1) C and, for I .GT. 1, LCC(I-1) + 3 .LE. LCC(I) C .LE. N-2. (Each constraint curve has at least C three nodes.) C (X(I),Y(I), I = 1,N) = Nodal coordinates with non- C constraint nodes followed by the NCC C sequences of constraint nodes (format C 2F13.8). C C Various options may be changed by altering the DATA C statements below. C C This driver must be linked to TRIPACK and SRFPACK. C INTEGER I, IER, IFLGS, IMAX, IST, J, K, L, LIN, LNEW, . LOUT, LPLT, LW, LWK, MAXITG, MAXITZ, N, N6, . NA, NB, NCC, NCMAX, NCON, NI, NITG, NITZ, NJ, . NLS, NMAX, NT LOGICAL DFLAG, SFLAG REAL A1, A2, A3, DGMAX, DGMX, DMAX, DSM, DZMAX, . DZMX, PLTSIZ, PX, PY, PZ, SVAL, TOL, XF, YF, . ZTRUE, ZX, ZY REAL AREAP, FT, VOLUME C PARAMETER (NMAX=5000, NCMAX=5, N6=6*NMAX, LWK=2*NMAX, . NI=13, NJ=13) C C Array storage: C INTEGER LCC(NCMAX), LEND(NMAX), LIST(N6), LPTR(N6), . NODES(LWK) REAL GRAD(2,NMAX), P(NI), SIGMA(N6), . X(NMAX), Y(NMAX), Z(NMAX), ZZ(NI,NJ) C C Plot size and number of contours for the contour plot. C DATA PLTSIZ/7.5/, NCON/8/ C C Default node set: C DATA N/12/, NCC/1/, LCC(1)/9/ DATA X(1), X(2), X(3), X(4), X(5), X(6), X(7) . / 0., 1., .5, .15, .85, .5, 0./, . X(8), X(9), X(10), X(11), X(12) . / 1., .35, .65, .65, .35/, . Y(1), Y(2), Y(3), Y(4), Y(5), Y(6), Y(7) . / 0., 0., .15, .5, .5, .85, 1./, . Y(8), Y(9), Y(10), Y(11), Y(12) . / 1., .35, .35, .65, .65/ C C Evaluation points: C DATA P(1), P(2), P(3), P(4), P(5), P(6), P(7) . /-0.1, 0.0, 0.1, 0.2, 0.3, 0.4, .5/, . P(8), P(9), P(10), P(11), P(12), P(13) . / 0.6, 0.7, 0.8, 0.9, 1.0, 1.1/ C C Test function for ZGRADG, GETSIG, and CRPLOT: C FT(XF,YF) = (TANH(9.*(YF-XF)) + 1.)/9. C C Logical unit numbers for I/O: C DATA LIN/1/, LOUT/2/, LPLT/3/ OPEN (LOUT,FILE='RES') OPEN (LPLT,FILE='RES.eps') C C *** Read triangulation parameters -- N, NCC, LCC, X, Y. C C# OPEN (LIN,FILE='srftest.dat',STATUS='OLD') C# READ (LIN,100,ERR=30) N, NCC IF (N .LT. 3 .OR. N .GT. NMAX .OR. NCC .LT. 0 . .OR. NCC .GT. NCMAX) GO TO 31 C# IF (NCC .GT. 0) READ (LIN,100,ERR=30) C# . (LCC(I), I = 1,NCC) C# READ (LIN,110,ERR=30) (X(I),Y(I), I = 1,N) C#100 FORMAT (I4) C#110 FORMAT (2F13.8) C C Set NLS to the number of nonconstraint nodes. C IF (NCC .EQ. 0) THEN NLS = N ELSE NLS = LCC(1) - 1 ENDIF C C Print a heading, and construct the triangulation. C NODES and GRAD are used as work space. C WRITE (LOUT,400) N CALL TRMESH (N,X,Y, LIST,LPTR,LEND,LNEW,NODES, . NODES(N+1),GRAD,IER) IF (IER .EQ. -2) THEN WRITE (LOUT,210) STOP ELSEIF (IER .EQ. -4) THEN WRITE (LOUT,212) STOP ELSEIF (IER .GT. 0) THEN WRITE (LOUT,214) STOP ENDIF C C *** Add the constraint curves. Note that edges and C triangles are not removed from constraint regions. C ADDCST forces the inclusion of triangulation edges C connecting the sequences of constraint nodes. If it C is necessary to alter the triangulation, the empty C circumcircle property is no longer satisfied. NODES C is used as a work space array. C LW = LWK CALL ADDCST (NCC,LCC,N,X,Y, LW,NODES,LIST,LPTR, . LEND, IER) IF (IER .NE. 0) THEN WRITE (LOUT,220) IER STOP ENDIF C C *** Test VOLUME by comparing the values returned by AREAP C applied to the triangulation boundary and C VOLUME applied to the constant function C Z = 1. C CALL BNODES (N,LIST,LPTR,LEND, NODES,NB,NA,NT) A1 = AREAP(X,Y,NB,NODES) DO 1 K = 1,N Z(K) = 1. GRAD(1,K) = 0. GRAD(2,K) = 0. 1 CONTINUE A2 = VOLUME(0,LCC,N,X,Y,Z,LIST,LPTR,LEND) A3 = VOLUME(NCC,LCC,N,X,Y,Z,LIST,LPTR,LEND) WRITE (LOUT,300) NB, NA, NT, A1, A2, A3 C C *** Test GRADG on the constant function with no tension. C IFLGS = 0 SIGMA(1) = 0. NITG = 1 DGMX = 0. CALL GRADG (NCC,LCC,N,X,Y,Z,LIST,LPTR,LEND, . IFLGS,SIGMA, NITG,DGMX,GRAD, IER) IF (IER .LT. 0) THEN WRITE (LOUT,230) STOP ENDIF DMAX = 0. DO 2 K = 1,N DMAX = MAX(DMAX,ABS(GRAD(1,K)),ABS(GRAD(2,K))) 2 CONTINUE WRITE (LOUT,310) DMAX C C *** Test INTRC0 on the constant function (for which both C interpolation and extrapolation are C exact). C DMAX = 0. IST = 1 DO 4 J = 1,NJ PY = P(J) DO 3 I = 1,NI PX = P(I) CALL INTRC0 (PX,PY,NCC,LCC,N,X,Y,Z,LIST, . LPTR,LEND, IST, PZ,IER) DMAX = MAX(DMAX,ABS(PZ-1.)) 3 CONTINUE 4 CONTINUE WRITE (LOUT,320) DMAX C C *** Test GRADL on a quadratic function. The true partial C derivatives are stored in ZX and ZY. C DO 5 K = 1,N Z(K) = ((X(K)+Y(K))**2)/4. 5 CONTINUE DMAX = 0. DO 6 K = 1,N ZX = (X(K)+Y(K))/2. ZY = ZX CALL GRADL (K,NCC,LCC,N,X,Y,Z,LIST,LPTR, . LEND, GRAD(1,K),GRAD(2,K),IER) DMAX = MAX(DMAX,ABS(ZX-GRAD(1,K)), . ABS(ZY-GRAD(2,K))) 6 CONTINUE WRITE (LOUT,330) DMAX C C *** Test UNIF (and INTRC1) on the linear function (for C which both interpolation and extrapolation C are exact). C DO 7 K = 1,N Z(K) = (X(K)+Y(K))/2. GRAD(1,K) = .5 GRAD(2,K) = .5 7 CONTINUE IFLGS = 0 SIGMA(1) = 0. SFLAG = .FALSE. CALL UNIF (NCC,LCC,N,X,Y,Z,GRAD,LIST,LPTR,LEND, . IFLGS,SIGMA,NI,NI,NJ,P,P,SFLAG,SVAL, ZZ, . IER) DMAX = 0. DO 9 J = 1,NJ DO 8 I = 1,NI ZTRUE = (P(I)+P(J))/2. DMAX = MAX(DMAX,ABS(ZTRUE-ZZ(I,J))) 8 CONTINUE 9 CONTINUE WRITE (LOUT,340) DMAX, IER, NI*NJ C C *** Test GETSIG, ZINIT, and ZGRADG by computing tension C factors SIGMA, global gradient estimates GRAD, and C constraint node values (Z(I), I = LCC(1) to N) using C data values (Z(I), I = 1 to NLS) from test function C FT. Gradient estimates are initialized to zeros. C DO 10 K = 1,NLS Z(K) = FT(X(K),Y(K)) 10 CONTINUE DO 11 K = 1,N GRAD(1,K) = 0. GRAD(2,K) = 0. 11 CONTINUE C C ZINIT is used to get good initial estimates of constraint- C node values for ZGRADG. C CALL ZINIT (NCC,LCC,N,X,Y,LIST,LPTR,LEND, Z, IER) C C Initialize tension factors to zeros. L is the length of C LIST, LPTR, and SIGMA. C L = LNEW - 1 DO 12 K = 1,L SIGMA(K) = 0. 12 CONTINUE C C The global gradient estimates depend on the tension C factors, and the optimal tension factors depend on the C gradient estimates. This requires an iteration around C the calls to ZGRADG/GRADG and GETSIG. An appropriate C convergence test would involve an upper bound on the C relative change in the gradients between iterations. C However, a few iterations is usually sufficient to C produce good results. C C Parameters: C C MAXITZ = Maximum number of Gauss-Seidel iterations to be C used by ZGRADG. C DZMAX = Tolerance defining convergence in ZGRADG. C MAXITG = Maximum number of Gauss-Seidel iterations to be C used by GRADG. C DGMAX = Tolerance defining convergence in GRADG. C TOL = Tolerance for GETSIG. C IMAX = Number of ZGRADG/GETSIG iterations. C IFLGS = Tension factor option: 0 means uniform tension. C MAXITZ = 10 DZMAX = 1.E-3 MAXITG = 10 DGMAX = 1.E-3 TOL = 1.E-3 IMAX = 3 DO 13 I = 1,IMAX IFLGS = 0 IF (I .GT. 1) IFLGS = 1 C C Compute global gradient estimates and constraint node C values (ZGRADG). C IF (NCC .GT. 0 .AND. MAXITZ .GT. 0) THEN NITZ = MAXITZ DZMX = DZMAX CALL ZGRADG (NCC,LCC,N,X,Y,LIST,LPTR,LEND, . IFLGS,SIGMA, NITZ,DZMX,Z, . GRAD, IER) IF (IER .LT. 0) THEN WRITE (LOUT,240) IER, NITZ STOP ENDIF ELSE NITZ = 0 ENDIF C C Improve the gradient estimates (GRADG). C IF (MAXITG .GT. 0) THEN NITG = MAXITG DGMX = DGMAX CALL GRADG (NCC,LCC,N,X,Y,Z,LIST,LPTR,LEND, . IFLGS,SIGMA, NITG,DGMX,GRAD, IER) ELSE NITG = 0 ENDIF C C Compute tension factors SIGMA. C CALL GETSIG (N,X,Y,Z,LIST,LPTR,LEND,GRAD, . TOL, SIGMA, DSM,IER) 13 CONTINUE C C Compute interpolated values and the maximum error on the C uniform grid. C IFLGS = 1 DFLAG = .FALSE. IST = 1 DMAX = 0. DO 15 J = 1,NJ PY = P(J) DO 14 I = 1,NI PX = P(I) CALL INTRC1 (PX,PY,NCC,LCC,N,X,Y,Z,LIST,LPTR,LEND, . IFLGS,SIGMA,GRAD,DFLAG, IST, ZZ(I,J), . ZX,ZY,IER) IF (IER .LT. 0) THEN WRITE (LOUT,250) IER STOP ENDIF C C Exclude extrapolation errors. C IF (IER .GT. 0) GO TO 14 ZTRUE = FT(PX,PY) DMAX = MAX(DMAX,ABS(ZTRUE-ZZ(I,J))) 14 CONTINUE 15 CONTINUE WRITE (LOUT,350) DMAX C C Print SIGMA. C CALL SGPRNT (N,LOUT,LIST,LPTR,LEND,SIGMA) C C Create a contour plot. NODES, X, and Y are used as work C space arrays. It is assumed that NMAX >= 2.5*NI*NJ. C CALL CRPLOT (LPLT,PLTSIZ,NI,NJ,P,P,ZZ,NCON,NODES, . X,Y, IER) IF (IER .NE. 0) THEN WRITE (LOUT,260) IER STOP ENDIF WRITE (LOUT,360) C C Successful test. C STOP C C Error reading the data set. C C# 30 WRITE (*,200) STOP C C Invalid value of N or NCC. C 31 WRITE (*,205) N, NCC STOP C C Heading: C 400 FORMAT (///1X,21X,'SRFPACK Test: N =',I4///) C C Error message formats: C C#200 FORMAT (//5X,'*** Input data set invalid ***'/) 205 FORMAT (//5X,'*** N or NCC is outside its valid ', . 'range: N =',I5,', NCC = ',I4,' ***'/) 210 FORMAT (//5X,'*** Error in TRMESH: the first three ', . 'nodes are collinear ***'/) 212 FORMAT (//5X,'*** Error in TRMESH: invalid ', . 'triangulation ***'/) 214 FORMAT (//5X,'*** Error in TRMESH: duplicate nodes ', . 'encountered ***'/) 220 FORMAT (//5X,'*** Error in ADDCST: IER = ',I1, . ' ***'/) 230 FORMAT (//5X,'*** Error in GRADG (invalid flag ', . 'returned): IER = ',I2,' ***'/) 240 FORMAT (//5X,'*** Error in ZGRADG: IER = ', . I2,', NIT = ',I3,' ***'/) 250 FORMAT (//5X,'*** Error in INTRC1 (invalid flag ', . 'returned): IER = ',I2,' ***'/) 260 FORMAT (//5X,'*** Error in CRPLOT: IER = ',I1, . ' ***'/) C C Test message formats: C 300 FORMAT (//5X,'Output from BNODES, AREA, and VOLUME'// . 5X,'BNODES: ',I4,' Boundary Nodes, ',I4, . ' Arcs, ',I4,' Triangles'/5X, . 'AREAP: Area of convex hull = ',E15.8/5X, . 'VOLUME: Area of convex hull = ',E15.8/5X, . 'VOLUME: Area excluding constraint regions', . ' = ',E15.8//) 310 FORMAT (5X,'GRADG Test: Maximum error = ', . E15.8//) 320 FORMAT (5X,'INTRC0 Test: Maximum error = ', . E15.8//) 330 FORMAT (5X,'GRADL Test: Maximum error = ', . E15.8//) 340 FORMAT (5X,'UNIF (INTRC1) Test: Maximum ', . 'error = ',E15.8/5X,I2,' of ',I3, . ' points required extrapolation'//) 350 FORMAT (5X,'ZGRADG/GETSIG Test: Maximum ', . 'interpolation error = ',E15.8//) 360 FORMAT (/5X,'A contour plot file was ', . 'successfully created.') END SHAR_EOF fi # end of overwriting check cd .. cd .. cd .. if test ! -d 'Src' then mkdir 'Src' fi cd 'Src' if test ! -d 'Fortran77' then mkdir 'Fortran77' fi cd 'Fortran77' if test ! -d 'Sp' then mkdir 'Sp' fi cd 'Sp' if test -f 'srfpack.f' then echo shar: will not over-write existing file "'srfpack.f'" else cat << SHAR_EOF > 'srfpack.f' SUBROUTINE ARCINT (B,X1,X2,Y1,Y2,H1,H2,HX1,HX2,HY1, . HY2,SIGMA,DFLAG, HP,HXP,HYP,IER) INTEGER IER LOGICAL DFLAG REAL B, X1, X2, Y1, Y2, H1, H2, HX1, HX2, HY1, . HY2, SIGMA, HP, HXP, HYP C C*********************************************************** C C From SRFPACK C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C renka@cs.unt.edu C 11/18/96 C C Gigen a line segment P1-P2 containing a point P = C (XP,YP), along with function values and partial deriva- C tives at the endpoints, this subroutine computes an C interpolated value and, optionally, a gradient at P. The C value and tangential gradient component at P are taken to C be the value and derivative of the Hermite interpolatory C tension spline H defined by the endpoint values and tan- C gential gradient components. The normal gradient compo- C nent at P is obtained by linear interpolation applied to C the normal components at the endpoints. C C On input: C C B = Local coordinate of P with respect to P1-P2: C P = B*P1 + (1-B)*P2. Note that B may be comput- C ed from the coordinates of P as / C . C C X1,X2,Y1,Y2 = Coordinates of a pair of distinct C points P1 and P2. C C H1,H2 = Values of the interpolant H at P1 and P2, C respectively. C C HX1,HX2,HY1,HY2 = x and y partial derivatives of H C at P1 and P2. C C SIGMA = Tension factor associated with P1-P2. C C DFLAG = Logical flag which specifies whether first C partial derivatives at P are to be computed: C DFLAG = .TRUE. if and only if partials are C to be returned. C C Input parameters are not altered by this routine. C C On output: C C HP = Interpolated value at P unless IER < 0, in C which case HP is not defined. C C HXP,HYP = x and y partial derivatives at P unless C DFLAG = FALSE or IER < 0. C C IER = Error indicator: C IER = 0 if no errors were encountered. C IER = 1 if B < 0 or B > 1 and thus HP is an C extrapolated value. C IER = -1 if P1 and P2 coincide. C C SRFPACK module required by ARCINT: SNHCSH C C Intrinsic functions called by ARCINT: ABS, EXP C C*********************************************************** C REAL B1, B2, CM, CM2, CMM, D1, D2, DS, DUMMY, DX, DY, . E, E1, E2, EMS, GN, GT, S, S1, S2, SB1, SB2, . SBIG, SIG, SINH2, SM, SM2, TM, TM1, TM2, TP1, . TP2, TS DATA SBIG/85./ C DX = X2 - X1 DY = Y2 - Y1 DS = DX*DX + DY*DY IF (DS .EQ. 0.) GO TO 1 IER = 0 C C Compute local coordinates B1 and B2, tangential deriva- C tives S1 and S2, slope S, and second differences D1 and C D2. S1, S2, S, D1, and D2 are scaled by the separation C D between P1 and P2. C B1 = B B2 = 1. - B1 IF (B1 .LT. 0. .OR. B2 .LT. 0.) IER = 1 S1 = HX1*DX + HY1*DY S2 = HX2*DX + HY2*DY S = H2 - H1 D1 = S - S1 D2 = S2 - S C C Compute HP and, if required, the scaled tangential grad- C ient component GT. C SIG = ABS(SIGMA) IF (SIG .LT. 1.E-9) THEN C C SIG = 0: use Hermite cubic interpolation. C HP = H1 + B2*(S1 + B2*(D1 + B1*(D1 - D2))) IF (.NOT. DFLAG) RETURN GT = S1 + B2*(D1 + D2 + 3.*B1*(D1 - D2)) ELSEIF (SIG .LE. .5) THEN C C 0 .LT. SIG .LE. .5: use approximations designed to avoid C cancellation error in the hyperbolic functions. C SB2 = SIG*B2 CALL SNHCSH (SIG, SM,CM,CMM) CALL SNHCSH (SB2, SM2,CM2,DUMMY) E = SIG*SM - CMM - CMM HP = H1 + B2*S1 + ((CM*SM2-SM*CM2)*(D1+D2) + SIG* . (CM*CM2-(SM+SIG)*SM2)*D1)/(SIG*E) IF (.NOT. DFLAG) RETURN SINH2 = SM2 + SB2 GT = S1 + ((CM*CM2-SM*SINH2)*(D1+D2) + SIG* . (CM*SINH2-(SM+SIG)*CM2)*D1)/E ELSE C C SIG > .5: use negative exponentials in order to avoid C overflow. Note that EMS = EXP(-SIG). In the case of C extrapolation (negative B1 or B2), H is approximated C by a linear function if -SIG*B1 or -SIG*B2 is large. C SB1 = SIG*B1 SB2 = SIG - SB1 IF (-SB1 .GT. SBIG .OR. -SB2 .GT. SBIG) THEN HP = H1 + B2*S IF (.NOT. DFLAG) RETURN GT = S ELSE E1 = EXP(-SB1) E2 = EXP(-SB2) EMS = E1*E2 TM = 1. - EMS TS = TM*TM TM1 = 1. - E1 TM2 = 1. - E2 E = TM*(SIG*(1.+EMS) - TM - TM) HP = H1 + B2*S + (TM*TM1*TM2*(D1+D2) + SIG* . ((E2*TM1*TM1-B1*TS)*D1 + . (E1*TM2*TM2-B2*TS)*D2))/(SIG*E) IF (.NOT. DFLAG) RETURN TP1 = 1. + E1 TP2 = 1. + E2 GT = S + (TM1*(TM*TP2-SIG*E2*TP1)*D1 - . TM2*(TM*TP1-SIG*E1*TP2)*D2)/E ENDIF ENDIF C C Compute the gradient at P, (HXP,HYP) = (GT/D)T + (GN/D)N, C where T = (DX,DY)/D (unit tangent vector), N = (-DY,DX)/ C D (unit normal), and the scaled normal component is GN = C B1<(HX1,HY1),N> + B2<(HX2,HY2),N>. C GN = B1*(HY1*DX-HX1*DY) + B2*(HY2*DX-HX2*DY) HXP = (GT*DX - GN*DY)/DS HYP = (GT*DY + GN*DX)/DS RETURN C C P1 and P2 coincide. C 1 IER = -1 RETURN END SUBROUTINE CNTOUR (NX,NY,X,Y,Z,CVAL,LC,NCMAX,IWK, XC, . YC,ILC,NC,IER) INTEGER NX, NY, LC, NCMAX, IWK(NX,*), ILC(NCMAX), NC, . IER REAL X(NX), Y(NY), Z(NX,NY), CVAL, XC(LC), YC(LC) C C*********************************************************** C C From SRFPACK C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C renka@cs.unt.edu C 04/28/90 C C Given a set of function values Z = F(X,Y) at the verti- C ces of an NX by NY rectangular grid, this subroutine de- C termines a set of contour lines associated with F = CVAL. C A contour line is specified by an ordered sequence of C points (XC,YC), each lying on a grid edge and computed C from the linear interpolant of the function values at the C endpoints of the edge. The accuracy of the contour lines C is thus directly related to the number of grid points. If C a contour line forms a closed curve, the first point coin- C cides with the last point. Otherwise, the first and last C points lie on the grid boundary. C C Note that the problem is ill-conditioned in the vicinity C of a double zero of F-CVAL. Thus, if a grid cell is C crossed by two contour lines (all four sides intersected), C three different configurations are possible, corresponding C to a local minimum, a local maximum, or a saddle point. C It is arbitrarily assumed in this case that the contour C lines intersect, representing a saddle point. Also, in C order to treat the case of F = CVAL at a vertex in a con- C sistent manner, this case is always treated as F > CVAL. C Hence, if F takes on the same value at both ends of an C edge, it is assumed that no contour line intersects that C edge. In particular, a constant function, including C F = CVAL, results in no contour lines. C C On input: C C NX = Number of grid points in the x direction. C NX .GE. 2. C C NY = Number of grid points in the y direction. C NY .GE. 2. C C X = Array of length NX containing a strictly in- C creasing sequence of values. C C Y = Array of length NY containing a strictly in- C creasing sequence of values. C C Z = Array of function values at the vertices of the C rectangular grid. Z(I,J) = F(X(I),Y(J)) for C I = 1,...,NX and J = 1,...,NY. C C CVAL = Constant function value defining a contour C line as the set of points (X,Y) such that C F(X,Y) = CVAL. C C LC = Length of arrays XC and YC, and maximum allow- C able number of points defining contour lines. C LC = 2(NX-1)(NY-1) + (NX*NY+1)/2 is (probably C more than) sufficient. LC .GE. 2. C C NCMAX = Length of array ILC, and maximum allowable C number of contour lines. NCMAX = (NX*NY+1)/ C 2 is sufficient. NCMAX .GE. 1. C C The above parameters are not altered by this routine. C C IWK = Integer array of length .GE. NX*(NY-1) to be C used as work space. C C XC,YC = Arrays of length LC. C C ILC = Integer array of length NCMAX. C C On output: C C XC,YC = Arrays containing the coordinates of NC con- C tour lines. For K = 1,...,NC, contour line C K is defined by the sequence of points with C indexes ILC(K-1)+1,...,ILC(K) where ILC(0) = C 0. C C ILC = Array containing the indexes (to XC and YC) C associated with the terminal point of contour C line K in position K for K = 1,...,NC (if NC C .GT. 0). C C NC = Number of contour lines whose points are stored C in XC and YC. C C IER = Error indicator: C IER = 0 if no errors were encountered and all C contour lines were found. C IER = 1 if NX, NY, LC, or NCMAX is outside its C valid range. NC = 0 and XC, YC, and C ILC are not altered in this case. C IER = 2 if X or Y is not strictly increasing. C NC = 0 and XC, YC, and ILC are not C altered in this case. C IER = K for K > LC, where K is the required C length of XC and YC, if more storage C space is required to complete the C specification of contour line NC and/ C or additional contour lines up to a C total of NCMAX. NC .GE. 1 and ILC(NC) C = LC in this case. C IER = -1 if more than NCMAX contour lines are C present (more space is required in C ILC). NC = NCMAX, and LC may or may C not be sufficient for the additional C contour lines in this case. (This is C not determined.) C C In the unlikely event of an internal failure, a message C is printed on logical unit LUN (specified in the DATA C statement below). IER may be 0 in this case. C C Modules required by CNTOUR: None C C*********************************************************** C INTEGER I, I1, I2, IB, IN, IND, ISID, ISIDB, ISIDN, . J, J1, J2, JB, JN, K, LCON, LMX, LUN, NCMX, . NCON, NI, NIM1, NJ, NJM1 LOGICAL BDRY REAL CV, W, XF, XN, XP, YF, YN, YP, Z1, Z2 DATA LUN/0/ C C Store parameters in local variables. C NI = NX NJ = NY NIM1 = NI - 1 NJM1 = NJ - 1 CV = CVAL LMX = LC NCMX = NCMAX NC = 0 C C Test for invalid input parameters. C IER = 1 IF (NI .LT. 2 .OR. NJ .LT. 2 .OR. LMX .LT. 2 .OR. . NCMX .LT. 1) RETURN C C Test for nonincreasing values of X or Y. C IER = 2 DO 1 I = 2,NI IF (X(I) .LE. X(I-1)) RETURN 1 CONTINUE DO 2 J = 2,NJ IF (Y(J) .LE. Y(J-1)) RETURN 2 CONTINUE C C Loop on grid cells, initializing edge indicators (stored C in IWK) to zeros. For each cell, the indicator IND is a C 4-bit integer with each bit corresponding to an edge of C the cell, and having value 1 iff the edge has been pro- C cessed. Note that two IND values must be adjusted when C an interior edge is processed. The cell sides (edges) C are numbered (1,2,4,8) in counterclockwise order start- C ing from the bottom. This corresponds to an ordering of C the weighted IND bits from low order to high order. C Grid cells are identified with their lower left corners. C DO 4 J = 1,NJM1 DO 3 I = 1,NIM1 IWK(I,J) = 0 3 CONTINUE 4 CONTINUE C C First determine open contours by looping on boundary edges C in counterclockwise order starting from the lower left. C For each unprocessed boundary edge intersected by a con- C tour line, the contour line is determined and IWK is up- C dated to reflect the edges intersected. The boundary C cell (lower left corner) is indexed by (IB,JB) and the C boundary edge is specified by ISIDB. NCON and LCON are C local variables containing the number of contour lines C encountered and the current length of XC and YC. C NCON = 0 LCON = 0 ISIDB = 1 IB = 1 JB = 1 C C Top of loop on boundary edges. The edge has been C processed iff IND/ISIDB is odd. C 5 IND = IWK(IB,JB) IF (IND/ISIDB .NE. 2*((IND/ISIDB)/2)) GO TO 9 C C Update the edge indicator and store the vertex indexes of C the endpoints of the edge. C IWK(IB,JB) = IND + ISIDB IF (ISIDB .EQ. 1) THEN I1 = IB J1 = JB I2 = IB + 1 J2 = JB ELSEIF (ISIDB .EQ. 2) THEN I1 = IB + 1 J1 = JB I2 = IB + 1 J2 = JB + 1 ELSEIF (ISIDB .EQ. 4) THEN I1 = IB + 1 J1 = JB + 1 I2 = IB J2 = JB + 1 ELSE I1 = IB J1 = JB + 1 I2 = IB J2 = JB ENDIF C C Proceed to the next edge if there is no intersection. C Z1 = Z(I1,J1) Z2 = Z(I2,J2) IF ((Z1 .LT. CV .AND. Z2 .LT. CV) .OR. . (Z1 .GE. CV .AND. Z2 .GE. CV)) GO TO 9 C C Store the zero of the linear interpolant of Z1-CV and C Z2-CV as the first point of an open contour unless C NCMAX contour lines have been found or there is in- C sufficient space reserved for XC and YC. C IF (NCON .EQ. NCMX) THEN IER = -1 GO TO 16 ENDIF NCON = NCON + 1 LCON = LCON + 1 W = (CV-Z1)/(Z2-Z1) XP = X(I1) + W*(X(I2)-X(I1)) YP = Y(J1) + W*(Y(J2)-Y(J1)) IF (LCON .LE. LMX) THEN XC(LCON) = XP YC(LCON) = YP ENDIF C C Initialize for loop on cells intersected by the open C contour line. C I = IB J = JB ISID = ISIDB C C Traverse the contour line. Cell (I,J) was entered on side C ISID = (I1,J1)->(I2,J2). Find an exit edge E (unproces- C sed edge intersected by the contour) by looping on the C remaining three sides, starting with the side opposite C ISID. C 6 IND = IWK(I,J) DO 7 K = 1,3 ISID = 2*ISID IF (K .NE. 2) ISID = 2*ISID IF (ISID .GT. 15) ISID = ISID/16 IF (ISID .EQ. 1) THEN I1 = I J1 = J I2 = I + 1 J2 = J ELSEIF (ISID .EQ. 2) THEN I1 = I + 1 J1 = J I2 = I + 1 J2 = J + 1 ELSEIF (ISID .EQ. 4) THEN I1 = I + 1 J1 = J + 1 I2 = I J2 = J + 1 ELSE I1 = I J1 = J + 1 I2 = I J2 = J ENDIF C C Test for a 1 in bit position ISID of cell (I,J) and bypass C the edge if it has been previously encountered. C IF (IND/ISID .NE. 2*((IND/ISID)/2)) GO TO 7 C C Update IWK for edge E = (I1,J1)->(I2,J2). (IN,JN) indexes C the cell which shares E with cell (I,J), and ISIDN is C the side number of E in (IN,JN). BDRY is true iff E is C a boundary edge (with no neighboring cell). C IWK(I,J) = IWK(I,J) + ISID IF (ISID .LE. 2) THEN IN = I1 JN = J2 - 1 ISIDN = 4*ISID ELSE IN = I1 - 1 JN = J2 ISIDN = ISID/4 ENDIF BDRY = IN .EQ. 0 .OR. IN .EQ. NI .OR. . JN .EQ. 0 .OR. JN .EQ. NJ IF (.NOT. BDRY) IWK(IN,JN) = IWK(IN,JN) + ISIDN C C Exit the loop on sides if E is intersected by the contour. C Z1 = Z(I1,J1) Z2 = Z(I2,J2) IF ((Z1 .LT. CV .AND. Z2 .GE. CV) .OR. . (Z1 .GE. CV .AND. Z2 .LT. CV)) GO TO 8 7 CONTINUE C* C Error -- No exit point found. Print a message and exit C the contour traversal loop. C WRITE (LUN,100) NCON 100 FORMAT (///5X,'Error in CNTOUR: Contour line L ', . 'begins on the boundary'/5X,'and terminates ', . 'in the interior for L =',I4/) ILC(NCON) = LCON GO TO 9 C* C Add the intersection point (XN,YN) to the list unless it C coincides with the previous point (XP,YP) or there is C not enough space in XC and YC. C 8 W = (CV-Z1)/(Z2-Z1) XN = X(I1) + W*(X(I2)-X(I1)) YN = Y(J1) + W*(Y(J2)-Y(J1)) IF (XN .NE. XP .OR. YN .NE. YP) THEN LCON = LCON + 1 XP = XN YP = YN IF (LCON .LE. LMX) THEN XC(LCON) = XN YC(LCON) = YN ENDIF ENDIF C C Bottom of contour traversal loop. If E is not a boundary C edge, reverse the edge direction (endpoint indexes) and C update the cell index and side number. C IF (.NOT. BDRY) THEN I = I1 J = J1 I1 = I2 J1 = J2 I2 = I J2 = J I = IN J = JN ISID = ISIDN GO TO 6 ENDIF C C Update ILC with a pointer to the end of the contour line. C ILC(NCON) = LCON C C Bottom of loop on boundary edges. Update the boundary C cell index and side number, and test for termination. C 9 IF (ISIDB .EQ. 1) THEN IF (IB .LT. NIM1) THEN IB = IB + 1 ELSE ISIDB = 2 ENDIF ELSEIF (ISIDB .EQ. 2) THEN IF (JB .LT. NJM1) THEN JB = JB + 1 ELSE ISIDB = 4 ENDIF ELSEIF (ISIDB .EQ. 4) THEN IF (IB .GT. 1) THEN IB = IB - 1 ELSE ISIDB = 8 ENDIF ELSE IF (JB .GT. 1) THEN JB = JB - 1 ELSE ISIDB = 16 ENDIF ENDIF IF (ISIDB .LT. 16) GO TO 5 C C Determine closed contours by looping on interior edges -- C the first two sides (bottom and right) of each cell, C excluding boundary edges. The beginning cell is indexed C by (IB,JB), and the beginning side number is ISIDB. C DO 15 JB = 1,NJM1 DO 14 IB = 1,NIM1 DO 13 ISIDB = 1,2 IF (JB .EQ. 1 .AND. ISIDB .EQ. 1) GO TO 13 IF (IB .EQ. NIM1 .AND. ISIDB .EQ. 2) GO TO 13 C C Bypass the edge if it was previously encountered C (IND/ISIDB odd). C IND = IWK(IB,JB) IF (IND/ISIDB .NE. 2*((IND/ISIDB)/2)) GO TO 13 C C Determine the endpoint indexes of the beginning edge E = C (I1,J1)->(I2,J2), find the index (I,J) and side number C ISID of the cell which shares E with (IB,JB), and up- C date IWK. C IF (ISIDB .EQ. 1) THEN I1 = IB J1 = JB I2 = IB + 1 J2 = JB I = IB J = JB - 1 ISID = 4 ELSE I1 = IB + 1 J1 = JB I2 = IB + 1 J2 = JB + 1 I = I1 J = J1 ISID = 8 ENDIF IWK(IB,JB) = IND + ISIDB IWK(I,J) = IWK(I,J) + ISID C C Proceed to the next interior edge if there is no C intersection. C Z1 = Z(I1,J1) Z2 = Z(I2,J2) IF ((Z1 .LT. CV .AND. Z2 .LT. CV) .OR. . (Z1 .GE. CV .AND. Z2 .GE. CV)) GO TO 13 C C Store the intersection point as the first point of a C closed contour unless NCMAX contour lines have been C found or there is insufficient space in XC and YC. C IF (NCON .EQ. NCMX) THEN IER = -1 GO TO 16 ENDIF NCON = NCON + 1 LCON = LCON + 1 W = (CV-Z1)/(Z2-Z1) XP = X(I1) + W*(X(I2)-X(I1)) YP = Y(J1) + W*(Y(J2)-Y(J1)) IF (LCON .LE. LMX) THEN XC(LCON) = XP YC(LCON) = YP ENDIF XF = XP YF = YP C C Traverse the contour line. Cell (I,J) was entered on side C ISID = edge (I2,J2)->(I1,J1). Reverse the edge direc- C tion. C 10 IN = I1 JN = J1 I1 = I2 J1 = J2 I2 = IN J2 = JN IND = IWK(I,J) C C Find an exit edge E by looping on the remaining three C sides, starting with the side opposite ISID. C DO 11 K = 1,3 ISID = 2*ISID IF (K .NE. 2) ISID = 2*ISID IF (ISID .GT. 15) ISID = ISID/16 IF (ISID .EQ. 1) THEN I1 = I J1 = J I2 = I + 1 J2 = J ELSEIF (ISID .EQ. 2) THEN I1 = I + 1 J1 = J I2 = I + 1 J2 = J + 1 ELSEIF (ISID .EQ. 4) THEN I1 = I + 1 J1 = J + 1 I2 = I J2 = J + 1 ELSE I1 = I J1 = J + 1 I2 = I J2 = J ENDIF C C Bypass the edge if it has been previously encountered. C IF (IND/ISID .NE. 2*((IND/ISID)/2)) GO TO 11 C C Determine the index (IN,JN) and side number ISIDN of the C cell which shares edge E = (I1,J1)->(I2,J2) with cell C (I,J), and update IWK. C IF (ISID .LE. 2) THEN IN = I1 JN = J2 - 1 ISIDN = 4*ISID ELSE IN = I1 - 1 JN = J2 ISIDN = ISID/4 ENDIF IWK(I,J) = IWK(I,J) + ISID IWK(IN,JN) = IWK(IN,JN) + ISIDN C C Exit the loop on sides if E is intersected. C Z1 = Z(I1,J1) Z2 = Z(I2,J2) IF ((Z1 .LT. CV .AND. Z2 .GE. CV) .OR. . (Z1 .GE. CV .AND. Z2 .LT. CV)) GO TO 12 11 CONTINUE C* C Error -- No exit point found. Print a message and exit C the contour traversal loop. C WRITE (LUN,110) NCON 110 FORMAT (///5X,'Error in CNTOUR: Contour line L ', . 'is open but'/5X,'does not intersect the ', . 'boundary for L =',I4/) ILC(NCON) = LCON GO TO 13 C* C Add the intersection point to the list unless it coincides C with the previous point or there is not enough space in C XC and YC. C 12 W = (CV-Z1)/(Z2-Z1) XN = X(I1) + W*(X(I2)-X(I1)) YN = Y(J1) + W*(Y(J2)-Y(J1)) IF (XN .NE. XP .OR. YN .NE. YP) THEN LCON = LCON + 1 XP = XN YP = YN IF (LCON .LE. LMX) THEN XC(LCON) = XN YC(LCON) = YN ENDIF ENDIF C C Bottom of contour traversal loop. If the next cell is not C the beginning cell, update the cell index and side num- C ber. C IF (IN .NE. IB .OR. JN .NE. JB) THEN I = IN J = JN ISID = ISIDN GO TO 10 ENDIF C C Add the first point as the last point (unless the first C and last points already coincide), and update ILC. C IF (XP .NE. XF .OR. YP .NE. YF) THEN LCON = LCON + 1 IF (LCON .LE. LMX) THEN XC(LCON) = XF YC(LCON) = YF ENDIF ENDIF ILC(NCON) = LCON C C Bottom of loop on interior edges. C 13 CONTINUE 14 CONTINUE 15 CONTINUE IER = 0 C C Test for insufficient storage reserved for XC and YC. C 16 IF (LCON .GT. LMX) IER = LCON NC = NCON RETURN END SUBROUTINE COORDS (XP,YP,X1,X2,X3,Y1,Y2,Y3, B1,B2, . B3, IER) INTEGER IER REAL XP, YP, X1, X2, X3, Y1, Y2, Y3, B1, B2, B3 C C*********************************************************** C C From SRFPACK C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C renka@cs.unt.edu C 09/01/88 C C This subroutine computes the barycentric (areal) coordi- C nates B1, B2, and B3 of a point P with respect to the tri- C angle with vertices P1, P2, and P3: the solution to the C linear system defined by B1 + B2 + B3 = 1 and B1*P1 + C B2*P2 + B3*P3 = P. Note that B1 is a linear function C of P which satisfies B1 = 1 at P = P1 and B1 = 0 on C the triangle side P2-P3. Also, B1 < 0 if and only if C P is to the right of P2->P3 (and thus exterior to the C triangle). B2 and B3 satisfy similar properties. C C On input: C C XP,YP = Cartesian coordinates of P. C C X1,X2,X3,Y1,Y2,Y3 = Coordinates of the vertices of C the triangle P1, P2, and P3. C C Input parameters are not altered by this routine. C C On output: C C B1,B2,B3 = Barycentric coordinates unless IER = 1, C in which case the coordinates are not C defined. C C IER = Error indicator: C IER = 0 if no errors were encountered. C IER = 1 if the vertices of the triangle are C collinear. C C Modules required by COORDS: None C C*********************************************************** C REAL A, PX, PY, XP1, XP2, XP3, YP1, YP2, YP3 C PX = XP PY = YP C C Compute components of the vectors P->P1, P->P2, and P->P3. C XP1 = X1 - PX YP1 = Y1 - PY XP2 = X2 - PX YP2 = Y2 - PY XP3 = X3 - PX YP3 = Y3 - PY C C Compute subtriangle areas B1 = P->P2 X P->P3, B2 = P->P3 X C P->P1, and B3 = P->P1 X P->P2. C B1 = XP2*YP3 - XP3*YP2 B2 = XP3*YP1 - XP1*YP3 B3 = XP1*YP2 - XP2*YP1 C C Compute twice the signed area of the triangle. C A = B1 + B2 + B3 IF (A .EQ. 0.) GO TO 1 C C Normalize the coordinates. C B1 = B1/A B2 = B2/A B3 = B3/A IER = 0 RETURN C C The vertices are collinear. C 1 IER = -1 RETURN END SUBROUTINE CRPLOT (LUN,PLTSIZ,NX,NY,PX,PY,PZ,NCON,IWK, . XC,YC, IER) INTEGER LUN, NX, NY, NCON, IWK(*), IER REAL PLTSIZ, PX(NX), PY(NY), PZ(NX,NY), . XC(*), YC(*) C C*********************************************************** C C From SRFPACK C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C renka@cs.unt.edu C 04/12/97 C C Given a set of function values PZ = F(X,Y) at the ver- C tices of an NX by NY rectangular grid, this subroutine C creates a level-2 Encapsulated PostScript (EPS) file C containing a contour plot of the piecewise bilinear inter- C polant of the function values. C C The accuracy of the contour lines increases with the C number of grid points. Refer to Subroutine CNTOUR for C further details. C C C On input: C C LUN = Logical unit number in the range 0 to 99. C The unit should be opened with an appropriate C file name before the call to this routine. C C PLTSIZ = Plot size in inches. A window containing C the plot is mapped, with aspect ratio C preserved, to a rectangular viewport with C maximum side-length PLTSIZ. The viewport C is centered on the 8.5 by 11 inch page, and C its boundary is drawn. 1.0 .LE. PLTSIZ C .LE. 7.5. C C NX = Number of grid points in the x direction. C NX .GE. 2. C C NY = Number of grid points in the y direction. C NY .GE. 2. C C PX = Array of length NX containing a strictly in- C creasing sequence of values. C C PY = Array of length NY containing a strictly in- C creasing sequence of values. C C PZ = Array of function values at the vertices of the C rectangular grid. PZ(I,J) = F(PX(I),PY(J)) for C I = 1,...,NX and J = 1,...,NY. C C NCON = Number of contour values. The contour values C are uniformly distributed over the range of C PZ values. NCON .GE. 1. C C The above parameters are not altered by this routine. C C IWK = Integer array of length at least 1.5*NX*NY to C be used as work space. C C XC,YC = Real arrays of length at least 2.5*NX*NY to C be used as work space. C C On output: C C IER = Error indicator: C IER = 0 if no errors were encountered. C IER = 1 if LUN, PLTSIZ, NX, NY, or NCON is C outside its valid range. C IER = 2 if PX or PY is not strictly C increasing. C IER = 3 if the range of PZ values has zero C width (F is constant). C IER = 4 if an error was encountered in writing C to unit LUN. C IER = 5 if an unexpected error flag was re- C turned by Subroutine CNTOUR. This C should not occur. C C In the unlikely event of an internal failure, a message C is printed on the standard output device. IER may be 0 C in this case. C C Module required by CRPLOT: CNTOUR C C Intrinsic functions called by CRPLOT: CHAR, REAL C C*********************************************************** C INTEGER I, IC, IERR, IH, IPX1, IPX2, IPY1, IPY2, IW, . J, K, KV, LC, NC, NCMAX REAL CVAL, DX, DY, DZ, PZIJ, R, SFX, SFY, T, . TX, TY, ZMAX, ZMIN C C Local parameters: C C CVAL = Contour value between ZMIN and ZMAX C DX = Window width PX(NX)-PX(1)N C DY = Window height PY(NY)-PY(1) C DZ = Interval between contour values: C (ZMAX-ZMIN)/(NCON+1) C I,J = Row and column indexes for PZ C IC = Index (for IWK) of a contour line associated C with contour value CVAL: 1 to NC C IERR = Error flag for calls to CNTOUR C IH = Height of the bounding box (viewport) in C points C IPX1,IPY1 = X and y coordinates (in points) of the lower C left corner of the bounding box C IPX2,IPY2 = X and y coordinates (in points) of the upper C right corner of the bounding box C IW = Width of the bounding box in points C K = Index (for XC and YC) of a point on a contour C line C KV = DO-loop index for loop on contour values C LC = Length of arrays XC and YC C NC = Number of contour lines associated with C contour value CVAL C NCMAX = Maximum allowable value of NC C PZIJ = PZ(I,J) C R = Aspect ratio DX/DY C SFX,SFY = Scale factors for mapping window coordinates C to viewport coordinates C T = Temporary variable C TX,TY = Translation vector for mapping window coordi- C nates to viewport coordinates C ZMIN,ZMAX = Minimum and maximum of the PZ values C C C Test for error 1. C IF (LUN .LT. 0 .OR. LUN .GT. 99 .OR. . PLTSIZ .LT. 1.0 .OR. PLTSIZ .GT. 7.5 .OR. . NX .LT. 2 .OR. NY .LT. 2 .OR. NCON .LT. 1) . GO TO 11 C C Compute the aspect ratio of the window. C DX = PX(NX) - PX(1) DY = PY(NY) - PY(1) IF (DX .EQ. 0.0 .OR. DY .EQ. 0.0) GO TO 12 R = DX/DY C C Compute the range of PZ values and the interval between C contour values. C ZMIN = PZ(1,1) ZMAX = ZMIN DO 2 J = 1,NY DO 1 I = 1,NX PZIJ = PZ(I,J) IF (PZIJ .LT. ZMIN) ZMIN = PZIJ IF (PZIJ .GT. ZMAX) ZMAX = PZIJ 1 CONTINUE 2 CONTINUE DZ = (ZMAX-ZMIN)/REAL(NCON+1) IF (DZ .LE. 0.0) GO TO 13 C C Compute the lower left (IPX1,IPY1) and upper right C (IPX2,IPY2) corner coordinates of the bounding box C (the viewport). The coordinates, specified in default C user space units (points, at 72 points/inch with origin C at the lower left corner of the page), are chosen to C preserve the aspect ratio R, and to center the plot on C the 8.5 by 11 inch page. The center of the page is C (306,396), and T = PLTSIZ/2 in points. C T = 36.0*PLTSIZ IF (R .GE. 1.0) THEN IPX1 = 306 - NINT(T) IPX2 = 306 + NINT(T) IPY1 = 396 - NINT(T/R) IPY2 = 396 + NINT(T/R) ELSE IPX1 = 306 - NINT(T*R) IPX2 = 306 + NINT(T*R) IPY1 = 396 - NINT(T) IPY2 = 396 + NINT(T) ENDIF C C Output header comments. C WRITE (LUN,100,ERR=14) IPX1, IPY1, IPX2, IPY2 100 FORMAT ('%!PS-Adobe-3.0 EPSF-3.0'/ . '%%BoundingBox:',4I4/ . '%%Title: Contour Plot'/ . '%%Creator: SRFPACK'/ . '%%EndComments') C C Draw the bounding box. C WRITE (LUN,110,ERR=14) IPX1, IPY1 WRITE (LUN,120,ERR=14) IPX1, IPY2 WRITE (LUN,120,ERR=14) IPX2, IPY2 WRITE (LUN,120,ERR=14) IPX2, IPY1 WRITE (LUN,130,ERR=14) WRITE (LUN,140,ERR=14) 110 FORMAT (2I4,' moveto') 120 FORMAT (2I4,' lineto') 130 FORMAT ('closepath') 140 FORMAT ('stroke') C C Set up a mapping from the window to the viewport. C IW = IPX2 - IPX1 IH = IPY2 - IPY1 SFX = REAL(IW)/DX SFY = REAL(IH)/DY TX = IPX1 - SFX*PX(1) TY = IPY1 - SFY*PY(1) WRITE (LUN,150,ERR=14) TX, TY, SFX, SFY 150 FORMAT (2F12.6,' translate'/ . 2F12.6,' scale') C C Set the line thickness to 2 points. (Since the scale C factors are applied to everything, the width must be C specified in world coordinates.) C T = 4.0/(SFX+SFY) WRITE (LUN,160,ERR=14) T 160 FORMAT (F12.6,' setlinewidth') C C Compute parameters for CNTOUR: C C NCMAX = Maximum allowable number of contour lines C associated with each contour value. C LC = Length of arrays XC and YC and maximum allowable C number of points defining all the contour lines C associated with a contour value. C NCMAX = (NX*NY+1)/2 LC = 2*(NX-1)*(NY-1) + NCMAX C C Loop on contour values CVAL uniformly spaced in the open C interval (ZMIN,ZMAX). C CVAL = ZMIN DO 5 KV = 1,NCON CVAL = CVAL + DZ C C Compute a sequence of NC contour lines associated with C F = CVAL. For IC = 1 to NC, IWK(IC) is the index (for C XC and YC) of the last point of contour IC. C CALL CNTOUR (NX,NY,PX,PY,PZ,CVAL,LC,NCMAX, . IWK(NCMAX+1), XC,YC,IWK,NC,IERR) IF (IERR .EQ. 2) GO TO 12 IF (IERR .NE. 0) GO TO 15 C C Draw the NC contours. C IC = 0 K = 0 3 IC = IC + 1 K = K + 1 C C Create a path consisting of contour IC. C WRITE (LUN,170,ERR=14) XC(K), YC(K) 170 FORMAT (2F12.6,' moveto') 4 K = K + 1 WRITE (LUN,180,ERR=14) XC(K), YC(K) 180 FORMAT (2F12.6,' lineto') IF (K .NE. IWK(IC)) GO TO 4 C C Paint the path. C WRITE (LUN,140,ERR=14) IF (IC .NE. NC) GO TO 3 5 CONTINUE C C Output the showpage command and end-of-file indicator. C WRITE (LUN,200,ERR=14) 200 FORMAT ('showpage'/ . '%%EOF') C C HP's interpreters require a one-byte End-of-PostScript-Job C indicator (to eliminate a timeout error message): C ASCII 4. C WRITE (LUN,210,ERR=14) CHAR(4) 210 FORMAT (A1) C C No error encountered. C IER = 0 RETURN C C Invalid input parameter. C 11 IER = 1 RETURN C C PX or PY is not strictly increasing. C 12 IER = 2 RETURN C C DZ = 0. C 13 IER = 3 RETURN C C Error writing to unit LUN. C 14 IER = 4 RETURN C C Error flag returned by CNTOUR. C 15 IER = 5 RETURN END SUBROUTINE FVAL (XP,YP,X1,X2,X3,Y1,Y2,Y3,F1,F2,F3, . FX1,FX2,FX3,FY1,FY2,FY3,SIG1,SIG2, . SIG3, FP,IER) INTEGER IER REAL XP, YP, X1, X2, X3, Y1, Y2, Y3, F1, F2, . F3, FX1, FX2, FX3, FY1, FY2, FY3, SIG1, . SIG2, SIG3, FP C C*********************************************************** C C From SRFPACK C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C renka@cs.unt.edu C 03/18/90 C C Given function values and gradients at the three ver- C tices of a triangle containing a point P, this routine C computes the value of F at P where F interpolates the ver- C tex data. Along the triangle arcs, the interpolatory C function F is the Hermite interpolatory tension spline de- C fined by the values and tangential gradient components at C the endpoints, and the derivative in the direction normal C to the arc varies linearly between the normal gradient C components at the endpoints. A first-order C-1 blending C method is used to extend F to the interior of the trian- C gle. Thus, since values and gradients on an arc depend C only on the vertex data, the method results in C-1 contin- C uity when used to interpolate over a triangulation. C C The blending method consists of taking F(P) to be the C weighted sum of the values at P of the three univariate C Hermite interpolatory tension splines defined on the line C segments which join the vertices to the opposite sides and C pass through P. The tension factors for these splines are C obtained by linear interpolation between the pair of ten- C sion factors associated with the triangle sides which join C at the appropriate vertex. C C A tension factor SIGMA associated with a Hermite interp- C olatory tension spline is a nonnegative parameter which C determines the curviness of the spline. SIGMA = 0 results C in a cubic spline, and the spline approaches the linear C interpolant as SIGMA increases. C C On input: C C XP,YP = Coordinates of a point P at which an interp- C olated value is to be computed. C C X1,X2,X3,Y1,Y2,Y3 = Coordinates of the vertices of a C triangle (V1,V2,V3) containing C P. V3 is strictly to the left C of V1->V2. C C F1,F2,F3 = Values of the interpolatory function at C the vertices. C C FX1,FX2,FX3 = x components of the gradients of F at C the vertices. C C FY1,FY2,FY3 = y components of the gradients of F at C the vertices. C C SIG1,SIG2,SIG3 = Tension factors associated with the C arcs opposite V1, V2, and V3, re- C spectively. C C Input parameters are not altered by this routine. C C On output: C C FP = Interpolated value at P unless IER < 0, in C which case FP is not defined. C C IER = Error indicator: C IER = 0 if no errors were encountered. C IER = 1 if P is not contained in the triangle. C This may result from roundoff error C when P lies near an arc, and the int- C erpolated value FP is valid in that C case. C IER = -1 if the triangle vertices are C collinear. C C SRFPACK modules required by FVAL: ARCINT, COORDS, SNHCSH C C*********************************************************** C INTEGER IERR REAL B, B1, B2, B3, C1, C2, C3, DUM, FQ, FXQ, FYQ, . H1, H2, H3, PX, PY, SIG, SUM, XQ, YQ C PX = XP PY = YP C C F(P) = C1*H1(P) + C2*H2(P) + C3*H3(P) where C1, C2, and C3 C are weight functions which sum to 1, and H1, H2, and H3 C are Hermite interpolatory tension splines on the line C segments which join vertices to opposite sides and con- C tain P. C C Compute barycentric coordinates of P with respect to the C triangle. C CALL COORDS (PX,PY,X1,X2,X3,Y1,Y2,Y3, B1,B2,B3,IER) IF (IER .NE. 0) RETURN IF (B1 .LT. 0. .OR. B2 .LT. 0. .OR. B3 .LT. 0.) . IER = 1 C C Compute the coefficients of the partial interpolants. C C1 = 1 on the side opposite V1, and C1 = 0 on the other C arcs. Similarly for C2 and C3. C C1 = B2*B3 C2 = B3*B1 C3 = B1*B2 SUM = C1 + C2 + C3 IF (SUM .EQ. 0.) THEN C C P coincides with a vertex. C FP = B1*F1 + B2*F2 + B3*F3 RETURN ENDIF C C Normalize the coefficients. C C1 = C1/SUM C2 = C2/SUM C3 = C3/SUM C C For each vertex Vi, compute the intersection Q of the side C opposite Vi with the line defined by Vi and P, the value C and gradient at Q, and the partial interpolant value Hi C at P. C C Side opposite V1: C B = B2/(B2+B3) XQ = B*X2 + (1.-B)*X3 YQ = B*Y2 + (1.-B)*Y3 SIG = B*SIG3 + (1.-B)*SIG2 CALL ARCINT (B,X2,X3,Y2,Y3,F2,F3,FX2,FX3,FY2,FY3,SIG1, . .TRUE., FQ,FXQ,FYQ,IERR) CALL ARCINT (B1,X1,XQ,Y1,YQ,F1,FQ,FX1,FXQ,FY1,FYQ,SIG, . .FALSE., H1,DUM,DUM,IERR) C C Side opposite V2: C B = B3/(B3+B1) XQ = B*X3 + (1.-B)*X1 YQ = B*Y3 + (1.-B)*Y1 SIG = B*SIG1 + (1.-B)*SIG3 CALL ARCINT (B,X3,X1,Y3,Y1,F3,F1,FX3,FX1,FY3,FY1,SIG2, . .TRUE., FQ,FXQ,FYQ,IERR) CALL ARCINT (B2,X2,XQ,Y2,YQ,F2,FQ,FX2,FXQ,FY2,FYQ,SIG, . .FALSE., H2,DUM,DUM,IERR) C C Side opposite V3: C B = B1/(B1+B2) XQ = B*X1 + (1.-B)*X2 YQ = B*Y1 + (1.-B)*Y2 SIG = B*SIG2 + (1.-B)*SIG1 CALL ARCINT (B,X1,X2,Y1,Y2,F1,F2,FX1,FX2,FY1,FY2,SIG3, . .TRUE., FQ,FXQ,FYQ,IERR) CALL ARCINT (B3,X3,XQ,Y3,YQ,F3,FQ,FX3,FXQ,FY3,FYQ,SIG, . .FALSE., H3,DUM,DUM,IERR) C C Accumulate the partial interpolant values. C FP = C1*H1 + C2*H2 + C3*H3 RETURN END SUBROUTINE GETSIG (N,X,Y,H,LIST,LPTR,LEND,HXHY, . TOL, SIGMA, DSMAX,IER) INTEGER N, LIST(*), LPTR(*), LEND(N), IER REAL X(N), Y(N), H(N), HXHY(2,N), TOL, SIGMA(*), . DSMAX C C*********************************************************** C C From SRFPACK C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C renka@cs.unt.edu C 07/03/98 C C Given a triangulation of a set of nodes in the plane, C along with data values H and gradients (HX,HY) at the C nodes, this subroutine determines, for each triangulation C arc, the smallest (nonnegative) tension factor SIGMA such C that the Hermite interpolatory tension spline H(T), de- C fined by SIGMA and the endpoint values and directional C derivatives, preserves local shape properties of the data. C In order to define the shape properties on an arc, it is C convenient to map the arc to an interval (T1,T2). Then, C denoting the endpoint data values by H1,H2 and the deriva- C tives by HP1,HP2, and letting S = (H2-H1)/(T2-T1), the C data properties are C C Monotonicity: S, HP1, and HP2 are nonnegative or C nonpositive, C and C C Convexity: HP1 .LE. S .LE. HP2 or HP1 .GE. S C .GE. HP2. C C The corresponding properties of H are constant sign of the C first and second derivatives, respectively. Note that, C unless HP1 = S = HP2, infinite tension is required (and H C is linear on the interval) if S = 0 in the case of mono- C tonicity, or if HP1 = S or HP2 = S in the case of C convexity. C C Note that if gradients are to be computed by Subroutine C GRADG or function values and gradients are computed by C SMSURF, it may be desirable to alternate those computa- C tions (which require tension factors) with calls to this C subroutine. This iterative procedure should terminate C with a call to GETSIG in order to ensure that the shape C properties are preserved, and convergence can be achieved C (at the cost of optimality) by allowing only increases in C tension factors (refer to the parameter descriptions for C SIGMA, DSMAX, and IER). C C Refer to functions SIG0, SIG1, and SIG2 for means of C selecting minimum tension factors to preserve more general C properties. C C On input: C C N = Number of nodes in the triangulation. N .GE. 3. C C X,Y = Arrays of length N containing the Cartesian C coordinates of the nodes. C C H = Array of length N containing data values at the C nodes. H(I) is associated with (X(I),Y(I)). C C LIST,LPTR,LEND = Data structure defining the trian- C gulation. Refer to TRIPACK C Subroutine TRMESH. C C HXHY = Array dimensioned 2 by N whose columns con- C tain partial derivatives at the nodes (X C partials in the first row). Refer to Subrou- C tines GRADC, GRADG, GRADL, and SMSURF. C C TOL = Tolerance whose magnitude determines how close C each tension factor is to its optimal value C when nonzero finite tension is necessary and C sufficient to satisfy the constraint -- C abs(TOL) is an upper bound on the magnitude C of the smallest (nonnegative) or largest (non- C positive) value of the first or second deriva- C tive of H in the interval. Thus, the con- C straint is satisfied, but possibly with more C tension than necessary. C C The above parameters are not altered by this routine. C C SIGMA = Array of length 2*NA = 6*(N-1)-2*NB, where C NA and NB are the numbers of arcs and boun- C dary nodes, respectively, containing minimum C values of the tension factors. The tension C factors are associated with arcs in one-to- C one correspondence with LIST entries. Note C that each arc N1-N2 has two LIST entries and C thus, the tension factor is stored in both C SIGMA(I) and SIGMA(J) where LIST(I) = N2 (in C the adjacency list for N1) and LIST(J) = N1 C (in the list associated with N2). SIGMA C should be set to all zeros if minimal ten- C sion is desired, and should be unchanged C from a previous call in order to ensure con- C vergence of the iterative procedure describ- C ed in the header comments. C C On output: C C SIGMA = Array containing tension factors for which C H(T) preserves the local data properties on C each triangulation arc, with the restriction C that SIGMA(I) .LE. 85 for all I (unless the C input value is larger). The factors are as C small as possible (within the tolerance) but C not less than their input values. If infin- C ite tension is required on an arc, the cor- C responding factor is SIGMA(I) = 85 (and H C is an approximation to the linear inter- C polant on the arc), and if neither property C is satisfied by the data, then SIGMA(I) = 0 C (assuming its input value is 0), and thus H C is cubic on the arc. C C DSMAX = Maximum increase in a component of SIGMA C from its input value. C C IER = Error indicator and information flag: C IER = I if no errors were encountered and I C components of SIGMA were altered from C their input values for I .GE. 0. C IER = -1 if N < 3. SIGMA is not altered in C this case. C IER = -2 if duplicate nodes were encountered. C C TRIPACK modules required by GETSIG: LSTPTR, STORE C C SRFPACK module required by GETSIG: SNHCSH C C Intrinsic functions called by GETSIG: ABS, EXP, MAX, MIN, C SIGN, SQRT C C*********************************************************** C INTEGER LSTPTR REAL STORE INTEGER ICNT, LP1, LP2, LPL, LUN, N1, N2, NIT, NM1 REAL A, C1, C2, COSHM, COSHMM, D0, D1, D1D2, D1PD2, . D2, DMAX, DSIG, DSM, DT, DX, DY, E, EMS, EMS2, . F, F0, FMAX, FNEG, FP, FTOL, RTOL, S, S1, S2, . SBIG, SCM, SGN, SIG, SIGIN, SINHM, SSINH, SSM, . STOL, T, T0, T1, T2, TM, TP1 C DATA SBIG/85./, LUN/-1/ NM1 = N - 1 IF (NM1 .LT. 2) GO TO 11 C C Compute an absolute tolerance FTOL = abs(TOL) and a C relative tolerance RTOL = 100*Macheps. C FTOL = ABS(TOL) RTOL = 1. 1 RTOL = RTOL/2. IF (STORE(RTOL+1.) .GT. 1.) GO TO 1 RTOL = RTOL*200. C C Print a heading. C IF (LUN .GE. 0) WRITE (LUN,100) N, FTOL 100 FORMAT (///13X,'GETSIG: N =',I4,', TOL = ',E10.3//) C C Initialize change counter ICNT and maximum change DSM for C the loop on arcs. C ICNT = 0 DSM = 0. C C Loop on arcs N1-N2 for which N2 > N1. LPL points to the C last neighbor of N1. C DO 10 N1 = 1,NM1 LPL = LEND(N1) LP1 = LPL C C Top of loop on neighbors N2 of N1. C 2 LP1 = LPTR(LP1) N2 = ABS(LIST(LP1)) IF (N2 .LE. N1) GO TO 9 C C Print a message and compute parameters for the arc: DT = C arc length and SIGIN = input SIGMA value. C IF (LUN .GE. 0) WRITE (LUN,110) N1, N2 110 FORMAT (/1X,'Arc',I4,' -',I4) DX = X(N2) - X(N1) DY = Y(N2) - Y(N1) DT = SQRT(DX*DX + DY*DY) IF (DT .EQ. 0.) GO TO 12 SIGIN = SIGMA(LP1) IF (SIGIN .GE. SBIG) GO TO 9 C C Compute scaled directional derivatives S1,S2 at the end- C points (for the direction N1->N2), first difference S, C and second differences D1,D2. C S1 = HXHY(1,N1)*DX + HXHY(2,N1)*DY S2 = HXHY(1,N2)*DX + HXHY(2,N2)*DY S = H(N2) - H(N1) D1 = S - S1 D2 = S2 - S D1D2 = D1*D2 C C Test for infinite tension required to satisfy either C property. C SIG = SBIG IF ((D1D2 .EQ. 0. .AND. S1 .NE. S2) .OR. . (S .EQ. 0. .AND. S1*S2 .GT. 0.)) GO TO 8 C C Test for SIGMA = 0 sufficient. The data satisfies convex- C ity iff D1D2 .GE. 0, and D1D2 = 0 implies S1 = S = S2. C SIG = 0. IF (D1D2 .LT. 0.) GO TO 4 IF (D1D2 .EQ. 0.) GO TO 8 T = MAX(D1/D2,D2/D1) IF (T .LE. 2.) GO TO 8 TP1 = T + 1. C C Convexity: find a zero of F(SIG) = SIG*COSHM(SIG)/ C SINHM(SIG) - TP1. C C F(0) = 2-T < 0, F(TP1) .GE. 0, the derivative of F C vanishes at SIG = 0, and the second derivative of F is C .2 at SIG = 0. A quadratic approximation is used to C obtain a starting point for the Newton method. C SIG = SQRT(10.*T-20.) NIT = 0 C C Top of loop: C 3 IF (SIG .LE. .5) THEN CALL SNHCSH (SIG, SINHM,COSHM,COSHMM) T1 = COSHM/SINHM FP = T1 + SIG*(SIG/SINHM - T1*T1 + 1.) ELSE C C Scale SINHM and COSHM by 2*exp(-SIG) in order to avoid C overflow with large SIG. C EMS = EXP(-SIG) SSM = 1. - EMS*(EMS+SIG+SIG) T1 = (1.-EMS)*(1.-EMS)/SSM FP = T1 + SIG*(2.*SIG*EMS/SSM - T1*T1 + 1.) ENDIF C F = SIG*T1 - TP1 IF (LUN .GE. 0) WRITE (LUN,120) SIG, F, FP 120 FORMAT (1X,'Convexity: SIG = ',E15.8, . ', F(SIG) = ',E15.8/1X,35X,'FP(SIG) = ', . E15.8) NIT = NIT + 1 C C Test for convergence. C IF (FP .LE. 0.) GO TO 8 DSIG = -F/FP IF (ABS(DSIG) .LE. RTOL*SIG .OR. (F .GE. 0. .AND. . F .LE. FTOL) .OR. ABS(F) .LE. RTOL) GO TO 8 C C Update SIG. C SIG = SIG + DSIG GO TO 3 C C Convexity cannot be satisfied. Monotonicity can be satis- C fied iff S1*S .GE. 0 and S2*S .GE. 0 since S .NE. 0. C 4 IF (S1*S .LT. 0. .OR. S2*S .LT. 0.) GO TO 8 T0 = 3.*S - S1 - S2 D0 = T0*T0 - S1*S2 C C SIGMA = 0 is sufficient for monotonicity iff S*T0 .GE. 0 C or D0 .LE. 0. C IF (D0 .LE. 0. .OR. S*T0 .GE. 0.) GO TO 8 C C Monotonicity: find a zero of F(SIG) = sign(S)*HP(R), C where HPP(R) = 0 and HP, HPP denote derivatives of H. C F has a unique zero, F(0) < 0, and F approaches C abs(S) as SIG increases. C C Initialize parameters for the secant method. The method C uses three points: (SG0,F0), (SIG,F), and C (SNEG,FNEG), where SG0 and SNEG are defined implicitly C by DSIG = SIG - SG0 and DMAX = SIG - SNEG. C SGN = SIGN(1.,S) SIG = SBIG FMAX = SGN*(SIG*S-S1-S2)/(SIG-2.) IF (FMAX .LE. 0.) GO TO 8 STOL = RTOL*SIG F = FMAX F0 = SGN*D0/(3.*(D1-D2)) FNEG = F0 DSIG = SIG DMAX = SIG D1PD2 = D1 + D2 NIT = 0 C C Top of loop: compute the change in SIG by linear C interpolation. C 5 DSIG = -F*DSIG/(F-F0) IF (LUN .GE. 0) WRITE (LUN,130) DSIG 130 FORMAT (1X,'Monotonicity: DSIG = ',E15.8) IF ( ABS(DSIG) .GT. ABS(DMAX) .OR. . DSIG*DMAX .GT. 0. ) GO TO 7 C C Restrict the step-size such that abs(DSIG) .GE. STOL/2. C Note that DSIG and DMAX have opposite signs. C IF (ABS(DSIG) .LT. STOL/2.) DSIG = -SIGN(STOL/2., . DMAX) C C Update SIG, F0, and F. C SIG = SIG + DSIG F0 = F IF (SIG .LE. .5) THEN C C Use approximations to the hyperbolic functions designed C to avoid cancellation error with small SIG. C CALL SNHCSH (SIG, SINHM,COSHM,COSHMM) C1 = SIG*COSHM*D2 - SINHM*D1PD2 C2 = SIG*(SINHM+SIG)*D2 - COSHM*D1PD2 A = C2 - C1 E = SIG*SINHM - COSHMM - COSHMM ELSE C C Scale SINHM and COSHM by 2*exp(-SIG) in order to avoid C overflow with large SIG. C EMS = EXP(-SIG) EMS2 = EMS + EMS TM = 1. - EMS SSINH = TM*(1.+EMS) SSM = SSINH - SIG*EMS2 SCM = TM*TM C1 = SIG*SCM*D2 - SSM*D1PD2 C2 = SIG*SSINH*D2 - SCM*D1PD2 C C R is in (0,1) and well-defined iff HPP(T1)*HPP(T2) < 0. C F = FMAX IF (C1*(SIG*SCM*D1 - SSM*D1PD2) .GE. 0.) GO TO 6 A = EMS2*(SIG*TM*D2 + (TM-SIG)*D1PD2) IF (A*(C2+C1) .LT. 0.) GO TO 6 E = SIG*SSINH - SCM - SCM ENDIF C F = (SGN*(E*S2-C2) + SQRT(A*(C2+C1)))/E C C Update the number of iterations NIT. C 6 NIT = NIT + 1 IF (LUN .GE. 0) WRITE (LUN,140) NIT, SIG, F 140 FORMAT (1X,11X,I2,' -- SIG = ',E15.8,', F = ', . E15.8) C C Test for convergence. C STOL = RTOL*SIG IF (ABS(DMAX) .LE. STOL .OR. (F .GE. 0. .AND. . F .LE. FTOL) .OR. ABS(F) .LE. RTOL) GO TO 8 DMAX = DMAX + DSIG IF (F0*F .GT. 0. .AND. ABS(F) .GE. ABS(F0)) . GO TO 7 IF (F0*F .LE. 0.) THEN C C F and F0 have opposite signs. Update (SNEG,FNEG) to C (SG0,F0) so that F and FNEG always have opposite C signs. If SIG is closer to SNEG than SG0 and abs(F) C < abs(FNEG), then swap (SNEG,FNEG) with (SG0,F0). C T1 = DMAX T2 = FNEG DMAX = DSIG FNEG = F0 IF ( ABS(DSIG) .GT. ABS(T1) .AND. . ABS(F) .LT. ABS(T2) ) THEN C DSIG = T1 F0 = T2 ENDIF ENDIF GO TO 5 C C Bottom of loop: F0*F > 0 and the new estimate would C be outside of the bracketing interval of length C abs(DMAX). Reset (SG0,F0) to (SNEG,FNEG). C 7 DSIG = DMAX F0 = FNEG GO TO 5 C C Update SIGMA, ICNT, and DSM if necessary. C 8 SIG = MIN(SIG,SBIG) IF (SIG .GT. SIGIN) THEN SIGMA(LP1) = SIG LP2 = LSTPTR(LEND(N2),N1,LIST,LPTR) SIGMA(LP2) = SIG ICNT = ICNT + 1 DSM = MAX(DSM,SIG-SIGIN) ENDIF C C Bottom of loop on neighbors N2 of N1. C 9 IF (LP1 .NE. LPL) GO TO 2 10 CONTINUE C C No errors encountered. C DSMAX = DSM IER = ICNT RETURN C C N < 3 C 11 DSMAX = 0. IER = -1 RETURN C C Nodes N1 and N2 coincide. C 12 DSMAX = DSM IER = -2 RETURN END SUBROUTINE GIVENS ( A,B, C,S) REAL A, B, C, S C C*********************************************************** C C From SRFPACK C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C renka@cs.unt.edu C 09/01/88 C C This subroutine constructs the Givens plane rotation, C C ( C S) C G = ( ) , where C*C + S*S = 1, C (-S C) C C which zeros the second component of the vector (A,B)**T C (transposed). Subroutine ROTATE may be called to apply C the transformation to a 2 by N matrix. C C This routine is identical to Subroutine SROTG from the C LINPACK BLAS (Basic Linear Algebra Subroutines). C C On input: C C A,B = Components of the vector defining the rota- C tion. These are overwritten by values R C and Z (described below) which define C and S. C C On output: C C A = Signed Euclidean norm R of the input vector: C R = +/-SQRT(A*A + B*B) C C B = Value Z such that: C C = SQRT(1-Z*Z) and S=Z if ABS(Z) .LE. 1, and C C = 1/Z and S = SQRT(1-C*C) if ABS(Z) > 1. C C C = +/-(A/R) or 1 if R = 0. C C S = +/-(B/R) or 0 if R = 0. C C Modules required by GIVENS: None C C Intrinsic functions called by GIVENS: ABS, SQRT C C*********************************************************** C REAL AA, BB, R, U, V C C Local parameters: C C AA,BB = Local copies of A and B C R = C*A + S*B = +/-SQRT(A*A+B*B) C U,V = Variables used to scale A and B for computing R C AA = A BB = B IF (ABS(AA) .LE. ABS(BB)) GO TO 1 C C ABS(A) > ABS(B). C U = AA + AA V = BB/U R = SQRT(.25 + V*V) * U C = AA/R S = V * (C + C) C C Note that R has the sign of A, C > 0, and S has C SIGN(A)*SIGN(B). C B = S A = R RETURN C C ABS(A) .LE. ABS(B). C 1 IF (BB .EQ. 0.) GO TO 2 U = BB + BB V = AA/U C C Store R in A. C A = SQRT(.25 + V*V) * U S = BB/A C = V * (S + S) C C Note that R has the sign of B, S > 0, and C has C SIGN(A)*SIGN(B). C B = 1. IF (C .NE. 0.) B = 1./C RETURN C C A = B = 0. C 2 C = 1. S = 0. RETURN END SUBROUTINE GRADC (K,NCC,LCC,N,X,Y,Z,LIST,LPTR, . LEND, DX,DY,DXX,DXY,DYY,IER) INTEGER K, NCC, LCC(*), N, LIST(*), LPTR(*), . LEND(N), IER REAL X(N), Y(N), Z(N), DX, DY, DXX, DXY, DYY C C*********************************************************** C C From SRFPACK C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C renka@cs.unt.edu C 02/22/97 C C Given a Delaunay triangulation of N points in the plane C with associated data values Z, this subroutine estimates C first and second partial derivatives at node K. The der- C ivatives are taken to be the partials at K of a cubic C function which interpolates Z(K) and fits the data values C at a set of nearby nodes in a weighted least squares C sense. A Marquardt stabilization factor is used if neces- C sary to ensure a well-conditioned system. Thus, a unique C solution exists if there are at least 10 noncollinear C nodes. C C The triangulation may include constraints introduced by C Subroutine ADDCST, in which case the derivative estimates C are influenced by the nonconvex geometry of the domain. C Refer to Subroutine GETNP. If data values at the con- C straint nodes are not known, Subroutine ZGRADL, which C computes approximate data values at constraint nodes along C with gradients, should be called in place of this routine. C C Subroutine GRADL uses a quadratic polynomial instead of C the cubic and may be more accurate if the nodal distribu- C tion is sparse. Another alternative routine, GRADG, C employs a global method to compute the first partial C derivatives at all of the nodes at once. That method C is usually more efficient (when all first partials are C needed) and may be more accurate, depending on the data. C C On input: C C K = Index of the node at which derivatives are to be C estimated. 1 .LE. K .LE. N. C C NCC = Number of constraint curves (refer to TRIPACK C Subroutine ADDCST). NCC .GE. 0. C C LCC = Array of length NCC (or dummy array of length C 1 if NCC = 0) containing the index of the C first node of constraint I in LCC(I). For I = C 1 to NCC, LCC(I+1)-LCC(I) .GE. 3, where C LCC(NCC+1) = N+1. C C N = Number of nodes in the triangulation. C N .GE. 10. C C X,Y = Arrays of length N containing the coordinates C of the nodes with non-constraint nodes in the C first LCC(1)-1 locations, followed by NCC se- C quences of constraint nodes. C C Z = Array of length N containing data values associ- C ated with the nodes. C C LIST,LPTR,LEND = Data structure defining the trian- C gulation. Refer to TRIPACK C Subroutine TRMESH. C C Input parameters are not altered by this routine. C C On output: C C DX,DY = Estimated first partial derivatives at node C K unless IER < 0. C C DXX,DXY,DYY = Estimated second partial derivatives C at node K unless IER < 0. C C IER = Error indicator: C IER = L > 0 if no errors were encountered and C L nodes (including node K) were C employed in the least squares fit. C IER = -1 if K, NCC, an LCC entry, or N is C outside its valid range on input. C IER = -2 if all nodes are collinear. C C TRIPACK modules required by GRADC: GETNP, INTSEC C C SRFPACK modules required by GRADC: GIVENS, ROTATE, SETRO3 C C Intrinsic functions called by GRADC: ABS, MIN, REAL, SQRT C C*********************************************************** C INTEGER LMN, LMX PARAMETER (LMN=14, LMX=30) INTEGER I, IERR, J, JP1, KK, L, LMAX, LMIN, LM1, . LNP, NP, NPTS(LMX) REAL A(10,10), C, DIST(LMX), DMIN, DS, DTOL, RIN, . RS, RTOL, S, SF, SFC, SFS, STF, SUM, W, XK, . YK, ZK DATA RTOL/1.E-5/, DTOL/.01/ C C Local parameters: C C A = Transpose of the augmented regression matrix C C = First component of the plane rotation deter- C mined by Subroutine GIVENS C DIST = Array containing the distances between K and C the elements of NPTS (refer to GETNP) C DMIN = Minimum of the magnitudes of the diagonal C elements of the regression matrix after C zeros are introduced below the diagonal C DS = Squared distance between nodes K and NPTS(LNP) C DTOL = Tolerance for detecting an ill-conditioned C system. The system is accepted when DMIN/W C .GE. DTOL. C I = DO-loop index C IERR = Error flag for calls to GETNP C J = DO-loop index C JP1 = J+1 C KK = Local copy of K C L = Number of columns of A**T to which a rotation C is applied C LMAX,LMIN = Min(LMX,N), Min(LMN,N) C LMN,LMX = Minimum and maximum values of LNP for N C sufficiently large. In most cases LMN-1 C nodes are used in the fit. 4 .LE. LMN .LE. C LMX. C LM1 = LMIN-1 or LNP-1 C LNP = Length of NPTS C NP = Element of NPTS to be added to the system C NPTS = Array containing the indexes of a sequence of C nodes ordered by distance from K. NPTS(1)=K C and the first LNP-1 elements of NPTS are C used in the least squares fit. Unless LNP C exceeds LMAX, NPTS(LNP) determines R. C RIN = Inverse of the distance R between node K and C NPTS(LNP) or some point further from K than C NPTS(LMAX) if NPTS(LMAX) is used in the fit. C R is a radius of influence which enters into C the weight W. C RS = R*R C RTOL = Tolerance for determining R. If the relative C change in DS between two elements of NPTS is C not greater than RTOL, they are treated as C being the same distance from node K. C S = Second component of the plane rotation deter- C mined by Subroutine GIVENS C SF = Scale factor for the linear terms (columns 8 C and 9) in the least squares fit -- inverse C of the root-mean-square distance between K C and the nodes (other than K) in the least C squares fit C SFS = Scale factor for the quadratic terms (columns C 5, 6, and 7) in the least squares fit -- C SF*SF C SFC = Scale factor for the cubic terms (first 4 C columns) in the least squares fit -- SF**3 C STF = Marquardt stabilization factor used to damp C out the first 4 solution components (third C partials of the cubic) when the system is C ill-conditioned. As STF increases, the C fitting function approaches a quadratic C polynomial. C SUM = Sum of squared distances between node K and C the nodes used in the least squares fit C W = Weight associated with a row of the augmented C regression matrix -- 1/D - 1/R, where D < R C and D is the distance between K and a node C entering into the least squares fit C XK,YK,ZK = Coordinates and data value associated with K C KK = K C C Test for errors and initialize LMIN and LMAX. C IF (KK .LT. 1 .OR. KK .GT. N .OR. NCC .LT. 0 . .OR. N .LT. 10) GO TO 13 LMIN = MIN(LMN,N) LMAX = MIN(LMX,N) C C Compute NPTS, DIST, LNP, SF, SFS, SFC, and RIN -- C C Set NPTS to the closest LMIN-1 nodes to K. C SUM = 0. NPTS(1) = KK DIST(1) = 0. LM1 = LMIN - 1 DO 1 LNP = 2,LM1 CALL GETNP (NCC,LCC,N,X,Y,LIST,LPTR,LEND, . LNP, NPTS,DIST, IERR) IF (IERR .NE. 0) GO TO 13 DS = DIST(LNP)**2 SUM = SUM + DS 1 CONTINUE C C Add additional nodes to NPTS until the relative increase C in DS is at least RTOL. C DO 3 LNP = LMIN,LMAX CALL GETNP (NCC,LCC,N,X,Y,LIST,LPTR,LEND, . LNP, NPTS,DIST, IERR) RS = DIST(LNP)**2 IF ((RS-DS)/DS .LE. RTOL) GO TO 2 IF (LNP .GT. 10) GO TO 4 2 SUM = SUM + RS 3 CONTINUE C C Use all LMAX nodes in the least squares fit. RS is C arbitrarily increased by 10 per cent. C RS = 1.1*RS LNP = LMAX + 1 C C There are LNP-2 equations corresponding to nodes NPTS(2), C ...,NPTS(LNP-1). C 4 SFS = REAL(LNP-2)/SUM SF = SQRT(SFS) SFC = SF*SFS RIN = 1./SQRT(RS) XK = X(KK) YK = Y(KK) ZK = Z(KK) C C A Q-R decomposition is used to solve the least squares C system. The transpose of the augmented regression C matrix is stored in A with columns (rows of A) defined C as follows: 1-4 are the cubic terms, 5-7 are the quad- C ratic terms with coefficients DXX/2, DXY, and DYY/2, C 8 and 9 are the linear terms with coefficients DX and C DY, and the last column is the right hand side. C C Set up the first 9 equations and zero out the lower tri- C angle with Givens rotations. C DO 6 I = 1,9 NP = NPTS(I+1) W = 1./DIST(I+1) - RIN CALL SETRO3 (XK,YK,ZK,X(NP),Y(NP),Z(NP),SF,SFS, . SFC,W, A(1,I)) IF (I .EQ. 1) GO TO 6 DO 5 J = 1,I-1 JP1 = J + 1 L = 10 - J CALL GIVENS (A(J,J),A(J,I),C,S) CALL ROTATE (L,C,S,A(JP1,J),A(JP1,I)) 5 CONTINUE 6 CONTINUE C C Add the additional equations to the system using C the last column of A. I .LE. LNP. C I = 11 7 IF (I .LT. LNP) THEN NP = NPTS(I) W = 1./DIST(I) - RIN CALL SETRO3 (XK,YK,ZK,X(NP),Y(NP),Z(NP),SF,SFS, . SFC,W, A(1,10)) DO 8 J = 1,9 JP1 = J + 1 L = 10 - J CALL GIVENS (A(J,J),A(J,10),C,S) CALL ROTATE (L,C,S,A(JP1,J),A(JP1,10)) 8 CONTINUE I = I + 1 GO TO 7 ENDIF C C Test the system for ill-conditioning. C DMIN = MIN( ABS(A(1,1)),ABS(A(2,2)),ABS(A(3,3)), . ABS(A(4,4)),ABS(A(5,5)),ABS(A(6,6)), . ABS(A(7,7)),ABS(A(8,8)),ABS(A(9,9)) ) IF (DMIN/W .GE. DTOL) GO TO 12 IF (LNP .LE. LMAX) THEN C C Add another node to the system and increase R. Note C that I = LNP. C LNP = LNP + 1 IF (LNP .LE. LMAX) THEN CALL GETNP (NCC,LCC,N,X,Y,LIST,LPTR,LEND, . LNP, NPTS,DIST, IERR) RS = DIST(LNP)**2 ENDIF RIN = 1./SQRT(1.1*RS) GO TO 7 ENDIF C C Stabilize the system by damping third partials -- add C multiples of the first four unit vectors to the first C four equations. C STF = W DO 11 I = 1,4 A(I,10) = STF DO 9 J = I+1,10 A(J,10) = 0. 9 CONTINUE DO 10 J = I,9 JP1 = J + 1 L = 10 - J CALL GIVENS (A(J,J),A(J,10),C,S) CALL ROTATE (L,C,S,A(JP1,J),A(JP1,10)) 10 CONTINUE 11 CONTINUE C C Test the damped system for ill-conditioning. C DMIN = MIN( ABS(A(5,5)),ABS(A(6,6)),ABS(A(7,7)), . ABS(A(8,8)),ABS(A(9,9)) ) IF (DMIN/W .LT. DTOL) GO TO 14 C C Solve the 9 by 9 triangular system for the last 5 C components (first and second partial derivatives). C 12 DY = A(10,9)/A(9,9) DX = (A(10,8)-A(9,8)*DY)/A(8,8) DYY = (A(10,7)-A(8,7)*DX-A(9,7)*DY)/A(7,7) DXY = (A(10,6)-A(7,6)*DYY-A(8,6)*DX-A(9,6)*DY)/A(6,6) DXX = (A(10,5)-A(6,5)*DXY-A(7,5)*DYY-A(8,5)*DX- . A(9,5)*DY)/A(5,5) C C Scale the solution components. C DX = SF*DX DY = SF*DY DXX = 2.*SFS*DXX DXY = SFS*DXY DYY = 2.*SFS*DYY IER = LNP - 1 RETURN C C Invalid input parameter. C 13 IER = -1 RETURN C C No unique solution due to collinear nodes. C 14 IER = -2 RETURN END SUBROUTINE GRADG (NCC,LCC,N,X,Y,Z,LIST,LPTR,LEND, . IFLGS,SIGMA, NIT,DGMAX,GRAD, IER) INTEGER NCC, LCC(*), N, LIST(*), LPTR(*), LEND(N), . IFLGS, NIT, IER REAL X(N), Y(N), Z(N), SIGMA(*), DGMAX, GRAD(2,N) C C*********************************************************** C C From SRFPACK C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C renka@cs.unt.edu C 11/12/94 C C Given a triangulation of N nodes in the plane, along C with data values at the nodes and tension factors associ- C ated with the arcs, this subroutine employs a global C method to compute estimated gradients at the nodes. The C method consists of minimizing a quadratic functional Q(G) C over vectors G of length 2N (gradient components), where C Q is an approximation to the linearized curvature over the C triangulation of a C-1 bivariate function F(X,Y) which C interpolates the nodal values and gradients. C C The restriction of F to an arc of the triangulation is C taken to be the Hermite interpolatory tension spline C defined by the data values and tangential gradient compo- C nents at the endpoints of the arc, and Q is the sum over C the triangulation arcs, excluding interior constraint C arcs, of the linearized curvatures of F along the arcs -- C the integrals over the arcs of D2F(T)**2, where D2F(T) is C the second derivative of F with respect to distance T C along the arc. C C Subroutines INTRC1 and UNIF may be called to evaluate F C at arbitrary points. The interpolant F is further de- C scribed in Subroutines FVAL and TVAL, and Q is identical C to the functional Q1 described in Subroutine SMSURF. C C The minimization problem corresponds to an order 2N C symmetric positive definite sparse linear system which is C solved for the X and Y partial derivatives by the block C Gauss-Seidel method with N blocks of order 2. C C If constraints are present and data values at the con- C straint nodes are not known, Subroutine ZGRADG, which C computes approximate data values at constraint nodes C along with the gradients, should be called in place of C this routine. C C An alternative method, Subroutine GRADC or GRADL, com- C putes a local approximation to the partials at a single C node and may be more accurate, depending on the data C values and distribution of nodes (neither method emerged C as superior in tests for accuracy). If all gradients are C required and a uniform tension factor SIGMA = 0 is used, C GRADG is significantly faster than either GRADC or GRADL. C C On input: C C NCC = Number of constraint curves (refer to TRIPACK C Subroutine ADDCST). NCC .GE. 0. C C LCC = Array of length NCC (or dummy array of length C 1 if NCC = 0) containing the index of the C first node of constraint I in LCC(I). For I = C 1 to NCC, LCC(I+1)-LCC(I) .GE. 3, where C LCC(NCC+1) = N+1. C C N = Number of nodes in the triangulation. N .GE. 3. C C X,Y = Arrays of length N containing the coordinates C of the nodes with non-constraint nodes in the C first LCC(1)-1 locations, followed by NCC se- C quences of constraint nodes. C C Z = Array of length N containing data values at the C nodes. Z(I) is associated with (X(I),Y(I)). C C LIST,LPTR,LEND = Data structure defining the trian- C gulation. Refer to TRIPACK C Subroutine TRMESH. C C IFLGS = Tension factor option: C IFLGS .LE. 0 if a single uniform tension C factor is to be used. C IFLGS .GE. 1 if variable tension is desired. C C SIGMA = Uniform tension factor (IFLGS .LE. 0), or C array containing tension factors associated C with arcs in one-to-one correspondence with C LIST entries (IFLGS .GE. 1). Refer to Sub- C routines GETSIG, SIG0, SIG1, and SIG2. C C The above parameters are not altered by this routine. C C NIT = Maximum number of Gauss-Seidel iterations to C be employed. This maximum will likely be C achieved if DGMAX is smaller than the machine C precision. Note that complete convergence is C not necessary to achieve maximum accuracy of C the interpolant. For SIGMA = 0, optimal ef- C ficiency was achieved in testing with DGMAX = C 0, and NIT = 3 or 4. NIT > 0. C C DGMAX = Nonnegative convergence criterion. The C method is terminated when the maximum change C in a gradient between iterations is at most C DGMAX. The change in a gradient is taken to C be the Euclidean norm of the difference rel- C ative to 1 plus the norm of the old value. C DGMAX = 1.E-3 is sufficient for effective C convergence. C C GRAD = 2 by N array whose columns contain initial C estimates of the partial derivatives. Zero C vectors are sufficient. C C On output: C C NIT = Number of Gauss-Seidel iterations employed. C C DGMAX = Maximum relative change in a gradient at the C last iteration. C C GRAD = Estimated X and Y partial derivatives at the C nodes with X partials in the first row. Grad C is not altered if IER = -1. C C IER = Error indicator: C IER = 0 if no errors were encountered and the C convergence criterion was achieved. C IER = 1 if no errors were encountered but con- C vergence was not achieved within NIT C iterations. C IER = -1 if NCC, an LCC entry, N, NIT, or C DGMAX is outside its valid range on C input. C IER = -2 if all nodes are collinear or the C triangulation data structure is in- C valid. C IER = -3 if duplicate nodes were encountered. C C SRFPACK modules required by GRADG: GRCOEF, SNHCSH C C Intrinsic functions called by GRADG: ABS, MAX, SQRT C C*********************************************************** C INTEGER I, IFL, IFRST, ILAST, ITER, J, K, KBAK, KFOR, . LCC1, LP, LPJ, LPL, MAXIT, NB, NN REAL A11, A12, A22, D, DCUB, DELX, DELXS, DELY, . DELYS, DET, DF, DGMX, DSQ, DZX, DZY, R1, R2, . SDF, SIG, T, TOL, XK, YK, ZK, ZXK, ZYK C NN = N IFL = IFLGS MAXIT = NIT TOL = DGMAX C C Test for errors in input parameters. C IF (NCC .LT. 0 .OR. MAXIT .LT. 1 .OR. TOL .LT. 0.) . GO TO 9 LCC1 = NN+1 IF (NCC .EQ. 0) THEN IF (NN .LT. 3) GO TO 9 ELSE DO 1 I = NCC,1,-1 IF (LCC1-LCC(I) .LT. 3) GO TO 9 LCC1 = LCC(I) 1 CONTINUE IF (LCC1 .LT. 1) GO TO 9 ENDIF C C Initialize iteration count and SIG (overwritten if C IFLGS > 0). C ITER = 0 SIG = SIGMA(1) C C Top of iteration loop: If K is a constraint node, I C indexes the constraint containing node K, IFRST and C ILAST are the first and last nodes of constraint I, C and (KBAK,K,KFOR) is a subsequence of constraint I. C 2 IF (ITER .EQ. MAXIT) GO TO 8 DGMX = 0. I = 0 IFRST = 1 ILAST = LCC1-1 KBAK = 0 KFOR = 0 C C Loop on nodes. C DO 7 K = 1,NN IF (K .GE. LCC1) THEN IF (K .GT. ILAST) THEN I = I + 1 IFRST = K IF (I .LT. NCC) THEN ILAST = LCC(I+1) - 1 ELSE ILAST = NN ENDIF KBAK = ILAST KFOR = K + 1 ELSE KBAK = K - 1 IF (K .LT. ILAST) THEN KFOR = K + 1 ELSE KFOR = IFRST ENDIF ENDIF ENDIF XK = X(K) YK = Y(K) ZK = Z(K) ZXK = GRAD(1,K) ZYK = GRAD(2,K) C C Initialize components of the order 2 system for the C change (DZX,DZY) in the K-th solution components C (symmetric matrix in A and residual in R). C A11 = 0. A12 = 0. A22 = 0. R1 = 0. R2 = 0. C C Loop on neighbors J of node K. C LPL = LEND(K) LPJ = LPL 3 LPJ = LPTR(LPJ) J = ABS(LIST(LPJ)) C C Arc K-J lies in a constraint region and is bypassed iff C K and J are nodes in the same constraint and J follows C KFOR and precedes KBAK as a neighbor of K. C IF (K .LT. LCC1 .OR. J .LT. IFRST .OR. . J .GT. ILAST) GO TO 5 IF (J .EQ. KBAK .OR. J .EQ. KFOR) GO TO 5 LP = LPJ C 4 LP = LPTR(LP) NB = ABS(LIST(LP)) IF (NB .EQ. KBAK) GO TO 6 IF (NB .NE. KFOR) GO TO 4 C C Compute parameters associated with edge C K->J, and test for duplicate nodes. C 5 DELX = X(J) - XK DELY = Y(J) - YK DELXS = DELX*DELX DELYS = DELY*DELY DSQ = DELXS + DELYS D = SQRT(DSQ) DCUB = D*DSQ IF (D .EQ. 0.) GO TO 11 IF (IFL .GE. 1) SIG = SIGMA(LPJ) CALL GRCOEF (SIG,DCUB, DF,SDF) C C Update the system components for node J. The contribu- C tion from edge K->J is weighted by 1/D, where D is C the arc length. C A11 = A11 + DF*DELXS/D A12 = A12 + DF*DELX*DELY/D A22 = A22 + DF*DELYS/D T = ((DF+SDF)*(Z(J)-ZK) - DF*(ZXK*DELX + ZYK*DELY) . - SDF*(GRAD(1,J)*DELX + GRAD(2,J)*DELY))/D R1 = R1 + T*DELX R2 = R2 + T*DELY C C Bottom of loop on neighbors. C 6 IF (LPJ .NE. LPL) GO TO 3 C C Solve the system associated with the K-th block. C DET = A11*A22 - A12*A12 IF (DET .EQ. 0. .OR. A11 .EQ. 0.) GO TO 10 DZY = (A11*R2 - A12*R1)/DET DZX = (R1 - A12*DZY)/A11 C C Update the partials at node K and the maximum relative C change DGMX. C GRAD(1,K) = ZXK + DZX GRAD(2,K) = ZYK + DZY DGMX = MAX(DGMX,SQRT(DZX*DZX+DZY*DZY)/ . (1.+SQRT(ZXK*ZXK+ZYK*ZYK))) 7 CONTINUE C C Increment ITER and test for convergence. C ITER = ITER + 1 IF (DGMX .GT. TOL) GO TO 2 C C Method converged. C NIT = ITER DGMAX = DGMX IER = 0 RETURN C C Method failed to converge within NIT iterations. C 8 DGMAX = DGMX IER = 1 RETURN C C Invalid input parameter. C 9 NIT = 0 DGMAX = 0. IER = -1 RETURN C C Node K and its neighbors are collinear, resulting in a C singular system. C 10 NIT = 0 DGMAX = DGMX IER = -2 RETURN C C Nodes K and J coincide. C 11 NIT = 0 DGMAX = DGMX IER = -3 RETURN END SUBROUTINE GRADL (K,NCC,LCC,N,X,Y,Z,LIST,LPTR, . LEND, DX,DY,IER) INTEGER K, NCC, LCC(*), N, LIST(*), LPTR(*), . LEND(N), IER REAL X(N), Y(N), Z(N), DX, DY C C*********************************************************** C C From SRFPACK C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C renka@cs.unt.edu C 02/22/97 C C Given a Delaunay triangulation of N points in the plane C with associated data values Z, this subroutine estimates C X and Y partial derivatives at node K. The derivatives C are taken to be the partials at K of a quadratic function C which interpolates Z(K) and fits the data values at a set C of nearby nodes in a weighted least squares sense. A Mar- C quardt stabilization factor is used if necessary to ensure C a well-conditioned system. Thus, a unique solution exists C if there are at least 6 noncollinear nodes. C C The triangulation may include constraints introduced by C Subroutine ADDCST, in which case the gradient estimates C are influenced by the nonconvex geometry of the domain. C Refer to Subroutine GETNP. If data values at the con- C straint nodes are not known, Subroutine ZGRADL, which C computes approximate data values at constraint nodes along C with gradients, should be called in place of this routine. C C Subroutine GRADC uses a cubic polynomial instead of the C quadratic and is generally more accurate than this routine C if the nodal distribution is sufficiently dense. Another C alternative routine, GRADG, employs a global method to C compute the partial derivatives at all of the nodes at C once. That method is usually more efficient (when all C partials are needed) and may be more accurate, depending C on the data. C C On input: C C K = Index of the node at which derivatives are to be C estimated. 1 .LE. K .LE. N. C C NCC = Number of constraint curves (refer to TRIPACK C Subroutine ADDCST). NCC .GE. 0. C C LCC = Array of length NCC (or dummy array of length C 1 if NCC = 0) containing the index of the C first node of constraint I in LCC(I). For I = C 1 to NCC, LCC(I+1)-LCC(I) .GE. 3, where C LCC(NCC+1) = N+1. C C N = Number of nodes in the triangulation. N .GE. 6. C C X,Y = Arrays of length N containing the coordinates C of the nodes with non-constraint nodes in the C first LCC(1)-1 locations, followed by NCC se- C quences of constraint nodes. C C Z = Array of length N containing data values associ- C ated with the nodes. C C LIST,LPTR,LEND = Data structure defining the trian- C gulation. Refer to TRIPACK C Subroutine TRMESH. C C Input parameters are not altered by this routine. C C On output: C C DX,DY = Estimated partial derivatives at node K C unless IER < 0. C C IER = Error indicator: C IER = L > 0 if no errors were encountered and C L nodes (including node K) were C employed in the least squares fit. C IER = -1 if K, NCC, an LCC entry, or N is C outside its valid range on input. C IER = -2 if all nodes are collinear. C C TRIPACK modules required by GRADL: GETNP, INTSEC C C SRFPACK modules required by GRADL: GIVENS, ROTATE, SETRO1 C C Intrinsic functions called by GRADL: ABS, MIN, REAL, SQRT C C*********************************************************** C INTEGER LMN, LMX PARAMETER (LMN=10, LMX=30) INTEGER I, IERR, J, JP1, KK, L, LMAX, LMIN, LM1, . LNP, NP, NPTS(LMX) REAL A(6,6), C, DIST(LMX), DMIN, DS, DTOL, RIN, . RS, RTOL, S, SF, SFS, STF, SUM, W, XK, YK, . ZK DATA RTOL/1.E-5/, DTOL/.01/ C C Local parameters: C C A = Transpose of the augmented regression matrix C C = First component of the plane rotation deter- C mined by Subroutine GIVENS C DIST = Array containing the distances between K and C the elements of NPTS (refer to GETNP) C DMIN = Minimum of the magnitudes of the diagonal C elements of the regression matrix after C zeros are introduced below the diagonal C DS = Squared distance between nodes K and NPTS(LNP) C DTOL = Tolerance for detecting an ill-conditioned C system. The system is accepted when DMIN/W C .GE. DTOL C I = DO-loop index C IERR = Error flag for calls to GETNP C J = DO-loop index C JP1 = J+1 C KK = Local copy of K C L = Number of columns of A**T to which a rotation C is applied C LMAX,LMIN = Min(LMX,N), Min(LMN,N) C LMN,LMX = Minimum and maximum values of LNP for N C sufficiently large. In most cases LMN-1 C nodes are used in the fit. 4 .LE. LMN .LE. C LMX. C LM1 = LMIN-1 or LNP-1 C LNP = Length of NPTS C NP = Element of NPTS to be added to the system C NPTS = Array containing the indexes of a sequence of C nodes ordered by distance from K. NPTS(1)=K C and the first LNP-1 elements of NPTS are C used in the least squares fit. Unless LNP C exceeds LMAX, NPTS(LNP) determines R. C RIN = Inverse of the distance R between node K and C NPTS(LNP) or some point further from K than C NPTS(LMAX) if NPTS(LMAX) is used in the fit. C R is a radius of influence which enters into C the weight W. C RS = R*R C RTOL = Tolerance for determining R. If the relative C change in DS between two elements of NPTS is C not greater than RTOL, they are treated as C being the same distance from node K C S = Second component of the plane rotation deter- C mined by Subroutine GIVENS C SF = Scale factor for the linear terms (columns 4 C and 5) in the least squares fit -- inverse C of the root-mean-square distance between K C and the nodes (other than K) in the least C squares fit. C SFS = Scale factor for the quadratic terms (first 3 C columns) in the least squares fit -- SF*SF. C STF = Marquardt stabilization factor used to damp C out the first 3 solution components (second C partials of the quadratic) when the system C is ill-conditioned. As STF increases, the C fitting function approaches a linear C SUM = Sum of squared distances between node K and C the nodes used in the least squares fit C W = Weight associated with a row of the augmented C regression matrix -- 1/R - 1/D, where D < R C and D is the distance between K and a node C entering into the least squares fit. C XK,YK,ZK = Coordinates and data value associated with K C KK = K C C Test for errors and initialize LMIN and LMAX. C IF (KK .LT. 1 .OR. KK .GT. N .OR. NCC .LT. 0 . .OR. N .LT. 6) GO TO 13 LMIN = MIN(LMN,N) LMAX = MIN(LMX,N) C C Compute NPTS, DIST, LNP, SF, SFS, and RIN -- C C Set NPTS to the closest LMIN-1 nodes to K. C SUM = 0. NPTS(1) = KK DIST(1) = 0. LM1 = LMIN - 1 DO 1 LNP = 2,LM1 CALL GETNP (NCC,LCC,N,X,Y,LIST,LPTR,LEND, . LNP, NPTS,DIST, IERR) IF (IERR .NE. 0) GO TO 13 DS = DIST(LNP)**2 SUM = SUM + DS 1 CONTINUE C C Add additional nodes to NPTS until the relative increase C in DS is at least RTOL. C DO 3 LNP = LMIN,LMAX CALL GETNP (NCC,LCC,N,X,Y,LIST,LPTR,LEND, . LNP, NPTS,DIST, IERR) RS = DIST(LNP)**2 IF ((RS-DS)/DS .LE. RTOL) GO TO 2 IF (LNP .GT. 6) GO TO 4 2 SUM = SUM + RS 3 CONTINUE C C Use all LMAX nodes in the least squares fit. RS is C arbitrarily increased by 10 per cent. C RS = 1.1*RS LNP = LMAX + 1 C C There are LNP-2 equations corresponding to nodes NPTS(2), C ...,NPTS(LNP-1). C 4 SFS = REAL(LNP-2)/SUM SF = SQRT(SFS) RIN = 1./SQRT(RS) XK = X(KK) YK = Y(KK) ZK = Z(KK) C C A Q-R decomposition is used to solve the least squares C system. The transpose of the augmented regression C matrix is stored in A with columns (rows of A) defined C as follows: 1-3 are the quadratic terms, 4 and 5 are C the linear terms with coefficients DX and DY, and the C last column is the right hand side. C C Set up the first 5 equations and zero out the lower tri- C angle with Givens rotations. C DO 6 I = 1,5 NP = NPTS(I+1) W = 1./DIST(I+1) - RIN CALL SETRO1 (XK,YK,ZK,X(NP),Y(NP),Z(NP),SF,SFS, . W, A(1,I)) IF (I .EQ. 1) GO TO 6 DO 5 J = 1,I-1 JP1 = J + 1 L = 6 - J CALL GIVENS (A(J,J),A(J,I),C,S) CALL ROTATE (L,C,S,A(JP1,J),A(JP1,I)) 5 CONTINUE 6 CONTINUE C C Add the additional equations to the system using C the last column of A. I .LE. LNP. C I = 7 7 IF (I .LT. LNP) THEN NP = NPTS(I) W = 1./DIST(I) - RIN CALL SETRO1 (XK,YK,ZK,X(NP),Y(NP),Z(NP),SF,SFS, . W, A(1,6)) DO 8 J = 1,5 JP1 = J + 1 L = 6 - J CALL GIVENS (A(J,J),A(J,6),C,S) CALL ROTATE (L,C,S,A(JP1,J),A(JP1,6)) 8 CONTINUE I = I + 1 GO TO 7 ENDIF C C Test the system for ill-conditioning. C DMIN = MIN( ABS(A(1,1)),ABS(A(2,2)),ABS(A(3,3)), . ABS(A(4,4)),ABS(A(5,5)) ) IF (DMIN/W .GE. DTOL) GO TO 12 IF (LNP .LE. LMAX) THEN C C Add another node to the system and increase R. Note C that I = LNP. C LNP = LNP + 1 IF (LNP .LE. LMAX) THEN CALL GETNP (NCC,LCC,N,X,Y,LIST,LPTR,LEND, . LNP, NPTS,DIST, IERR) RS = DIST(LNP)**2 ENDIF RIN = 1./SQRT(1.1*RS) GO TO 7 ENDIF C C Stabilize the system by damping second partials -- add C multiples of the first three unit vectors to the first C three equations. C STF = W DO 11 I = 1,3 A(I,6) = STF DO 9 J = I+1,6 A(J,6) = 0. 9 CONTINUE DO 10 J = I,5 JP1 = J + 1 L = 6 - J CALL GIVENS (A(J,J),A(J,6),C,S) CALL ROTATE (L,C,S,A(JP1,J),A(JP1,6)) 10 CONTINUE 11 CONTINUE C C Test the damped system for ill-conditioning. C DMIN = MIN( ABS(A(4,4)),ABS(A(5,5)) ) IF (DMIN/W .LT. DTOL) GO TO 14 C C Solve the 2 by 2 triangular system for the partial C derivatives. C 12 DY = A(6,5)/A(5,5) DX = SF*(A(6,4) - A(5,4)*DY)/A(4,4) DY = SF*DY IER = LNP - 1 RETURN C C Invalid input parameter. C 13 IER = -1 RETURN C C No unique solution due to collinear nodes. C 14 IER = -2 RETURN END SUBROUTINE GRCOEF (SIGMA,DCUB, D,SD) REAL SIGMA, DCUB, D, SD C C*********************************************************** C C From SRFPACK C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C renka@cs.unt.edu C 11/18/96 C C This subroutine computes factors involved in the linear C system solved by Subroutines GRADG and SMSGS. C C On input: C C SIGMA = Nonnegative tension factor associated with a C triangulation arc. C C DCUB = Cube of the positive arc length. C C Input parameters are not altered by this routine. C C On output: C C D = Diagonal factor. D = SIG*(SIG*COSHM(SIG) - C SINHM(SIG))/(E*DCUB), where E = SIG*SINH(SIG) - C 2*COSHM(SIG). D > 0. C C SD = Off-diagonal factor. SD = SIG*SINHM(SIG)/ C (E*DCUB). SD > 0. C C SRFPACK module required by GRCOEF: SNHCSH C C Intrinsic function called by GRCOEF: EXP C C*********************************************************** C REAL COSHM, COSHMM, E, EMS, SCM, SIG, SINHM, SSINH, . SSM C SIG = SIGMA IF (SIG .LT. 1.E-9) THEN C C SIG = 0: cubic interpolant. C D = 4./DCUB SD = 2./DCUB ELSEIF (SIG .LE. .5) THEN C C 0 .LT. SIG .LE. .5: use approximations designed to avoid C cancellation error in the hyperbolic C functions when SIGMA is small. C CALL SNHCSH (SIG, SINHM,COSHM,COSHMM) E = (SIG*SINHM - COSHMM - COSHMM)*DCUB D = SIG*(SIG*COSHM-SINHM)/E SD = SIG*SINHM/E ELSE C C SIG > .5: scale SINHM and COSHM by 2*EXP(-SIG) in order C to avoid overflow when SIGMA is large. C EMS = EXP(-SIG) SSINH = 1. - EMS*EMS SSM = SSINH - 2.*SIG*EMS SCM = (1.-EMS)*(1.-EMS) E = (SIG*SSINH - SCM - SCM)*DCUB D = SIG*(SIG*SCM-SSM)/E SD = SIG*SSM/E ENDIF RETURN END SUBROUTINE INTRC0 (PX,PY,NCC,LCC,N,X,Y,Z,LIST,LPTR, . LEND, IST, PZ,IER) INTEGER NCC, LCC(*), N, LIST(*), LPTR(*), LEND(N), . IST, IER REAL PX, PY, X(N), Y(N), Z(N), PZ C C*********************************************************** C C From SRFPACK C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C renka@cs.unt.edu C 07/03/98 C C Given a triangulation of a set of nodes in the plane, C along with data values at the nodes, this subroutine com- C putes the value at P = (PX,PY) of the piecewise linear C function which interpolates the data values. The surface C is extended in a continuous fashion beyond the boundary of C the triangulation, allowing extrapolation. C C On input: C C PX,PY = Coordinates of the point P at which the sur- C face is to be evaluated. C C NCC = Number of constraint curves (refer to TRIPACK C Subroutine ADDCST). NCC .GE. 0. C C LCC = Array of length NCC (or dummy array of length C 1 if NCC = 0) containing the index of the C first node of constraint I in LCC(I). For I = C 1 to NCC, LCC(I+1)-LCC(I) .GE. 3, where C LCC(NCC+1) = N+1. C C N = Number of nodes in the triangulation. N .GE. 3. C C X,Y = Arrays of length N containing the coordinates C of the nodes with non-constraint nodes in the C first LCC(1)-1 locations, followed by NCC se- C quences of constraint nodes. C C Z = Array of length N containing data values at the C nodes. Refer to Subroutine ZGRADL. C C LIST,LPTR,LEND = Data structure defining the trian- C gulation. Refer to TRIPACK C Subroutine TRMESH. C C The above parameters are not altered by this routine. C C IST = Index of the starting node in the search for a C triangle containing P. 1 .LE. IST .LE. N. C The output value of IST from a previous call C may be a good choice. C C On output: C C IST = Index of one of the vertices of the triangle C containing P (or a boundary node which is vis- C ible from P) unless IER < 0. C C PZ = Value of the interpolatory surface at P, or C zero if IER < 0. C C IER = Error indicator: C IER = 0 if no errors were encountered and P is C contained in a triangle but not in a C constraint region. C IER = 1 if no errors were encountered and P C lies in a constraint region triangle. C PZ is effectively an extrapolated C value in this case. C IER = 2 if no errors were encountered and P is C exterior to the triangulation. PZ is C an extrapolated value in this case. C IER = -1 if NCC, N, or IST is outside its C valid range on input. LCC is not C tested for validity. C IER = -2 if the nodes are collinear. C C TRIPACK modules required by INTRC0: CRTRI, JRAND, LEFT, C LSTPTR, TRFIND C C SRFPACK module required by INTRC0: COORDS C C*********************************************************** C LOGICAL CRTRI INTEGER I1, I2, I3, IERR, LPL, N1, N2 REAL B1, B2, B3, DP, X1, X2, XP, Y1, Y2, YP C XP = PX YP = PY PZ = 0. C C Test for invalid input parameters. C IF (NCC .LT. 0 .OR. N .LT. 3 .OR. IST .LT. 1 . .OR. IST .GT. N) THEN IER = -1 RETURN ENDIF C C Find a triangle (I1,I2,I3) containing P, or a pair of C visible boundary nodes I1 and I2. C CALL TRFIND (IST,XP,YP,N,X,Y,LIST,LPTR,LEND, I1,I2,I3) IF (I1 .EQ. 0) THEN IER = -2 RETURN ENDIF IST = I1 IF (I3 .EQ. 0) GO TO 1 C C P is in a triangle. Compute its barycentric coordinates. C CALL COORDS (XP,YP,X(I1),X(I2),X(I3),Y(I1),Y(I2), . Y(I3), B1,B2,B3,IERR) IF (IERR .NE. 0) THEN IER = -2 RETURN ENDIF C C Compute an interpolated value. C PZ = B1*Z(I1) + B2*Z(I2) + B3*Z(I3) IER = 0 C IF (CRTRI(NCC,LCC,I1,I2,I3)) THEN IER = 1 ELSE IER = 0 ENDIF RETURN C C P is exterior to the triangulation. Extrapolate to P by C extending the interpolatory surface as a constant C beyond the boundary: PZ is the function value at Q C where Q is the closest boundary point to P. C C Determine Q by traversing the boundary starting from the C rightmost visible node I1. C 1 IER = 2 N2 = I1 C C Top of loop: C C Set N1 to the last neighbor of N2, and compute the dot C product DP = (N2->N1,N2->P). P FORWARD N2->N1 iff C DP > 0. C 2 LPL = LEND(N2) N1 = -LIST(LPL) X1 = X(N1) Y1 = Y(N1) X2 = X(N2) Y2 = Y(N2) DP = (X1-X2)*(XP-X2) + (Y1-Y2)*(YP-Y2) IF (DP .LE. 0.) THEN C C N2 is the closest boundary point to P. C PZ = Z(N2) RETURN ENDIF C C P FORWARD N2->N1. Test for P FORWARD N1->N2. C IF ((XP-X1)*(X2-X1) + (YP-Y1)*(Y2-Y1) .GT. 0.) THEN C C The closest boundary point to P lies on N2-N1. Compute C its local coordinates with respect to N2-N1. C B1 = DP/( (X2-X1)**2 + (Y2-Y1)**2 ) B2 = 1. - B1 PZ = B1*Z(N1) + B2*Z(N2) RETURN ENDIF C C Bottom of boundary traversal loop. C N2 = N1 GO TO 2 END SUBROUTINE INTRC1 (PX,PY,NCC,LCC,N,X,Y,Z,LIST,LPTR, . LEND,IFLGS,SIGMA,GRAD, . DFLAG, IST, PZ,PZX,PZY,IER) INTEGER NCC, LCC(*), N, LIST(*), LPTR(*), LEND(N), . IFLGS, IST, IER LOGICAL DFLAG REAL PX, PY, X(N), Y(N), Z(N), SIGMA(*), GRAD(2,N), . PZ, PZX, PZY C C*********************************************************** C C From SRFPACK C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C renka@cs.unt.edu C 07/03/98 C C Given a triangulation of a set of nodes in the plane, C along with data values and estimated gradients at the C nodes, this subroutine computes the value and, optionally, C the first partial derivatives at P = (PX,PY) of a C-1 C (once-continuously differentiable) function F which int- C erpolates the data values and gradients. Extrapolation to C a point exterior to the triangulation is accomplished by C extending the surface in such a way that F is C-1 over the C entire plane. C C Subroutine FVAL is used to evaluate an interpolatory C surface under tension, while Subroutine TVAL is called in C the case of no tension (IFLGS .LE. 0 and SIGMA(1) = 0). C However, the surface under tension is well-defined with C SIGMA = 0, and, in this case, the two interpolants are C identical on triangulation arcs and outside the triangula- C tion. The use of FVAL with no tension can be forced (at C a cost in efficiency) by setting IFLGS = 1 and storing C zeros in all components of SIGMA. Note, however, that C first partial derivatives are only available from TVAL C (and at points outside the triangulation); i.e., a proce- C dure for differentiating the surface under tension has not C been implemented. C C A set of interpolated values at the vertices of a rec- C tangular grid can be obtained by a single call to C Subroutine UNIF. Subroutine INTRC0 provides for evalua- C tion of the piecewise linear interpolatory surface. C C On input: C C PX,PY = Coordinates of the point P at which the sur- C face is to be evaluated. C C NCC = Number of constraint curves (refer to TRIPACK C Subroutine ADDCST). NCC .GE. 0. C C LCC = Array of length NCC (or dummy array of length C 1 if NCC = 0) containing the index of the C first node of constraint I in LCC(I). For I = C 1 to NCC, LCC(I+1)-LCC(I) .GE. 3, where C LCC(NCC+1) = N+1. C C N = Number of nodes in the triangulation. N .GE. 3. C C X,Y = Arrays of length N containing the coordinates C of the nodes with non-constraint nodes in the C first LCC(1)-1 locations, followed by NCC se- C quences of constraint nodes. C C Z = Array of length N containing data values at the C nodes. Refer to Subroutines ZGRADG and ZGRADL. C C LIST,LPTR,LEND = Data structure defining the trian- C gulation. Refer to TRIPACK C Subroutine TRMESH. C C IFLGS = Tension factor option: C IFLGS .LE. 0 if a single uniform tension C factor is to be used. C IFLGS .GE. 1 if variable tension is desired. C C SIGMA = Uniform tension factor (IFLGS .LE. 0), or C array containing tension factors associated C with arcs in one-to-one correspondence with C LIST entries (IFLGS .GE. 1). Refer to Sub- C routines FVAL, GETSIG, SIG0, SIG1, and SIG2. C C GRAD = 2 by N array whose columns contain estimated C gradients at the nodes with X partial deriva- C tives in the first row and Y partials in the C second. Refer to Subroutines GRADC, GRADG, C GRADL, SMSURF, ZGRADG, and ZGRADL. C C DFLAG = Logical flag which specifies whether first C partial derivatives at P are to be computed: C DFLAG = TRUE if and only if partials are C to be computed by TVAL. This option is only C valid for IFLGS .LE. 0 and SIGMA(1) = 0 (and C for points outside the triangulation). C C The above parameters are not altered by this routine. C C IST = Index of the starting node in the search for a C triangle containing P. 1 .LE. IST .LE. N. C The output value of IST from a previous call C may be a