SUBROUTINE DTGC1 (NEWTRI, S, ZOUT, WANTDZ, DZOUT) c Copyright (c) 1996 California Institute of Technology, Pasadena, CA. c ALL RIGHTS RESERVED. c Based on Government Sponsored Research NAS7-03001. C>> 1996-02-02 DTGC1 CLL C>> 1995-09-26 DTGC1 CLL Editing for inclusion into MATH77. c C.L.Lawson, JPL, 1976 Dec 7 c This subr interpolates over a triangle using the piecewise c bicubic formulas of Clough and Tocher (1965) as formulated for c computation by C.L.Lawson in JPL Tech. Memorandum 33-770, May,1976 c This method gives C1 continuity with neighboring triangles, i.e., c continuity of value and first partial derivatives. c Generally there will be jumps in second and higher order c derivatives across edges between triangles, and along implicit c edges that partition each triangle into three smaller triangles. c Optionally the subr also computes first partial derivatives. c c NEWTRI [in] This argument allows certain computations to be c skipped, saving 15 multiplies, 3 divides, and 9 adds, c when the user has made no changes to the contents of c columns 2, 3, 5, 6, 7, 8, & 9 of the S() array since a c previous call. (Thus it is allowable to have changed the c contents of columns 1 and 4.) c Setting NEWTRI = .TRUE. causes these computations to be done, c while .FALSE. cause them to be skipped. C c (S(1:3, 1:9) [inout] Columns 1-6 contain data set by the user to c specify the interpolation problem. Col 1 contains c unnormalized barycentric coordinates of the interpolation c point. The other cols contain data depending only on the c triangle and its vertex data and not on the interpolation c point. Cols 7-9 are computed by this subroutine, and in c appropriate cases can be saved from one call to the next. c See NEWTRI above. c c S( ,1) = UNNORMALIZED BARYCENTRIC COORDS OF INTERP POINT. c S( ,2) = U = X COORD OF EDGE VECTOR c S( ,3) = V = Y COORD OF EDGE VECTOR c S( ,4) = Z = FCN VALUE AT VERTEX c S( ,5) = ZX = PARTIAL DERIV OF Z W.R.T. X c S( ,6) = ZY = PARTIAL DERIV OF Z W.R.T. Y c S( ,7) = HTILDA = SCALED P.D. TANGENTIAL TO EDGE AT LEFT END. c S( ,8) = KTILDA = SCALED P.D. TANGENTIAL TO EDGE AT RIGHT END. c S( ,9) = LFAC = FACTOR INVOLVING LENGTHS OF EDGES. c c ZOUT [out] INTERPOLATED VALUE COMPUTED BY SUBR. C c WANTDZ [in] =.TRUE. MEANS COMPUTE DZOUT() AS WELL AS ZOUT. c =.FALSE. MEANS COMPUTE ONLY ZOUT AND NOT DZOUT(). C c DZOUT(1:2) [out] FIRST PARTIAL DERIVS W.R.T. X AND Y OF THE c INTERPOLATED SURFACE AT THE INTERPOLATION POINT. c ------------------------------------------------------------------ c Details of the contents of TRI(1:7) and S(1:3, 1:6). c c We assume the point Q is in the triangle indexed by INDTRI. c The indices of the vertices of this triangle, in counter- c clockwise, order are P(1), P(2), and P(3) which may be obtained as c P(i) = TRI(3+i) for i = 1, 2, and 3. TRI(7) contains the same c value as TRI(4). The triangle adjacent to this triangle across c the edge from P(i) to P(i+1) is indexed by TRI(i) for i = 1, 2, c and 3. If there is no adjacent triangle across this edge then c TRI(i) = 0. c c For descriptive convenience we regard the subscript of P() and the c 1st subsubscript of S(,) as always being reduced modulo 3 to 1, 2, c or 3. Also for convenience we shall write P(i) to mean the vertex c indexed by P(i). c c The unnormalized barycentric coordinate that is zero along the edge c from P(i) to P(i+1) and has a positive value at P(i+2) is stored in c S(i,1). The (x,y) coordinates of the vector from P(i) to P(i+1) are c stored in (S(i,2), S(i,3)). c The function value Z at P(i+2) is stored in S(i,4) for i = 1, 2, & 3. c ZX, the partial deriv of Z w.r.t. X at P(i+2) is stored in S(i,5) c for i = 1, 2, & 3. c ZY, the partial deriv of Z w.r.t. Y at P(i+2) is stored in S(i,6) c for i = 1, 2, & 3. c ------------------------------------------------------------------ c--D replaces "?": ?TGC1 C ------------------------------------------------------------------ INTEGER ADD1(3), i, im1, ip1, j, m, mm1, mp1, SUB1(3) DOUBLE PRECISION A, B, C double precision DA, DB, DC, DGTILD(3), DPHI(3) double precision DR(3), DRHO(3), DZOUT(2) double precision fac, GTILDA(3), LSQ(3) DOUBLE PRECISION r(3), RHO(3), PHI(3), S(3,9), sum DOUBLE PRECISION zout LOGICAL WANTDZ,NEWTRI DATA ADD1(1),ADD1(2),ADD1(3)/ 2, 3, 1 / DATA SUB1(1),SUB1(2),SUB1(3)/ 3, 1, 2 / C ------------------------------------------------------------------ c COMPUTE R() = NORMALIZED BARYCENTRIC COORDINATES. c AND LSQ() = SQUARED EDGE LENGTHS C FAC = 1.D00/(S(1,1)+S(2,1)+S(3,1)) DO 10 I=1,3 R(I) = FAC*S(I,1) LSQ(I)= S(I,2)**2 + S(I,3)**2 10 continue c COMPUTE PHI(I) = R(I+1) * R(I-1) PHI(1) = R(2)*R(3) PHI(2) = R(3)*R(1) PHI(3) = R(1)*R(2) C c COMPUTE RHO() = CLOUGH-TOCHER PIECEWISE BICUBIC c CORRECTION FUNCTION C IF ( R(1) .LE. R(2)) THEN IF ( R(1) .LE. R(3)) THEN M = 1 ELSE M = 3 END IF ELSE IF ( R(2) .LE. R(3)) THEN M = 2 ELSE M = 3 END IF END IF C A = 0.5D0 * R(M)**2 B = (1.0D0/3.0D0)*R(M) C = PHI(M) + (5.0D0/3.0D0)*A RHO(M) = R(M) * C - A MP1 = ADD1(M) MM1 = SUB1(M) RHO(MP1) = A*(R(MM1) - B) RHO(MM1) = A*(R(MP1) - B) C SUM = 0.0d0 DO 20 I=1,3 IP1 = ADD1(I) IM1 = SUB1(I) C IF (NEWTRI ) THEN c H TILDA S(I,7) = S(I,2)*S(IP1,5) + S(I,3)*S(IP1,6) c K TILDA S(I,8) = S(I,2)*S(IM1,5) + S(I,3)*S(IM1,6) c LFAC S(I,9) = 3.0D0*(LSQ(IP1)-LSQ(IM1)) / LSQ(I) END IF c G TILDA GTILDA(I)= (R(IP1)-R(IM1))*PHI(I) + S(I,9)*RHO(I) * -RHO(IP1)+RHO(IM1) SUM = SUM + S(I,7)*(GTILDA(I)+PHI(I)) * + S(I,8)*(GTILDA(I)-PHI(I)) 20 continue ZOUT = 0.5D0 * SUM C DO 30 I=1,3 IP1 = ADD1(I) IM1 = SUB1(I) ZOUT=ZOUT + S(I,4)*(R(I)+GTILDA(IM1)-GTILDA(IP1)) 30 continue c FINISHED COMPUTING ZOUT. C c NOW COMPUTE DZOUT() IF REQUESTED. IF ( WANTDZ) THEN DO 100 J=1,2 if(j .eq. 1) then DO 40 I=1,3 DR(I) = -FAC*S(I,3) 40 continue else DO 50 I=1,3 DR(I) = +FAC*S(I,2) 50 continue endif C DPHI(1) = R(2)*DR(3) + DR(2)*R(3) DPHI(2) = R(3)*DR(1) + DR(3)*R(1) DPHI(3) = R(1)*DR(2) + DR(1)*R(2) C DA = R(M)*DR(M) DB = (1.D0/3.D0)*DR(M) DC = DPHI(M) + (5.D0/3.D0)*DA DRHO(M) = R(M)*DC + DR(M)*C - DA DRHO(MP1) = A*(DR(MM1) - DB) + DA*(R(MM1) - B) DRHO(MM1) = A*(DR(MP1) - DB) + DA*(R(MP1) - B) C SUM = 0.0d0 DO 60 I=1,3 IP1=ADD1(I) IM1=SUB1(I) DGTILD(I) =(R(IP1)-R(IM1))*DPHI(I) * +(DR(IP1) - DR(IM1))*PHI(I) * +S(I,9)*DRHO(I) - DRHO(IP1) + DRHO(IM1) SUM = SUM +S(I,7)*(DGTILD(I)+DPHI(I)) * +S(I,8)*(DGTILD(I)-DPHI(I)) 60 continue SUM = 0.5D0 * SUM DO 70 I=1,3 IP1=ADD1(I) IM1=SUB1(I) SUM = SUM + S(I,4)*(DR(I) + DGTILD(IM1) - DGTILD(IP1)) 70 continue DZOUT(J) = SUM 100 continue END IF RETURN END