35 #include <grass/interpf.h> 45 static double *A =
NULL;
51 fprintf(stderr,
"Cannot allocate memory for A\n");
82 int n1, k1, k2, k, i1,
l, m, i, j;
83 double fstar2 = params->
fi * params->
fi / 4.;
85 double rsin = 0, rcos = 0, teta, scale = 0;
104 for (k = 1; k <= n_points; k++) {
113 for (k = 1; k <= n_points; k++) {
117 if (params->
rsm < 0.) {
118 A[i1] = -points[k - 1].
sm;
127 for (l = k2; l <= n_points; l++) {
128 xx = points[k - 1].
x - points[l - 1].
x;
129 yy = points[k - 1].
y - points[l - 1].
y;
133 xxr = xx * rcos + yy * rsin;
134 yyr = yy * rcos - xx * rsin;
137 r = scale * xx * xx + yy * yy;
138 rfsta2 = fstar2 * (scale * xx * xx + yy * yy);
141 r = xx * xx + yy * yy;
142 rfsta2 = fstar2 * (xx * xx + yy * yy);
146 fprintf(stderr,
"ident. points in segm.\n");
147 fprintf(stderr,
"x[%d]=%f, x[%d]=%f, y[%d]=%f, y[%d]=%f\n",
148 k - 1, points[k - 1].
x, l - 1, points[l - 1].
x, k - 1,
149 points[k - 1].y, l - 1, points[l - 1].y);
153 A[i1] = params->
interp(r, params->
fi);
159 for (k = 1; k <= n1; k++) {
162 for (l = k2; l <= n1; l++) {
163 m = (l - 1) * n1 + k;
165 amaxa =
amax1(A[m], amaxa);
169 for (i = 0; i <= n_points; i++) {
170 for (j = 0; j <= n_points; j++) {
176 G_debug(3,
"calling G_ludcmp() n=%d indx=%d", n_points, *indx);
177 if (
G_ludcmp(matrix, n_points + 1, indx, &d) <= 0) {
179 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)
double amax1(double, double)
double * G_alloc_vector(size_t)
Vector matrix memory allocation.
int G_ludcmp(double **, int, int *, double *)
LU decomposition.
int IL_matrix_create_alloc(struct interp_params *params, struct triple *points, int n_points, double **matrix, int *indx, double *A)
Creates system of linear equations from interpolated points.
int G_debug(int, const char *,...) __attribute__((format(printf