43 #include <grass/gis.h>
44 #include <grass/raster.h>
45 #include <grass/glocale.h>
46 #include <grass/bitmap.h>
48 #include <grass/interpf.h>
51 #define CEULER .57721566
68 double zmin,
double zmax,
69 double *zminac,
double *zmaxac,
70 double *gmin,
double *gmax,
71 double *c1min,
double *c1max,
72 double *c2min,
double *c2max,
83 double x_or = data->
x_orig;
84 double y_or = data->
y_orig;
89 static double *w2 =
NULL;
90 static double *w =
NULL;
93 double stepix, stepiy, xx, xg, yg, xx2;
94 double rfsta2, wm, dx, dy, dxx, dyy, dxy, h, bmgd1,
98 int ngstc, nszc, ngstr, nszr;
101 static int first_time_z = 1;
102 off_t offset, offset2;
103 double fstar2 = params->
fi * params->
fi / 4.;
104 double tfsta2, tfstad;
105 double ns_res, ew_res;
106 double rsin = 0, rcos = 0, teta, scale = 0;
110 teta = params->
theta / 57.295779;
119 ew_res = (((
struct quaddata *)(data))->xmax -
123 tfsta2 = (fstar2 * 2.) / dnorm;
124 tfstad = tfsta2 / dnorm;
130 stepix = ew_res / dnorm;
131 stepiy = ns_res / dnorm;
135 cond1 = ((params->
adx !=
NULL) || (params->
ady !=
NULL) || cond2);
138 if (!(w = (
double *)G_malloc(
sizeof(
double) * (params->
KMAX2 + 9)))) {
144 if (!(w2 = (
double *)G_malloc(
sizeof(
double) * (params->
KMAX2 + 9)))) {
153 ngstc = (int)(x_or / ew_res + 0.5) + 1;
154 nszc = ngstc + n_cols - 1;
155 ngstr = (int)(y_or / ns_res + 0.5) + 1;
156 nszr = ngstr + n_rows - 1;
159 for (k = ngstr; k <= nszr; k++) {
160 offset = offset1 * (k - 1);
161 yg = (k - ngstr) * stepiy + stepiy / 2.;
163 wm = yg - points[m - 1].
y;
167 for (l = ngstc; l <= nszc; l++) {
170 bmask =
BM_get(bitmask, l - 1, k - 1);
172 xg = (l - ngstc) * stepix + stepix / 2.;
183 xx = xg - points[m - 1].
x;
186 xxr = xx * rcos + w[m] * rsin;
187 yyr = w[m] * rcos - xx * rsin;
190 r2 = scale * xx2 + w2[m];
192 rfsta2 = scale * xx2 + w2[m];
198 rfsta2 = xx2 + w2[m];
201 h = h + b[m] * params->
interp(r, params->
fi);
206 dx = dx + bmgd1 * xx;
207 dy = dy + bmgd1 * w[m];
210 dxx = dxx + bmgd2 * xx2 + bmgd1;
211 dyy = dyy + bmgd2 * w2[m] + bmgd1;
212 dxy = dxy + bmgd2 * xx * w[m];
222 *zmaxac = *zminac = zz;
224 *zmaxac =
amax1(zz, *zmaxac);
225 *zminac =
amin1(zz, *zminac);
226 if ((zz > zmax + 0.1 * (zmax - zmin))
227 || (zz < zmin - 0.1 * (zmax - zmin))) {
232 G_warning(_(
"Overshoot - increase in tension suggested. "
233 "Overshoot occures at (%d,%d) cell. "
234 "Z-value %f, zmin %f, zmax %f."),
235 l, k, zz, zmin, zmax);
239 params->
az[
l] = (FCELL) zz;
242 params->
adx[
l] = (FCELL) (-dx * tfsta2);
243 params->
ady[
l] = (FCELL) (-dy * tfsta2);
245 params->
adxx[
l] = (FCELL) (-dxx * tfstad);
246 params->
adyy[
l] = (FCELL) (-dyy * tfstad);
247 params->
adxy[
l] = (FCELL) (-dxy * tfstad);
253 Rast_set_d_null_value(params->
az + l, 1);
257 Rast_set_d_null_value(params->
adx + l, 1);
258 Rast_set_d_null_value(params->
ady + l, 1);
260 Rast_set_d_null_value(params->
adxx + l, 1);
261 Rast_set_d_null_value(params->
adyy + l, 1);
262 Rast_set_d_null_value(params->
adxy + l, 1);
268 if (cond1 && (params->
deriv != 1)) {
269 if (params->
secpar(params, ngstc, nszc, k, bitmask,
270 gmin, gmax, c1min, c1max, c2min, c2max, cond1,
275 offset2 = (offset + ngstc - 1) *
sizeof(FCELL);
276 if (params->
wr_temp(params, ngstc, nszc, offset2) < 0)
int BM_get(struct BM *map, int x, int y)
Gets 'val' from the bitmap.
double amax1(double, double)
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, double *b, off_t offset1, double dnorm)
double amin1(double, double)
void G_warning(const char *msg,...)
Print a warning message to stderr.