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;
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 occurs at (%d,%d) cell. " 234 "Z-value %f, zmin %f, zmax %f."),
235 l, k, zz, zmin, zmax);
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)
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)
int BM_get(struct BM *, int, int)
Gets 'val' from the bitmap.
double amin1(double, double)
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.