34 #include <grass/gis.h>
35 #include <grass/interpf.h>
36 #include <grass/gmath.h>
62 int n1, k1, k2, k, i1,
l, m, i, j;
63 double fstar2 = params->
fi * params->
fi / 4.;
64 static double *A =
NULL;
66 double rsin = 0, rcos = 0, teta, scale = 0;
70 teta = params->
theta * (M_PI / 180);
84 fprintf(stderr,
"Cannot allocate memory for A\n");
94 for (k = 1; k <= n_points; k++) {
103 for (k = 1; k <= n_points; k++) {
107 if (params->
rsm < 0.) {
108 A[i1] = -points[k - 1].
sm;
117 for (l = k2; l <= n_points; l++) {
118 xx = points[k - 1].
x - points[l - 1].
x;
119 yy = points[k - 1].
y - points[l - 1].
y;
123 xxr = xx * rcos + yy * rsin;
124 yyr = yy * rcos - xx * rsin;
127 r = scale * xx * xx + yy * yy;
128 rfsta2 = fstar2 * (scale * xx * xx + yy * yy);
131 r = xx * xx + yy * yy;
132 rfsta2 = fstar2 * (xx * xx + yy * yy);
136 fprintf(stderr,
"ident. points in segm.\n");
137 fprintf(stderr,
"x[%d]=%f, x[%d]=%f, y[%d]=%f, y[%d]=%f\n",
138 k - 1, points[k - 1].x, l - 1, points[l - 1].x, k - 1,
139 points[k - 1].y, l - 1, points[l - 1].y);
143 A[i1] = params->
interp(r, params->
fi);
149 for (k = 1; k <= n1; k++) {
152 for (l = k2; l <= n1; l++) {
153 m = (l - 1) * n1 + k;
155 amaxa =
amax1(A[m], amaxa);
159 for (i = 0; i <= n_points; i++) {
160 for (j = 0; j <= n_points; j++) {
166 G_debug(3,
"calling G_ludcmp() n=%d indx=%d", n_points, *indx);
167 if (
G_ludcmp(matrix, n_points + 1, indx, &d) <= 0) {
169 fprintf(stderr,
"G_ludcmp() failed! n=%d d=%.2f\n", n_points, d);
int IL_matrix_create(struct interp_params *params, struct triple *points, int n_points, double **matrix, int *indx)
Creates system of linear equations from interpolated points.
double amax1(double, double)
int G_ludcmp(double **a, int n, int *indx, double *d)
LU decomposition.
int G_debug(int level, const char *msg,...)
Print debugging message.
double * G_alloc_vector(size_t n)
Vector matrix memory allocation.