/*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 */