19 #include <grass/N_pde.h>
37 "N_calc_gradient_field_2d_stats: compute gradient field stats");
52 field->
sum = sumx + sumy;
53 field->
nonull = nonullx + nonully;
114 double dx, dy, p1, p2, r1, r2, mean, grad, res;
120 (
"N_compute_gradient_field_2d: the arrays are not of equal size");
124 (
"N_compute_gradient_field_2d: the arrays are not of equal size");
128 (
"N_compute_gradient_field_2d: array sizes and geometry data are different");
131 G_debug(3,
"N_compute_gradient_field_2d: compute gradient field");
144 (
"N_compute_gradient_field_2d: gradient field sizes and geometry data are different");
148 for (j = 0; j < rows; j++)
149 for (i = 0; i < cols - 1; i++) {
158 grad = (p1 - p2) / dx;
173 for (j = 0; j < rows - 1; j++)
174 for (i = 0; i <
cols; i++) {
183 grad = (p1 - p2) / dy;
192 res = -1 * mean * grad;
256 G_fatal_error(
"N_compute_gradient_components_2d: x array is empty");
258 G_fatal_error(
"N_compute_gradient_components_2d: y array is empty");
264 if (x->
cols != cols || x->
rows != rows)
266 (
"N_compute_gradient_components_2d: the size of the x array doesn't fit the gradient field size");
267 if (y->
cols != cols || y->
rows != rows)
269 (
"N_compute_gradient_components_2d: the size of the y array doesn't fit the gradient field size");
271 for (j = 0; j < rows; j++)
272 for (i = 0; i <
cols; i++) {
276 if (grad.
WC == 0.0 || grad.
EC == 0.0)
277 vx = (grad.
WC + grad.
EC);
279 vx = (grad.
WC + grad.
EC) / 2;
280 if (grad.
NC == 0.0 || grad.
SC == 0.0)
281 vy = (grad.
NC + grad.
SC);
283 vy = (grad.
NC + grad.
SC) / 2;
302 double minx, miny, minz;
303 double maxx, maxy, maxz;
304 double sumx, sumy, sumz;
305 int nonullx, nonully, nonullz;
308 "N_calc_gradient_field_3d_stats: compute gradient field stats");
314 if (minx <= minz && minx <= miny)
316 if (miny <= minz && miny <= minx)
318 if (minz <= minx && minz <= miny)
321 if (maxx >= maxz && maxx >= maxy)
323 if (maxy >= maxz && maxy >= maxx)
325 if (maxz >= maxx && maxz >= maxy)
328 field->
sum = sumx + sumy + sumz;
329 field->
nonull = nonullx + nonully + nonullz;
395 int cols, rows, depths;
396 double dx, dy, dz, p1, p2, r1, r2, mean, grad, res;
403 (
"N_compute_gradient_field_3d: the arrays are not of equal size");
408 (
"N_compute_gradient_field_3d: the arrays are not of equal size");
413 (
"N_compute_gradient_field_3d: the arrays are not of equal size");
418 (
"N_compute_gradient_field_3d: array sizes and geometry data are different");
420 G_debug(3,
"N_compute_gradient_field_3d: compute gradient field");
429 if (gradfield ==
NULL) {
436 (
"N_compute_gradient_field_3d: gradient field sizes and geometry data are different");
439 for (k = 0; k < depths; k++)
440 for (j = 0; j < rows; j++)
441 for (i = 0; i < cols - 1; i++) {
450 grad = (p1 - p2) / dx;
462 "N_compute_gradient_field_3d: X-direction insert value %6.5g at %i %i %i ",
469 for (k = 0; k < depths; k++)
470 for (j = 0; j < rows - 1; j++)
471 for (i = 0; i <
cols; i++) {
480 grad = (p1 - p2) / dy;
489 res = -1 * mean * grad;
493 "N_compute_gradient_field_3d: Y-direction insert value %6.5g at %i %i %i ",
500 for (k = 0; k < depths - 1; k++)
501 for (j = 0; j < rows; j++)
502 for (i = 0; i <
cols; i++) {
511 grad = (p1 - p2) / dz;
523 "N_compute_gradient_field_3d: Z-direction insert value %6.5g at %i %i %i ",
585 int rows,
cols, depths;
594 G_fatal_error(
"N_compute_gradient_components_3d: x array is empty");
596 G_fatal_error(
"N_compute_gradient_components_3d: y array is empty");
598 G_fatal_error(
"N_compute_gradient_components_3d: z array is empty");
607 (
"N_compute_gradient_components_3d: the size of the x array doesn't fit the gradient field size");
610 (
"N_compute_gradient_components_3d: the size of the y array doesn't fit the gradient field size");
613 (
"N_compute_gradient_components_3d: the size of the z array doesn't fit the gradient field size");
615 for (k = 0; k < depths; k++)
616 for (j = 0; j < rows; j++)
617 for (i = 0; i <
cols; i++) {
620 if (grad.
WC == 0.0 || grad.
EC == 0.0)
621 vx = (grad.
WC + grad.
EC);
623 vx = (grad.
WC + grad.
EC) / 2;
624 if (grad.
NC == 0.0 || grad.
SC == 0.0)
625 vy = (grad.
NC + grad.
SC);
627 vy = (grad.
NC + grad.
SC) / 2;
628 if (grad.
TC == 0.0 || grad.
BC == 0.0)
629 vz = (grad.
TC + grad.
BC);
631 vz = (grad.
TC + grad.
BC) / 2;