46 #include <grass/interpf.h>
48 #define CEULER .57721566
64 double zmin,
double zmax,
65 double *zminac,
double *zmaxac,
66 double *gmin,
double *gmax,
67 double *c1min,
double *c1max,
68 double *c2min,
double *c2max,
79 double x_or = data->
x_orig;
80 double y_or = data->
y_orig;
85 static double *w2 =
NULL;
86 static double *w =
NULL;
89 double stepix, stepiy, xx, xg, yg, xx2;
90 double wm, dx, dy, dxx, dyy, dxy, h, bmgd1,
94 int ngstc, nszc, ngstr, nszr;
97 static int first_time_z = 1;
98 off_t offset, offset2;
99 double fstar2 = params->
fi * params->
fi / 4.;
100 double tfsta2, tfstad;
101 double ns_res, ew_res;
102 double rsin = 0, rcos = 0, teta,
115 ((
struct quaddata *)(data))->y_orig) /
117 ew_res = (((
struct quaddata *)(data))->xmax -
122 tfsta2 = (fstar2 * 2.) / dnorm;
123 tfstad = tfsta2 / dnorm;
129 stepix = ew_res / dnorm;
130 stepiy = ns_res / dnorm;
134 cond1 = ((params->
adx !=
NULL) || (params->
ady !=
NULL) || cond2);
137 if (!(w = (
double *)
G_malloc(
sizeof(
double) * (params->
KMAX2 + 9)))) {
143 if (!(w2 = (
double *)
G_malloc(
sizeof(
double) * (params->
KMAX2 + 9)))) {
152 ngstc = (int)(x_or / ew_res + 0.5) + 1;
153 nszc = ngstc +
n_cols - 1;
154 ngstr = (int)(y_or / ns_res + 0.5) + 1;
155 nszr = ngstr +
n_rows - 1;
157 for (k = ngstr; k <= nszr; k++) {
158 offset = offset1 * (k - 1);
159 yg = (k - ngstr) * stepiy + stepiy / 2.;
165 for (
l = ngstc;
l <= nszc;
l++) {
172 xg = (
l - ngstc) * stepix +
187 xxr = xx * rcos + w[m] * rsin;
188 yyr = w[m] * rcos - xx * rsin;
191 r2 = scale * xx2 + w2[m];
207 dx = dx + bmgd1 * xx;
208 dy = dy + bmgd1 * w[m];
211 dxx = dxx + bmgd2 * xx2 + bmgd1;
212 dyy = dyy + bmgd2 * w2[m] + bmgd1;
213 dxy = dxy + bmgd2 * xx * w[m];
223 *zmaxac = *zminac = zz;
225 *zmaxac =
amax1(zz, *zmaxac);
226 *zminac =
amin1(zz, *zminac);
227 if ((zz > zmax + 0.1 * (zmax - zmin)) ||
228 (zz < zmin - 0.1 * (zmax - zmin))) {
234 _(
"Overshoot - increase in tension suggested. "
235 "Overshoot occurs at (%d,%d) cell. "
236 "Z-value %f, zmin %f, zmax %f."),
237 l, k, zz, zmin, zmax);
269 if (cond1 && (params->
deriv != 1)) {
270 if (params->
secpar(params, ngstc, nszc, k, bitmask, gmin, gmax,
271 c1min, c1max, c2min, c2max, cond1, cond2) < 0)
275 offset2 = (offset + ngstc - 1) *
sizeof(
FCELL);
276 if (params->
wr_temp(params, ngstc, nszc, offset2) < 0)
int BM_get(struct BM *, int, int)
Gets 'val' from the bitmap.
void G_warning(const char *,...) __attribute__((format(printf
void Rast_set_d_null_value(DCELL *, int)
To set a number of DCELL raster values to NULL.
#define UNUSED
A macro for an attribute, if attached to a variable, indicating that the variable is not used.
int IL_grid_calc_2d(struct interp_params *params, struct quaddata *data, struct BM *bitmask, double zmin, double zmax, double *zminac, double *zmaxac, double *gmin, double *gmax, double *c1min, double *c1max, double *c2min, double *c2max, double *ertot UNUSED, double *b, off_t offset1, double dnorm)
double amin1(double, double)
double amax1(double, double)