📄 shepard.cpp
字号:
/* RSQ = ARRAY OF LENGTH N CONTAINING THE SQUARED RADII */
/* WHICH ENTER INTO THE WEIGHTS DEFINING Q. */
/* */
/* A = 5 BY N ARRAY CONTAINING THE COEFFICIENTS FOR THE */
/* NODAL FUNCTIONS DEFINING Q. */
/* */
/* INPUT PARAMETERS ARE NOT ALTERED BY THIS FUNCTION. THE */
/* PARAMETERS OTHER THAN PX AND PY SHOULD BE INPUT UNALTERED */
/* FROM THEIR VALUES ON OUTPUT FROM QSHEP2. THIS FUNCTION */
/* SHOULD NOT BE CALLED IF A NONZERO ERROR FLAG WAS RETURNED */
/* BY QSHEP2. */
/* */
/* ON OUTPUT -- */
/* */
/* QS2VAL = FUNCTION VALUE Q(PX, PY) UNLESS N, NR, DX, */
/* DY, OR RMAX IS INVALID, IN WHICH CASE NO */
/* VALUE IS RETURNED. */
/* */
/* MODULES REQUIRED BY QS2VAL -- NONE */
/* */
/* INTRINSIC FUNCTIONS CALLED BY QS2VAL -- IFIX, SQRT */
/* */
/****************************************************************/
/* Parameter adjustments */
a -= 6;
--rsq;
--lnext;
--f;
--y;
--x;
lcell_dim1 = *nr;
lcell_offset = 1 + lcell_dim1;
lcell -= lcell_offset;
/* Function Body */
xp = *px;
yp = *py;
if (*n < 6 || *nr < 1 || *dx <= 0.f || *dy <= 0.f || *rmax < 0.f)
{
return _missing_;
}
/* SET IMIN, IMAX, JMIN, AND JMAX TO CELL INDICES DEFINING */
/* THE RANGE OF THE SEARCH FOR NODES WHOSE RADII INCLUDE */
/* P. THE CELLS WHICH MUST BE SEARCHED ARE THOSE INTER- */
/* SECTED BY(OR CONTAINED IN) A CIRCLE OF RADIUS RMAX */
/* CENTERED AT P. */
imin = (int)((xp - *xmin - *rmax) / *dx) + 1;
imax = (int)((xp - *xmin + *rmax) / *dx) + 1;
if (imin < 1)
{
imin = 1;
}
if (imax > *nr)
{
imax = *nr;
}
jmin = (int)((yp - *ymin - *rmax) / *dy) + 1;
jmax = (int)((yp - *ymin + *rmax) / *dy) + 1;
if (jmin < 1)
{
jmin = 1;
}
if (jmax > *nr)
{
jmax = *nr;
}
/* THE FOLLOWING IS A TEST FOR NO CELLS WITHIN THE CIRCLE */
/* OF RADIUS RMAX. */
if (imin > imax || jmin > jmax)
{
goto L5;
}
/* ACCUMULATE WEIGHT VALUES IN SW AND WEIGHTED NODAL FUNCTION */
/* VALUES IN SWQ. THE WEIGHTS ARE W(K) =((R - D)+/(R*D))**2 */
/* FOR R**2 = RSQ(K) AND D = DISTANCE BETWEEN P AND NODE K. */
sw = 0.f;
swq = 0.f;
/* OUTER LOOP ON CELLS(I, J). */
i__1 = jmax;
for (j = jmin; j <= i__1; ++j)
{
i__2 = imax;
for (i__ = imin; i__ <= i__2; ++i__)
{
k = lcell[i__ + j * lcell_dim1];
if (k == 0)
{
goto L3;
}
/* INNER LOOP ON NODES K. */
L1:
delx = xp - x[k];
dely = yp - y[k];
dxsq = delx * delx;
dysq = dely * dely;
ds = dxsq + dysq;
rs = rsq[k];
if (ds >= rs)
{
goto L2;
}
if (ds == 0.f)
{
goto L4;
}
rds = rs * ds;
rd = sqrt(rds);
w = (rs + ds - rd - rd) / rds;
sw += w;
swq += w * (a[k * 5 + 1] * dxsq + a[k * 5 + 2] * delx * dely + a[
k * 5 + 3] * dysq + a[k * 5 + 4] * delx + a[k * 5 + 5] *
dely + f[k]);
/* BOTTOM OF LOOP ON NODES IN CELL(I, J). */
L2:
kp = k;
k = lnext[kp];
if (k != kp)
{
goto L1;
}
L3:
;
}
}
/* SW = 0 IFF P IS NOT WITHIN THE RADIUS R(K) FOR ANY NODE K. */
if (sw == 0.f)
{
goto L5;
}
ret_val = swq / sw;
return ret_val;
/*(PX, PY) =(X(K), Y(K)) */
L4:
ret_val = f[k];
return ret_val;
/* ALL WEIGHTS ARE 0 AT P. */
L5:
ret_val = _missing_;
return ret_val;
}
int qs2grd_(double *px, double *py, int *n, double *x, double *y,
double *f, int *nr, int *lcell, int *lnext, double *xmin,
double *ymin, double *dx, double *dy, double *rmax, double *rsq, double *a, double *
q, double *qx, double *qy, int *ier)
{
int lcell_dim1, lcell_offset, i__1, i__2;
// double sqrt(double);
static int i__, j, k;
static double t, w, rd, ds;
static int kp;
static double qk, rs, xp, yp, sw, wx, wy, rds, qkx, qky, swq, sws, swx, swy;
static int imin, jmin, imax, jmax;
static double delx, dely, dxsq, dysq, swqx, swqy;
/****************************************************************/
/* */
/* ROBERT RENKA */
/* UNIV. OF NORTH TEXAS */
/*(817) 565 - 2767 */
/* 10/28/87 */
/* */
/* THIS SUBROUTINE COMPUTES THE VALUE AND GRADIENT AT */
/*(PX, PY) OF THE INTERPOLATORY FUNCTION Q DEFINED IN SUB- */
/* ROUTINE QSHEP2. Q(X, Y) IS A WEIGHTED SUM OF QUADRATIC */
/* NODAL FUNCTIONS. */
/* */
/* ON INPUT -- */
/* */
/* PX, PY = CARTESIAN COORDINATES OF THE POINT AT WHICH */
/* Q AND ITS PARTIALS ARE TO BE EVALUATED. */
/* */
/* N = NUMBER OF NODES AND DATA VALUES DEFINING Q. */
/* N .GE. 6. */
/* */
/* X, Y, F = ARRAYS OF LENGTH N CONTAINING THE NODES AND */
/* DATA VALUES INTERPOLATED BY Q. */
/* */
/* NR = NUMBER OF ROWS AND COLUMNS IN THE CELL GRID. */
/* REFER TO STORE2. NR .GE. 1. */
/* */
/* LCELL = NR BY NR ARRAY OF NODAL INDICES ASSOCIATED */
/* WITH CELLS. REFER TO STORE2. */
/* */
/* LNEXT = ARRAY OF LENGTH N CONTAINING NEXT - NODE INDI- */
/* CES. REFER TO STORE2. */
/* */
/* XMIN, YMIN, DX, DY = MINIMUM NODAL COORDINATES AND CELL */
/* DIMENSIONS. DX AND DY MUST BE */
/* POSITIVE. REFER TO STORE2. */
/* */
/* RMAX = SQUARE ROOT OF THE LARGEST ELEMENT IN RSQ -- */
/* MAXIMUM RADIUS. */
/* */
/* RSQ = ARRAY OF LENGTH N CONTAINING THE SQUARED RADII */
/* WHICH ENTER INTO THE WEIGHTS DEFINING Q. */
/* */
/* A = 5 BY N ARRAY CONTAINING THE COEFFICIENTS FOR THE */
/* NODAL FUNCTIONS DEFINING Q. */
/* */
/* INPUT PARAMETERS ARE NOT ALTERED BY THIS SUBROUTINE. */
/* THE PARAMETERS OTHER THAN PX AND PY SHOULD BE INPUT UNAL- */
/* TERED FROM THEIR VALUES ON OUTPUT FROM QSHEP2. THIS SUB- */
/* ROUTINE SHOULD NOT BE CALLED IF A NONZERO ERROR FLAG WAS */
/* RETURNED BY QSHEP2. */
/* */
/* ON OUTPUT -- */
/* */
/* Q = VALUE OF Q AT(PX, PY) UNLESS IER .EQ. 1, IN */
/* WHICH CASE NO VALUES ARE RETURNED. */
/* */
/* QX, QY = FIRST PARTIAL DERIVATIVES OF Q AT(PX, PY) */
/* UNLESS IER .EQ. 1. */
/* */
/* IER = ERROR INDICATOR */
/* IER = 0 IF NO ERRORS WERE ENCOUNTERED. */
/* IER = 1 IF N, NR, DX, DY OR RMAX IS INVALID. */
/* IER = 2 IF NO ERRORS WERE ENCOUNTERED BUT */
/*(PX, PY) IS NOT WITHIN THE RADIUS R(K) */
/* FOR ANY NODE K(AND THUS Q = QX = QY = 0). */
/* */
/* MODULES REQUIRED BY QS2GRD -- NONE */
/* */
/* INTRINSIC FUNCTIONS CALLED BY QS2GRD -- IFIX, SQRT */
/* */
/* **************************************************************/
a -= 6;
--rsq;
--lnext;
--f;
--y;
--x;
lcell_dim1 = *nr;
lcell_offset = 1 + lcell_dim1;
lcell -= lcell_offset;
xp = *px;
yp = *py;
if (*n < 6 || *nr < 1 || *dx <= 0.f || *dy <= 0.f || *rmax < 0.f)
{
goto L5;
}
/* SET IMIN, IMAX, JMIN, AND JMAX TO CELL INDICES DEFINING */
/* THE RANGE OF THE SEARCH FOR NODES WHOSE RADII INCLUDE */
/* P. THE CELLS WHICH MUST BE SEARCHED ARE THOSE INTER- */
/* SECTED BY(OR CONTAINED IN) A CIRCLE OF RADIUS RMAX */
/* CENTERED AT P. */
imin = (int)((xp - *xmin - *rmax) / *dx) + 1;
imax = (int)((xp - *xmin + *rmax) / *dx) + 1;
if (imin < 1)
{
imin = 1;
}
if (imax > *nr)
{
imax = *nr;
}
jmin = (int)((yp - *ymin - *rmax) / *dy) + 1;
jmax = (int)((yp - *ymin + *rmax) / *dy) + 1;
if (jmin < 1)
{
jmin = 1;
}
if (jmax > *nr)
{
jmax = *nr;
}
/* THE FOLLOWING IS A TEST FOR NO CELLS WITHIN THE CIRCLE */
/* OF RADIUS RMAX. */
if (imin > imax || jmin > jmax)
{
goto L6;
}
/* Q = SWQ/SW = SUM(W(K)*Q(K))/SUM(W(K)) WHERE THE SUM IS */
/* FROM K = 1 TO N, Q(K) IS THE QUADRATIC NODAL FUNCTION, */
/* AND W(K) =((R - D)+/(R*D))**2 FOR RADIUS R(K) AND DIST- */
/* ANCE D(K). THUS */
/* QX =(SWQX*SW - SWQ*SWX)/SW**2 AND */
/* QY =(SWQY*SW - SWQ*SWY)/SW**2 */
/* WHERE SWQX AND SWX ARE PARTIAL DERIVATIVES WITH RESPECT */
/* TO X OF SWQ AND SW, RESPECTIVELY. SWQY AND SWY ARE DE- */
/* FINED SIMILARLY. */
sw = 0.f;
swx = 0.f;
swy = 0.f;
swq = 0.f;
swqx = 0.f;
swqy = 0.f;
/* OUTER LOOP ON CELLS(I, J). */
i__1 = jmax;
for (j = jmin; j <= i__1; ++j)
{
i__2 = imax;
for (i__ = imin; i__ <= i__2; ++i__)
{
k = lcell[i__ + j * lcell_dim1];
if (k == 0)
{
goto L3;
}
/* INNER LOOP ON NODES K. */
L1:
delx = xp - x[k];
dely = yp - y[k];
dxsq = delx * delx;
dysq = dely * dely;
ds = dxsq + dysq;
rs = rsq[k];
if (ds >= rs)
{
goto L2;
}
if (ds == 0.f)
{
goto L4;
}
rds = rs * ds;
rd = sqrt(rds);
w = (rs + ds - rd - rd) / rds;
t = (rd - rs) * 2.f / (ds * rds);
wx = delx * t;
wy = dely * t;
qkx = a[k * 5 + 1] * 2.f * delx + a[k * 5 + 2] * dely;
qky = a[k * 5 + 2] * delx + a[k * 5 + 3] * 2.f * dely;
qk = (qkx * delx + qky * dely) / 2.f;
qkx += a[k * 5 + 4];
qky += a[k * 5 + 5];
qk = qk + a[k * 5 + 4] * delx + a[k * 5 + 5] * dely + f[k];
sw += w;
swx += wx;
swy += wy;
swq += w * qk;
swqx = swqx + wx * qk + w * qkx;
swqy = swqy + wy * qk + w * qky;
/* BOTTOM OF LOOP ON NODES IN CELL(I, J). */
L2:
kp = k;
k = lnext[kp];
if (k != kp)
{
goto L1;
}
L3:
;
}
}
/* SW = 0 IFF P IS NOT WITHIN THE RADIUS R(K) FOR ANY NODE K. */
if (sw == 0.f)
{
goto L6;
}
*q = swq / sw;
sws = sw * sw;
*qx = (swqx * sw - swq * swx) / sws;
*qy = (swqy * sw - swq * swy) / sws;
*ier = 0;
return 0;
/*(PX, PY) =(X(K), Y(K)) */
L4:
*q = f[k];
*qx = a[k * 5 + 4];
*qy = a[k * 5 + 5];
*ier = 0;
return 0;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -