/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:49 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "stgext.h" #include /* PARAMETER translations */ #define ONE 1.0e0 #define SIX 6.0e0 #define TWO 2.0e0 #define ZERO 0.0e0 /* end of PARAMETER translations */ void /*FUNCTION*/ stgext( float x[], float y[], float z[], float dz[][2], long triang[], long bdry[][4], long mb, long ncont, float q[], long indtri, long mode, float *zout, LOGICAL32 wantdz, float dzout[]) { long int i, ib, p0, p1, p2, tri[8]; float alpha, b, beta, bsq, c, derk0s, derk1s, fu0, fuu0, fuuv, fvu0, fvuv, g0, g1, h0s, h1s, k0s, k1s, mean, p0test, p1test, q0[2], q2[2], recipb, s, sarray[4][3], ss, st, t, tt, ugrad0, ugrad1, vgrad0, vgrad1, vq, w0[2], w0norm, w2[2], w2norm, xloc[3], yloc[3], z0, z2; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const Dzout = &dzout[0] - 1; float *const Q = &q[0] - 1; float *const Q0 = &q0[0] - 1; float *const Q2 = &q2[0] - 1; long *const Tri = &tri[0] - 1; long *const Triang = &triang[0] - 1; float *const W0 = &w0[0] - 1; float *const W2 = &w2[0] - 1; float *const X = &x[0] - 1; float *const Xloc = &xloc[0] - 1; float *const Y = &y[0] - 1; float *const Yloc = &yloc[0] - 1; float *const Z = &z[0] - 1; /* end of OFFSET VECTORS */ /* Copyright (c) 1996 California Institute of Technology, Pasadena, CA. * ALL RIGHTS RESERVED. * Based on Government Sponsored Research NAS7-03001. *>> 2000-12-01 STGEXT Krogh Removed unused parameters three & four. *>> 1997-07-01 STGEXT Krogh Reversed subscripts in B (CLL suggestion) *>> 1997-06-23 STGEXT Krogh Removed "_" in names, set up for M77CON. *>> 1997-06-22 cll *>> 1997-06-18 STGEXT CLL *>> 1996-10-15 STGEXT CLL * * This subroutine handles extrapolation outside the convex hull * of given data points. * * We partition the (infinite) space exterior to the convex hull by * constructing two rays into the exterior space rooted at each boundary * vertex, one ray perpendicular to each of the two boundary edges that * meet at the vertex. The boundary has NB vertices. This partitions * the exterior into 2*NB (infinite) cells. We call these cells wedges * and blocks. A wedge is a cell bounded by two boundary rays rooted at * a single vertex point. A block is a cell bounded by a boundary edge * and the two rays perpendicular to each end of the edge. * * We first determine which of these cells contains the specified point, * q, and then use an extrapolation formula in that cell that is * designed to maintain C0 continuity if Ncont = 0, or C1 continuity if * Ncont = 1. * * The C0 extrapolation is exact for a linear function but not for * polynomials of degree > 1. * The C1 extrapolation is exact for a linear function but not for * polynomials of degree > 1. * ------------------------------------------------------------------ * Subroutine arguments * * wantdz, dzout(1:2) * X(), Y() [in] (x,y) coordinates of vertices of the * triangular grid. * * Z() [in] Z(i) is the value at (X(i),Y(i)) of the data to * be interpolated. * * DZ(1:2, *) [in] DZ(1:2, i) are the values of the partial * derivatives of the interpolation function at (X(i),Y(i)) with * respect to x and y, respectively. * * Triang() [integer, in] Array of integer pointers defining the * connectivity of the triangular grid. */ /* Bdry(1:4, MB) [integer, in] Array containing pointers defining the * boundary of the (convex) triangular grid. * Bdry(1, K) = FWD POINTER. Points to next vertex in * counterclockwise order. * Bdry(2, K) = BACKWARD POINTER. Points to next vertex in * clockwise order. * Bdry(3, K) = A BOUNDARY POINT * Bdry(4, K) = A BOUNDARY TRIANGLE * The triangle Bdry(4, K) has a boundary edge that * connects the points Bdry(3, K) and Bdry(3, K+1). * In general not all elements of the array Bdry() are members of * the linked list defining the boundary. * The entry with K = 1 is a member of the boundary list. To scan * the boundary, start at K = 1 and follow the forward or backward * pointers. * * MB [in] * * NCONT [in] = 0 or 1 to request either C0 or C1 continuity. * * Q(1:2) [floating, in] The (x,y) coordinates of the point for which * this subr will compute an extrapolated value. * * indtri [integer, inout] On entry, this is the index of a boundary * triangle, relative to which Q is outside a boundary edge. * On return, this will be the index of a boundary triangle, * containing a boundary vertex that is a closest vertex on the * boundary to Q. This is being returned for use as a starting * point for the next search, if desired. * * MODE [integer, in] * * ZOUT [out] Extrapolated function value computed by this subroutine. * * WANTDZ [in] =.TRUE. means compute DZOUT() as well as ZOUT. * =.FALSE. means compute only ZOUT and not DZOUT(). * * DZOUT(1:2) [out] First partial derivs w.r.t. x and y of the * extrapolated surface at the point, Q(). * ------------------------------------------------------------------ * Internal Variables * * TRI(1:8) Integer array of pointers defining one triangle. * TRI(1:3) are pointers to neighboring triangles in counterclockwise * order. * TRI(4:6) are pointers to vertex points in counterclockwise order. * Subroutine _TGGET sets TRI(1:6). When needed in this subroutine * we set TRI(7) = TRI(4) and TRI(8) = TRI(5). * For i = 1, 2, or 3, the points TRI(i+3) and TRI(i+4) are endpoints * of the edge that is shared between this triangle and triangle * TRI(i). * If there is no adjacent triangle across this edge then TRI(i) = 0. */ /* p0, p1, p2, ib [integers] P0, p1 and p2 are indices of points. * Ib is the index into Bdry() such that p1 = Bdry(3, ib). * * If Q is found to be in a "block", the block is rooted on the * boundary segment connecting vertices p0 and p1 in the clockwise * direction. p0 = Bdry(3, Bdry(1, ib)), p1 = Bdry(3, ib). The * index of the boundary triangle containing points p0 and p1 is * Bdry(4, ib). P2 is not used in this case. * * If Q is found to be in a "wedge", the wedge is rooted at p1. * P0, p1, and ib have the same relations as described above for Q * in a block. In addition p2 is the next boundary vertex clockwise * from p1. Thus p2 = Bdry(3, Bdry(2, ib)). The index * of the boundary triangle containing points p1 and p2 is * Bdry(4, Bdry(2, ib)). * ------------------------------------------------------------------ *--S replaces "?": ?TGEXT, ?TGGET, ?TGQS, ?TGC0 * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ */ stgget( indtri, tri, triang ); Tri[7] = Tri[4]; p1 = Tri[mode + 3]; p0 = Tri[mode + 4]; /* Set ib to be the index of p1 in the boundary list. * */ ib = 1; for (i = 1; i <= mb; i++) { if (bdry[ib - 1][2] == p1) goto L_10; ib = bdry[ib - 1][0]; } /* call IERM1(SUBNAM,INDIC,LEVEL,MSG,LABEL,VALUE,FLAG) */ ierm1( "STGEXT", 1, 2, "Index p0 not found in boundary list.", "P0", p0, '.' ); *zout = ZERO; return; L_10: ; /* Define a (u,v) coordinate system, with p0 as the origin and the * positive u axis in the direction of p1. The point q is above * this u axis. * * Is q to the right of p0 in the (u,v) system? * */ p0test = (Q[1] - X[p0])*(X[p1] - X[p0]) + (Q[2] - Y[p0])*(Y[p1] - Y[p0]); if (p0test >= ZERO) { /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * * Searching to right (clockwise around boundary) * * Is q to the left of p1 in the (u,v) system? * */ L_20: ; p1test = (Q[1] - X[p1])*(X[p1] - X[p0]) + (Q[2] - Y[p1])*(Y[p1] - Y[p0]); if (p1test <= ZERO) { /* Q is in the block rooted on the segment from p0 to p1. * */ bsq = powif(X[p1] - X[p0],2) + powif(Y[p1] - Y[p0],2); s = ONE + p1test/bsq; goto L_200; } /* Move one segment to right around boundary. * New (u,v) system is based on new p0 and p1. * */ p0 = p1; ib = bdry[ib - 1][1]; p1 = bdry[ib - 1][2]; /* Is q to the left of p0 in the (u,v) system? * */ p0test = (Q[1] - X[p0])*(X[p1] - X[p0]) + (Q[2] - Y[p0])*(Y[p1] - Y[p0]); if (p0test <= ZERO) { /* Q is in the wedge rooted at p0. * Change names so wedge is rooted at p1. * */ p2 = p1; p1 = p0; ib = bdry[ib - 1][0]; p0 = bdry[bdry[ib - 1][0] - 1][2]; goto L_100; } goto L_20; } else { /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * * Searching to left (counterclockwise around boundary) * */ L_30: ; /* Move one segment to left around boundary. * New (u,v) system is based on new p0 and p1. * */ p1 = p0; ib = bdry[ib - 1][0]; p0 = bdry[bdry[ib - 1][0] - 1][2]; /* Is q to the right of p1 in the (u,v) system? * */ p1test = (Q[1] - X[p1])*(X[p1] - X[p0]) + (Q[2] - Y[p1])*(Y[p1] - Y[p0]); if (p1test >= ZERO) { /* Q is in the wedge rooted at p1. * */ p2 = bdry[bdry[ib - 1][1] - 1][2]; goto L_100; } /* Is q to the right of p0 in the (u,v) system? * */ p0test = (Q[1] - X[p0])*(X[p1] - X[p0]) + (Q[2] - Y[p0])*(Y[p1] - Y[p0]); if (p0test >= ZERO) { /* Q is in the block rooted on the segment from p0 to p1. * */ bsq = powif(X[p1] - X[p0],2) + powif(Y[p1] - Y[p0],2); s = p0test/bsq; goto L_200; /* Compute value ... */ } goto L_30; } /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ L_100: ; /* Here the point q is in a wedge rooted at p1 and * p1 = Bdry(3, ib). We also have p0 and p2 as the next boundary * vertices in the counterclockwise and clockwise directions, * respectively. * * Compute Zout, and optionally Dzout(1:2). * */ if (ncont == 0) { /* Let d0 be the vector from p1 to p0. Let w0 be the vector * resulting from rotating d0 90 degrees clockwise. * Let d2 be the vector from p1 to p2. Let w2 be the vector * resulting from rotating d2 90 degrees counterclockwise. * * Rescale w0 and w2 so they have the same length, say the * geometric mean of their original lengths. * Set q0 = p1+w0 and q2 = p1+w2. * * Let z0 be the z value obtained by using STGC0 to do C0 * extrapolation at q0 using the data of the boundary * triangle containing vertices p0 and p1. * * Let z2 be the z value obtained by using STGC0 to do C0 * extrapolation at q2 using the data of the boundary * triangle containing vertices p1 and p2. * * Construct a triangle with vertices p1, q2, and q0, * and give it vertex z values z(p1), z2 and z0, respectively. * * Use STGC0 with this triangle to interpolate/extrapolate * for a z value and/or for partial derivatives at Q. * */ W0[1] = Y[p0] - Y[p1]; W0[2] = -(X[p0] - X[p1]); w0norm = sqrtf( SQ(W0[1]) + SQ(W0[2]) ); W2[1] = -(Y[p2] - Y[p1]); W2[2] = X[p2] - X[p1]; w2norm = sqrtf( SQ(W2[1]) + SQ(W2[2]) ); mean = sqrtf( w0norm*w2norm ); Q0[1] = X[p1] + W0[1]*mean/w0norm; Q0[2] = Y[p1] + W0[2]*mean/w0norm; Q2[1] = X[p1] + W2[1]*mean/w2norm; Q2[2] = Y[p1] + W2[2]*mean/w2norm; stgget( bdry[ib - 1][3], tri, triang ); Tri[7] = Tri[4]; Tri[8] = Tri[5]; stgqs( q0, tri, x, y, sarray ); for (i = 1; i <= 3; i++) { sarray[3][i - 1] = Z[Tri[5 + i]]; } stgc0( sarray, &z0, FALSE, dzout ); stgget( bdry[bdry[ib - 1][1] - 1][3], tri, triang ); Tri[7] = Tri[4]; Tri[8] = Tri[5]; stgqs( q2, tri, x, y, sarray ); for (i = 1; i <= 3; i++) { sarray[3][i - 1] = Z[Tri[5 + i]]; } stgc0( sarray, &z2, FALSE, dzout ); Xloc[1] = X[p1]; Yloc[1] = Y[p1]; Xloc[2] = Q2[1]; Yloc[2] = Q2[2]; Xloc[3] = Q0[1]; Yloc[3] = Q0[2]; Tri[4] = 1; Tri[5] = 2; Tri[6] = 3; Tri[7] = 1; stgqs( q, tri, xloc, yloc, sarray ); sarray[3][0] = z0; sarray[3][1] = Z[p1]; sarray[3][2] = z2; stgc0( sarray, zout, wantdz, dzout ); } else { /* Here we assume Ncont = 1 * */ *zout = Z[p1] + dz[p1 - 1][0]*(Q[1] - X[p1]) + dz[p1 - 1][1]* (Q[2] - Y[p1]); if (wantdz) { Dzout[1] = dz[p1 - 1][0]; Dzout[2] = dz[p1 - 1][1]; } } goto L_800; /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ L_200: ; /* Here the point q is in a block rooted at the boundary edge * connecting p0 and p1. * Also, p0 = Bdry(3, Bdry(1, ib)), p1 = Bdry(3, ib) * Our goal is to compute Zout, and optionally Dzout(1:2). * * P1 is in the clockwise direction around the boundary from p0. * The (x,y) coordinates of p0 and p1 are (x(p0), y(p0)) and * (x(p1), y(p1)) respectively. * We use an (x',y') system which is the translate of the * (x,y) system to have p0 as its origin. Thus x' = x - x(p0) and * y' = y - y(p0). * We use a (u,v) coordinate system, which is a rotation of the * (x',y') system having its origin at p0 and the positive * u axis in the direction of p1. * * Let b denote the (positive) distance from p0 to p1. The (u,v) * coordinates of p0 and p1 are (0,0) and (0,b) respectively. * We already have bsq (= b**2) computed so we can compute * b = sqrt(bsq). * * Let uq and vq denote the u and v coordinates of q. Define * s = uq/b. We already have s computed, so we can compute * uq = b*s. * * In cases where needed we will compute alpha and/or beta. * These allow construction of the (orthogonal) matrix for * transformations between the (x',y') and (u,v) coordinate * systems. Alpha and beta are respectively the x' and y' * coordinates of p1. * * Let J denote the 2x2 orthogonal matrix: * J = |alpha -beta| * | beta alpha| * Then * J * |u| = |x'| * |v| |y'| * * J_transposed * |x'| = |u| * |y'| |v| * * J * |df/du| = |df/dx| * |df/dv| = |df/dy| * * J_transposed * |df/dx| = |df/du| * |df/dy| |df/dv| * */ if (ncont == 0) { stgget( bdry[ib - 1][3], tri, triang ); Tri[7] = Tri[4]; Tri[8] = Tri[5]; stgqs( q, tri, x, y, sarray ); for (i = 1; i <= 3; i++) { sarray[3][i - 1] = Z[Tri[5 + i]]; } stgc0( sarray, zout, wantdz, dzout ); /* zout = z(p0)*t + z(p1)*s * if(wantdz) then * ugrads = recipb * (z(p1) - z(p0)) * alpha = (x(p1) - x(p0)) * recipb * beta = (y(p1) - y(p0)) * recipb * dzout(1) = alpha*ugrads * dzout(2) = beta*ugrads * endif */ } else { /* Here we assume Ncont = 1 * */ b = sqrtf( bsq ); recipb = ONE/b; t = ONE - s; alpha = (X[p1] - X[p0])*recipb; beta = (Y[p1] - Y[p0])*recipb; vq = -beta*(Q[1] - X[p0]) + alpha*(Q[2] - Y[p0]); /* Test for and treat special cases of q at p0 or p1. * When q is at p0 we have vq and s both zero. * When q is at p1 we have vq and t both zero. * */ if (vq == ZERO) { if (s == ZERO) { *zout = Z[p0]; if (wantdz) { Dzout[1] = dz[p0 - 1][0]; Dzout[2] = dz[p0 - 1][1]; } goto L_800; } if (t == ZERO) { *zout = Z[p1]; if (wantdz) { Dzout[1] = dz[p1 - 1][0]; Dzout[2] = dz[p1 - 1][1]; } goto L_800; } } /* Partial derivs w.r.t. u and v at p0. * */ ugrad0 = alpha*dz[p0 - 1][0] + beta*dz[p0 - 1][1]; vgrad0 = -beta*dz[p0 - 1][0] + alpha*dz[p0 - 1][1]; /* Partial derivs w.r.t. u and v at p1. * */ ugrad1 = alpha*dz[p1 - 1][0] + beta*dz[p1 - 1][1]; vgrad1 = -beta*dz[p1 - 1][0] + alpha*dz[p1 - 1][1]; ss = SQ(s); tt = SQ(t); h0s = (ONE + TWO*s)*tt; h1s = (ONE + TWO*t)*ss; k0s = s*tt; k1s = -t*ss; fu0 = Z[p0]*h0s + Z[p1]*h1s + b*(ugrad0*k0s + ugrad1*k1s); fvu0 = t*vgrad0 + s*vgrad1; c = (vgrad1 - vgrad0)*recipb; g0 = c*vq/(b*s + vq); g1 = c*vq/(b*t + vq); *zout = fu0 + vq*fvu0 - b*vq*(k0s*g0 + k1s*g1); if (wantdz) { st = s*t; derk0s = tt - TWO*st; derk1s = ss - TWO*st; fuu0 = recipb*SIX*st*(Z[p1] - Z[p0]) + ugrad0*derk0s + ugrad1*derk1s; fuuv = fuu0 + c*vq - b*vq*(recipb*derk0s*g0 - k0s*g0/(b* s + vq) + recipb*derk1s*g1 + k1s*g1/(b*t + vq)); fvuv = fvu0 - b*(k0s*g0 + k1s*g1) - b*c*vq*(k0s*b*s/powif(b* s + vq,2) + k1s*b*t/powif(b*t + vq,2)); Dzout[1] = alpha*fuuv - beta*fvuv; Dzout[2] = beta*fuuv + alpha*fvuv; } } L_800: ; return; } /* end of function */ /* ================================================================== */ void /*FUNCTION*/ stgqs( float q[], long tri[], float x[], float y[], float s[][3]) { long int j, j1, j2; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const Q = &q[0] - 1; long *const Tri = &tri[0] - 1; float *const X = &x[0] - 1; float *const Y = &y[0] - 1; /* end of OFFSET VECTORS */ /*>> 1997-06-21 cll * Given point q() and the triangle described by tri(), this * subroutine puts values into s(1:3,1:3) to prepare for use of * subroutine STGC0 to evaluate the C0 interpolation/extrapolation * formulas. * The point q() is not required to be inside the triangle tri(). * This subroutine uses tri(4:7) * See comments in subroutines STGC0 or DTGFND for specifications of the * elements of tri() and s(). * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ */ for (j = 1; j <= 3; j++) { j1 = Tri[j + 3]; j2 = Tri[j + 4]; s[1][j - 1] = X[j2] - X[j1]; s[2][j - 1] = Y[j2] - Y[j1]; s[0][j - 1] = s[1][j - 1]*(Q[2] - Y[j1]) - s[2][j - 1]*(Q[1] - X[j1]); } return; } /* end of function */