/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:32:04 */
/*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */
#include <math.h>
#include "fcrt.h"
#include "stgrec.h"
#include <stdlib.h>
/* PARAMETER translations */
#define MXWRK 28
#define ZERO 0.0e0
/* end of PARAMETER translations */
void /*FUNCTION*/ stgrec(
float x[],
float y[],
float z[],
float dz[][2],
long np,
long triang[],
long nt,
long bdry[][4],
long mb,
float xylim[],
long nx,
long ny,
float zfill,
float *zvals,
long mx,
long my,
long ncont,
LOGICAL32 wantpd,
float *dzvals)
{
#define ZVALS(I_,J_) (*(zvals+(I_)*(mx)+(J_)))
#define DZVALS(I_,J_,K_) (*(dzvals+(I_)*(my)*(mx)+(J_)*(mx)+(K_)))
LOGICAL32 first;
long int i, ix, iy, mode;
float delx, dely, dzout[2], svwrk0[MXWRK], svwrk1[MXWRK], xy[2];
/* OFFSET Vectors w/subscript range: 1 to dimension */
float *const Dzout = &dzout[0] - 1;
float *const Svwrk0 = &svwrk0[0] - 1;
float *const Svwrk1 = &svwrk1[0] - 1;
float *const X = &x[0] - 1;
float *const Xy = &xy[0] - 1;
float *const Xylim = &xylim[0] - 1;
float *const Y = &y[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.
*>> 1997-07-01 STGREC Krogh Reversed subscripts in B (CLL suggestion)
*>> 1997-06-18 STGREC CLL Added arguments Bdry and MB.
*>> 1996-02-02 STGREC CLL
*>> 1996-01-11 STGREC CLL
*>> 1995-11-03 STGREC CLL Added argument, NCONT.
*>> 1995-09-26 STGREC CLL Editing for inclusion into MATH77.
* This subr builds a rectangular array of Z values by interpo-
* lation in a triangular grid. The subr will use either C0 or C1
* interpolation methods, depending on the setting of NCONT = 0 or 1.
* Optionally this subr also builds a rectangular array of values of
* first partial derivatives.
* CLL, 1979 JUNE 12, MINOR CHANGES.
* C.L.LAWSON, JPL, 1976 DEC 8. ADDED PARTIAL DERIVS 1977 MAR 3.
* ------------------------------------------------------------------
* Subroutine arguments
*
* All args are of intent IN except ZVALS() and DZVALS() which are
* of intent OUT. The contents of DZ() are only used when NCONT = 1.
*
* X(),Y(),Z() [floating,in] NP sets of (X,Y,Z) data.
*
* (DZ(I,J),I=1,2),J=1,NP) [floating,in] Partial derivs of Z w.r.t.
* X and Y. The contents of this array are only used when NCONT = 1,
* not when NCONT = 0.
*
* NP [integer,in] No. of data points.
*
* TRIANG(1:6*NT),NT [integer,in] TRIANG() contains pointers (integers)
* describing the triangular grid. NT is the number of triangles.
* Each triangle requires 6 pointers.
*
* Bdry(1:4, MB) [integer, in] Array containing pointers defining the
* boundary of the (convex) triangular grid.
*
* MB [integer, in] Dimension parameter for Bdry() array.
*
* XYLIM(1:4) [floating,in] Contains XMIN, XMAX, YMIN, and YMAX
* specifying boundaries of the rectangular grid.
*
* NX,NY [integer, in] NO. OF X AND Y GRID LINES IN RECTANGULAR GRID.
*
* ZFILL [floating, in] Selects action to be taken when a node of the
* rectangular grid is outside the (convex) triangular grid.
* If ZFILL = 0.0, extrapolated values will be provided by this
* package.
* If ZFILL .ne. 0.0, this subroutine will store ZFILL into ZVALS()
* and DZVALS() for locations outside the triangular grid.
*
* ZVALS(,) [floating,out] MX BY MY ARRAY IN WHICH THE NX BY NY SET OF
* INTERPOLATED VALUES WILL BE STORED.
*
* MX,MY [integer, in] DIMENSIONS OF THE ARRAYS ZVALS( , ) AND
* DZVALS( , ,2)
*
* NCONT [integer,in] Order of continuity of interpolation method to be
* used.
* 0 selects a C0 method. This method uses linear interpolation over
* each triangle. The function value is continuous across the whole
* grid but generally the 1st partial derivatives are discontinuous
* across boundaries between triangles.
* 1 selects a C1 method. This method partitions each triangle into
* 3 smaller triangles and constructs a bicubic patch over each of
* the smaller triangles. The function value and 1st partial deriv-
* atives are continuous across the whole grid but generally the 2nd
* and higher order partial derivatives are discontinuous across
* boundaries between the original triangles and between the small
* triangles.
*
* WANTPD [logical, in] User sets = .true. if partial derivs are
* desired, and .false. otherwise.
*
* DZVALS(,,) [floating,out] MX by NY by 2 array in which the
* NX by NY by 2 set of partial derivatives will be stored if
* WANTPD = .TRUE.
* In the case of NCONT = 0, the interpolation function does not
* generally have a partial derivative at points on the boundary
* between two triangles. For such points this subr will return the
* 1st partial derivatives associated with one of the adjoining
* triangles.
* ------------------------------------------------------------------
*--S replaces "?": ?TGREC, ?TGFI
* ------------------------------------------------------------------ */
/* ------------------------------------------------------------------ */
delx = (Xylim[2] - Xylim[1])/(nx - 1);
dely = (Xylim[4] - Xylim[3])/(ny - 1);
for (i = 1; i <= MXWRK; i++)
{
Svwrk0[i] = ZERO;
}
/* XY(1) AND XY(2) WILL BE SET SUCCESSIVELY TO THE (X,Y) COORDS
* OF EACH POINT OF THE RECTANGULAR GRID.
* When NCONT .eq. 0, the values of the 1st derivatives are
* constant within each triangle. Thus in the case of
* NCONT .eq. 0 and WANTDZ being true, we only need to compute the
* partial derivatives when NEWTRI is true.
* */
Xy[2] = Xylim[3];
for (iy = 1; iy <= ny; iy++)
{
for (i = 1; i <= MXWRK; i++)
{
Svwrk1[i] = Svwrk0[i];
}
first = TRUE;
Xy[1] = Xylim[1];
for (ix = 1; ix <= nx; ix++)
{
stgfi( x, y, z, dz, triang, nt, bdry, mb, ncont, xy, &ZVALS(iy - 1,ix - 1),
wantpd, dzout, &mode, svwrk1 );
if (first)
{
for (i = 1; i <= MXWRK; i++)
{
Svwrk0[i] = Svwrk1[i];
}
first = FALSE;
}
if (mode == 0)
{
/* 0 means XY was found within the triangular grid.
* */
if (wantpd)
{
DZVALS(0,iy - 1,ix - 1) = Dzout[1];
DZVALS(1,iy - 1,ix - 1) = Dzout[2];
}
}
else if (mode == 1)
{
/* 1 MEANS XY IS OUTSIDE THE TRIANGULAR GRID.
* Will use extrapolated value if ZFILL .eq. zero.
* Otherwise set results to the ZFILL
* */
if (zfill == ZERO)
{
if (wantpd)
{
DZVALS(0,iy - 1,ix - 1) = Dzout[1];
DZVALS(1,iy - 1,ix - 1) = Dzout[2];
}
}
else
{
ZVALS(iy - 1,ix - 1) = zfill;
if (wantpd)
{
DZVALS(0,iy - 1,ix - 1) = zfill;
DZVALS(1,iy - 1,ix - 1) = zfill;
}
}
}
else
{
/* Mode = 2 MEANS CYCLING HAS HAPPENED IN the lookup
* process. This should not happen.
* ERROR MESSAGE WILL BE PRINTED FROM _TGFI.
* */
ZVALS(iy - 1,ix - 1) = zfill;
if (wantpd)
{
DZVALS(0,iy - 1,ix - 1) = zfill;
DZVALS(1,iy - 1,ix - 1) = zfill;
}
}
Xy[1] += delx;
}
Xy[2] += dely;
}
return;
#undef ZVALS
#undef DZVALS
} /* end of function */