GRASS Programmer's Manual  6.5.svn(2014)-r66266
matrix.c
Go to the documentation of this file.
1 /*
2  *
3  * Original program and various modifications:
4  * Lubos Mitas
5  *
6  * GRASS4.1 version of the program and GRASS4.2 modifications:
7  * H. Mitasova,
8  * I. Kosinovsky, D. Gerdes
9  * D. McCauley
10  *
12  * L. Mitas ,
13  * H. Mitasova ,
14  * I. Kosinovsky,
15  * D.Gerdes
16  * D. McCauley
17  *
18  * modified by McCauley in August 1995
19  * modified by Mitasova in August 1995, Nov. 1996
20  *
21  */
22
23 #include <stdio.h>
24 #include <math.h>
25 #include <unistd.h>
26 #include <grass/gis.h>
27 #include <grass/interpf.h>
28 #include <grass/gmath.h>
29
30 int IL_matrix_create(struct interp_params *params, struct triple *points, /* points for interpolation */
31  int n_points, /* number of points */
32  double **matrix, /* matrix */
33  int *indx)
34 /*
35  Creates system of linear equations represented by matrix using given points
36  and interpolating function interp()
37  */
38 {
39  double xx, yy;
40  double rfsta2, r;
41  double d;
42  int n1, k1, k2, k, i1, l, m, i, j;
43  double fstar2 = params->fi * params->fi / 4.;
44  static double *A = NULL;
45  double RO, amaxa;
46  double rsin = 0, rcos = 0, teta, scale = 0; /*anisotropy parameters - added by JH 2002 */
47  double xxr, yyr;
48
49  if (params->theta) {
50  teta = params->theta / 57.295779; /* deg to rad */
51  rsin = sin(teta);
52  rcos = cos(teta);
53  }
54  if (params->scalex)
55  scale = params->scalex;
56
57
58  n1 = n_points + 1;
59
60  if (!A) {
61  if (!
62  (A =
63  G_alloc_vector((params->KMAX2 + 2) * (params->KMAX2 + 2) + 1))) {
64  fprintf(stderr, "Cannot allocate memory for A\n");
65  return -1;
66  }
67  }
68
69  /*
70  C
71  C GENERATION OF MATRIX
72  C
73  C FIRST COLUMN
74  C
75  */
76  A[1] = 0.;
77  for (k = 1; k <= n_points; k++) {
78  i1 = k + 1;
79  A[i1] = 1.;
80  }
81  /*
82  C
83  C OTHER COLUMNS
84  C
85  */
86  RO = -params->rsm;
87  /* fprintf (stderr,"sm[%d]=%f,ro=%f\n",1,points[1].smooth,RO); */
88  for (k = 1; k <= n_points; k++) {
89  k1 = k * n1 + 1;
90  k2 = k + 1;
91  i1 = k1 + k;
92  if (params->rsm < 0.) { /*indicates variable smoothing */
93  A[i1] = -points[k - 1].sm; /* added by Mitasova nov. 96 */
94  /* fprintf (stderr,"sm[%d]=%f,a=%f\n",k,points[k-1].sm,A[i1]); */
95  }
96  else {
97  A[i1] = RO; /* constant smoothing */
98  }
99  /* if (i1 == 100) fprintf (stderr,"A[%d]=%f\n",i1,A[i1]); */
100
101  /* A[i1] = RO; */
102  for (l = k2; l <= n_points; l++) {
103  xx = points[k - 1].x - points[l - 1].x;
104  yy = points[k - 1].y - points[l - 1].y;
105
106  if ((params->theta) && (params->scalex)) {
107  /* re run anisotropy */
108  xxr = xx * rcos + yy * rsin;
109  yyr = yy * rcos - xx * rsin;
110  xx = xxr;
111  yy = yyr;
112  r = scale * xx * xx + yy * yy;
113  rfsta2 = fstar2 * (scale * xx * xx + yy * yy);
114  }
115  else {
116  r = xx * xx + yy * yy;
117  rfsta2 = fstar2 * (xx * xx + yy * yy);
118  }
119
120  if (rfsta2 == 0.) {
121  fprintf(stderr, "ident. points in segm. \n");
122  fprintf(stderr, "x[%d]=%f,x[%d]=%f,y[%d]=%f,y[%d]=%f\n",
123  k - 1, points[k - 1].x, l - 1, points[l - 1].x, k - 1,
124  points[k - 1].y, l - 1, points[l - 1].y);
125  return -1;
126  }
127  i1 = k1 + l;
128  A[i1] = params->interp(r, params->fi);
129  }
130  }
131  /*
132  C
133  C SYMMETRISATION
134  C
135  */
136  amaxa = 1.;
137  for (k = 1; k <= n1; k++) {
138  k1 = (k - 1) * n1;
139  k2 = k + 1;
140  for (l = k2; l <= n1; l++) {
141  m = (l - 1) * n1 + k;
142  A[m] = A[k1 + l];
143  amaxa = amax1(A[m], amaxa);
144  }
145  }
146  m = 0;
147  for (i = 0; i <= n_points; i++) {
148  for (j = 0; j <= n_points; j++) {
149  m++;
150  matrix[i][j] = A[m];
151  }
152  }
153
154  if (G_ludcmp(matrix, n_points + 1, indx, &d) <= 0) { /* find the inverse of the mat
155  rix */
156  fprintf(stderr, "G_ludcmp() failed! n=%d\n", n_points);
157  return -1;
158  }
159  /*
160  G_free_vector(A);
161  */
162  return 1;
163 }
int l
double rsm
Definition: interpf.h:52
double scalex
Definition: interpf.h:58
float r
Definition: named_colr.c:8
double theta
Definition: interpf.h:57
int IL_matrix_create(struct interp_params *, struct triple *, int, double **, int *)
Definition: matrix.c:30
double amax1(double, double)
Definition: minmax.c:54
double(* interp)()
Definition: interpf.h:68
int y
Definition: plot.c:34
int G_ludcmp(double **a, int n, int *indx, double *d)
Definition: lu.c:9
double x