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,
107 teta = params->
theta / M_R2D;
114 ns_res = (((
struct quaddata *)(data))->ymax -
115 ((
struct quaddata *)(data))->y_orig) /
117 ew_res = (((
struct quaddata *)(data))->xmax -
118 ((
struct quaddata *)(data))->x_orig) /
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);
241 params->
az[
l] = (FCELL)zz;
244 params->
adx[
l] = (FCELL)(-dx * tfsta2);
245 params->
ady[
l] = (FCELL)(-dy * tfsta2);
247 params->
adxx[
l] = (FCELL)(-dxx * tfstad);
248 params->
adyy[
l] = (FCELL)(-dyy * tfstad);
249 params->
adxy[
l] = (FCELL)(-dxy * tfstad);
254 Rast_set_d_null_value(params->
az +
l, 1);
259 Rast_set_d_null_value(params->
adx +
l, 1);
260 Rast_set_d_null_value(params->
ady +
l, 1);
262 Rast_set_d_null_value(params->
adxx +
l, 1);
263 Rast_set_d_null_value(params->
adyy +
l, 1);
264 Rast_set_d_null_value(params->
adxy +
l, 1);
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 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)