/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:30:11 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "stgc1.h" #include void /*FUNCTION*/ stgc1( LOGICAL32 newtri, float s[][3], float *zout, LOGICAL32 wantdz, float dzout[]) { long int i, im1, ip1, j, m, mm1, mp1, _i, _r; static long int add1[3], sub1[3]; float a, b, c, da, db, dc, dgtild[3], dphi[3], dr[3], drho[3], fac, gtilda[3], lsq[3], phi[3], r[3], rho[3], sum; static int _aini = 1; /* OFFSET Vectors w/subscript range: 1 to dimension */ long *const Add1 = &add1[0] - 1; float *const Dgtild = &dgtild[0] - 1; float *const Dphi = &dphi[0] - 1; float *const Dr = &dr[0] - 1; float *const Drho = &drho[0] - 1; float *const Dzout = &dzout[0] - 1; float *const Gtilda = >ilda[0] - 1; float *const Lsq = &lsq[0] - 1; float *const Phi = &phi[0] - 1; float *const R = &r[0] - 1; float *const Rho = &rho[0] - 1; long *const Sub1 = &sub1[0] - 1; /* end of OFFSET VECTORS */ if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */ Add1[1] = 2; Add1[2] = 3; Add1[3] = 1; Sub1[1] = 3; Sub1[2] = 1; Sub1[3] = 2; _aini = 0; } /* Copyright (c) 1996 California Institute of Technology, Pasadena, CA. * ALL RIGHTS RESERVED. * Based on Government Sponsored Research NAS7-03001. *>> 1996-02-02 STGC1 CLL *>> 1995-09-26 STGC1 CLL Editing for inclusion into MATH77. * C.L.Lawson, JPL, 1976 Dec 7 * This subr interpolates over a triangle using the piecewise * bicubic formulas of Clough and Tocher (1965) as formulated for * computation by C.L.Lawson in JPL Tech. Memorandum 33-770, May,1976 * This method gives C1 continuity with neighboring triangles, i.e., * continuity of value and first partial derivatives. * Generally there will be jumps in second and higher order * derivatives across edges between triangles, and along implicit * edges that partition each triangle into three smaller triangles. * Optionally the subr also computes first partial derivatives. * * NEWTRI [in] This argument allows certain computations to be * skipped, saving 15 multiplies, 3 divides, and 9 adds, * when the user has made no changes to the contents of * columns 2, 3, 5, 6, 7, 8, & 9 of the S() array since a * previous call. (Thus it is allowable to have changed the * contents of columns 1 and 4.) * Setting NEWTRI = .TRUE. causes these computations to be done, * while .FALSE. cause them to be skipped. * * (S(1:3, 1:9) [inout] Columns 1-6 contain data set by the user to * specify the interpolation problem. Col 1 contains * unnormalized barycentric coordinates of the interpolation * point. The other cols contain data depending only on the * triangle and its vertex data and not on the interpolation * point. Cols 7-9 are computed by this subroutine, and in * appropriate cases can be saved from one call to the next. * See NEWTRI above. * * S( ,1) = UNNORMALIZED BARYCENTRIC COORDS OF INTERP POINT. * S( ,2) = U = X COORD OF EDGE VECTOR * S( ,3) = V = Y COORD OF EDGE VECTOR * S( ,4) = Z = FCN VALUE AT VERTEX * S( ,5) = ZX = PARTIAL DERIV OF Z W.R.T. X * S( ,6) = ZY = PARTIAL DERIV OF Z W.R.T. Y * S( ,7) = HTILDA = SCALED P.D. TANGENTIAL TO EDGE AT LEFT END. * S( ,8) = KTILDA = SCALED P.D. TANGENTIAL TO EDGE AT RIGHT END. * S( ,9) = LFAC = FACTOR INVOLVING LENGTHS OF EDGES. * * ZOUT [out] INTERPOLATED VALUE COMPUTED BY SUBR. * * 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 * INTERPOLATED SURFACE AT THE INTERPOLATION POINT. * ------------------------------------------------------------------ * Details of the contents of TRI(1:7) and S(1:3, 1:6). * * We assume the point Q is in the triangle indexed by INDTRI. * The indices of the vertices of this triangle, in counter- * clockwise, order are P(1), P(2), and P(3) which may be obtained as * P(i) = TRI(3+i) for i = 1, 2, and 3. TRI(7) contains the same * value as TRI(4). The triangle adjacent to this triangle across * the edge from P(i) to P(i+1) is indexed by TRI(i) for i = 1, 2, * and 3. If there is no adjacent triangle across this edge then * TRI(i) = 0. * * For descriptive convenience we regard the subscript of P() and the * 1st subsubscript of S(,) as always being reduced modulo 3 to 1, 2, * or 3. Also for convenience we shall write P(i) to mean the vertex * indexed by P(i). * * The unnormalized barycentric coordinate that is zero along the edge * from P(i) to P(i+1) and has a positive value at P(i+2) is stored in * S(i,1). The (x,y) coordinates of the vector from P(i) to P(i+1) are * stored in (S(i,2), S(i,3)). * The function value Z at P(i+2) is stored in S(i,4) for i = 1, 2, & 3. * ZX, the partial deriv of Z w.r.t. X at P(i+2) is stored in S(i,5) * for i = 1, 2, & 3. * ZY, the partial deriv of Z w.r.t. Y at P(i+2) is stored in S(i,6) * for i = 1, 2, & 3. * ------------------------------------------------------------------ *--S replaces "?": ?TGC1 * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ * COMPUTE R() = NORMALIZED BARYCENTRIC COORDINATES. * AND LSQ() = SQUARED EDGE LENGTHS * */ fac = 1.e00/(s[0][0] + s[0][1] + s[0][2]); for (i = 1; i <= 3; i++) { R[i] = fac*s[0][i - 1]; Lsq[i] = SQ(s[1][i - 1]) + SQ(s[2][i - 1]); } /* 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]; /* COMPUTE RHO() = CLOUGH-TOCHER PIECEWISE BICUBIC * CORRECTION FUNCTION * */ if (R[1] <= R[2]) { if (R[1] <= R[3]) { m = 1; } else { m = 3; } } else { if (R[2] <= R[3]) { m = 2; } else { m = 3; } } a = 0.5e0*SQ(R[m]); b = (1.0e0/3.0e0)*R[m]; c = Phi[m] + (5.0e0/3.0e0)*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); sum = 0.0e0; for (i = 1; i <= 3; i++) { ip1 = Add1[i]; im1 = Sub1[i]; if (newtri) { /* H TILDA */ s[6][i - 1] = s[1][i - 1]*s[4][ip1 - 1] + s[2][i - 1]* s[5][ip1 - 1]; /* K TILDA */ s[7][i - 1] = s[1][i - 1]*s[4][im1 - 1] + s[2][i - 1]* s[5][im1 - 1]; /* LFAC */ s[8][i - 1] = 3.0e0*(Lsq[ip1] - Lsq[im1])/Lsq[i]; } /* G TILDA */ Gtilda[i] = (R[ip1] - R[im1])*Phi[i] + s[8][i - 1]*Rho[i] - Rho[ip1] + Rho[im1]; sum += s[6][i - 1]*(Gtilda[i] + Phi[i]) + s[7][i - 1]*(Gtilda[i] - Phi[i]); } *zout = 0.5e0*sum; for (i = 1; i <= 3; i++) { ip1 = Add1[i]; im1 = Sub1[i]; *zout += s[3][i - 1]*(R[i] + Gtilda[im1] - Gtilda[ip1]); } /* FINISHED COMPUTING ZOUT. * * NOW COMPUTE DZOUT() IF REQUESTED. */ if (wantdz) { for (j = 1; j <= 2; j++) { if (j == 1) { for (i = 1; i <= 3; i++) { Dr[i] = -fac*s[2][i - 1]; } } else { for (i = 1; i <= 3; i++) { Dr[i] = fac*s[1][i - 1]; } } 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]; da = R[m]*Dr[m]; db = (1.e0/3.e0)*Dr[m]; dc = Dphi[m] + (5.e0/3.e0)*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); sum = 0.0e0; for (i = 1; i <= 3; i++) { ip1 = Add1[i]; im1 = Sub1[i]; Dgtild[i] = (R[ip1] - R[im1])*Dphi[i] + (Dr[ip1] - Dr[im1])*Phi[i] + s[8][i - 1]*Drho[i] - Drho[ip1] + Drho[im1]; sum += s[6][i - 1]*(Dgtild[i] + Dphi[i]) + s[7][i - 1]* (Dgtild[i] - Dphi[i]); } sum *= 0.5e0; for (i = 1; i <= 3; i++) { ip1 = Add1[i]; im1 = Sub1[i]; sum += s[3][i - 1]*(Dr[i] + Dgtild[im1] - Dgtild[ip1]); } Dzout[j] = sum; } } return; } /* end of function */