GRASS GIS 7 Programmer's Manual  7.0.2(2015)-r00000
matrix.c
Go to the documentation of this file.
1 
31 #include <stdio.h>
32 #include <math.h>
33 #include <unistd.h>
34 #include <grass/gis.h>
35 #include <grass/interpf.h>
36 #include <grass/gmath.h>
37 
38 
53 int IL_matrix_create(struct interp_params *params,
54  struct triple *points, /* points for interpolation */
55  int n_points, /* number of points */
56  double **matrix, /* matrix */
57  int *indx)
58 {
59  double xx, yy;
60  double rfsta2, r;
61  double d;
62  int n1, k1, k2, k, i1, l, m, i, j;
63  double fstar2 = params->fi * params->fi / 4.;
64  static double *A = NULL;
65  double RO, amaxa;
66  double rsin = 0, rcos = 0, teta, scale = 0; /*anisotropy parameters - added by JH 2002 */
67  double xxr, yyr;
68 
69  if (params->theta) {
70  teta = params->theta * (M_PI / 180); /* deg to rad */
71  rsin = sin(teta);
72  rcos = cos(teta);
73  }
74  if (params->scalex)
75  scale = params->scalex;
76 
77 
78  n1 = n_points + 1;
79 
80  if (!A) {
81  if (!
82  (A =
83  G_alloc_vector((params->KMAX2 + 2) * (params->KMAX2 + 2) + 1))) {
84  fprintf(stderr, "Cannot allocate memory for A\n");
85  return -1;
86  }
87  }
88 
89  /*
90  C GENERATION OF MATRIX
91  C FIRST COLUMN
92  */
93  A[1] = 0.;
94  for (k = 1; k <= n_points; k++) {
95  i1 = k + 1;
96  A[i1] = 1.;
97  }
98  /*
99  C OTHER COLUMNS
100  */
101  RO = -params->rsm;
102  /* fprintf (stderr, "sm[%d] = %f, ro=%f\n", 1, points[1].smooth, RO); */
103  for (k = 1; k <= n_points; k++) {
104  k1 = k * n1 + 1;
105  k2 = k + 1;
106  i1 = k1 + k;
107  if (params->rsm < 0.) { /*indicates variable smoothing */
108  A[i1] = -points[k - 1].sm; /* added by Mitasova nov. 96 */
109  /* G_debug(5, "sm[%d]=%f, a=%f", k, points[k-1].sm, A[i1]); */
110  }
111  else {
112  A[i1] = RO; /* constant smoothing */
113  }
114  /* if (i1 == 100) fprintf (stderr,i "A[%d] = %f\n", i1, A[i1]); */
115 
116  /* A[i1] = RO; */
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;
120 
121  if ((params->theta) && (params->scalex)) {
122  /* re run anisotropy */
123  xxr = xx * rcos + yy * rsin;
124  yyr = yy * rcos - xx * rsin;
125  xx = xxr;
126  yy = yyr;
127  r = scale * xx * xx + yy * yy;
128  rfsta2 = fstar2 * (scale * xx * xx + yy * yy);
129  }
130  else {
131  r = xx * xx + yy * yy;
132  rfsta2 = fstar2 * (xx * xx + yy * yy);
133  }
134 
135  if (rfsta2 == 0.) {
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);
140  return -1;
141  }
142  i1 = k1 + l;
143  A[i1] = params->interp(r, params->fi);
144  }
145  }
146 
147  /* C SYMMETRISATION */
148  amaxa = 1.;
149  for (k = 1; k <= n1; k++) {
150  k1 = (k - 1) * n1;
151  k2 = k + 1;
152  for (l = k2; l <= n1; l++) {
153  m = (l - 1) * n1 + k;
154  A[m] = A[k1 + l];
155  amaxa = amax1(A[m], amaxa);
156  }
157  }
158  m = 0;
159  for (i = 0; i <= n_points; i++) {
160  for (j = 0; j <= n_points; j++) {
161  m++;
162  matrix[i][j] = A[m];
163  }
164  }
165 
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) {
168  /* find the inverse of the matrix */
169  fprintf(stderr, "G_ludcmp() failed! n=%d d=%.2f\n", n_points, d);
170  return -1;
171  }
172 
173  /* G_free_vector(A); */
174  return 1;
175 }
interp_fn * interp
Definition: interpf.h:101
double rsm
Definition: interpf.h:83
double scalex
Definition: interpf.h:90
double theta
Definition: interpf.h:89
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.
Definition: matrix.c:53
double amax1(double, double)
Definition: minmax.c:54
double l
#define NULL
Definition: ccmath.h:32
int G_ludcmp(double **a, int n, int *indx, double *d)
LU decomposition.
Definition: lu.c:18
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition: debug.c:65
double x
Definition: dataquad.h:42
double sm
Definition: dataquad.h:45
double fi
Definition: interpf.h:80
double * G_alloc_vector(size_t n)
Vector matrix memory allocation.
Definition: dalloc.c:41
double r
double y
Definition: dataquad.h:43