GRASS Programmer's Manual
6.4.2(2012)
|
00001 /* 00002 * 00003 * Original program and various modifications: 00004 * Lubos Mitas 00005 * 00006 * GRASS4.1 version of the program and GRASS4.2 modifications: 00007 * H. Mitasova, 00008 * I. Kosinovsky, D. Gerdes 00009 * D. McCauley 00010 * 00011 * Copyright 1993, 1995: 00012 * L. Mitas , 00013 * H. Mitasova , 00014 * I. Kosinovsky, 00015 * D.Gerdes 00016 * D. McCauley 00017 * 00018 * modified by McCauley in August 1995 00019 * modified by Mitasova in August 1995, Nov. 1996 00020 * 00021 */ 00022 00023 #include <stdio.h> 00024 #include <math.h> 00025 #include <unistd.h> 00026 #include <grass/gis.h> 00027 #include <grass/interpf.h> 00028 #include <grass/gmath.h> 00029 00030 int IL_matrix_create(struct interp_params *params, struct triple *points, /* points for interpolation */ 00031 int n_points, /* number of points */ 00032 double **matrix, /* matrix */ 00033 int *indx) 00034 /* 00035 Creates system of linear equations represented by matrix using given points 00036 and interpolating function interp() 00037 */ 00038 { 00039 double xx, yy; 00040 double rfsta2, r; 00041 double d; 00042 int n1, k1, k2, k, i1, l, m, i, j; 00043 double fstar2 = params->fi * params->fi / 4.; 00044 static double *A = NULL; 00045 double RO, amaxa; 00046 double rsin = 0, rcos = 0, teta, scale = 0; /*anisotropy parameters - added by JH 2002 */ 00047 double xxr, yyr; 00048 00049 if (params->theta) { 00050 teta = params->theta / 57.295779; /* deg to rad */ 00051 rsin = sin(teta); 00052 rcos = cos(teta); 00053 } 00054 if (params->scalex) 00055 scale = params->scalex; 00056 00057 00058 n1 = n_points + 1; 00059 00060 if (!A) { 00061 if (! 00062 (A = 00063 G_alloc_vector((params->KMAX2 + 2) * (params->KMAX2 + 2) + 1))) { 00064 fprintf(stderr, "Cannot allocate memory for A\n"); 00065 return -1; 00066 } 00067 } 00068 00069 /* 00070 C 00071 C GENERATION OF MATRIX 00072 C 00073 C FIRST COLUMN 00074 C 00075 */ 00076 A[1] = 0.; 00077 for (k = 1; k <= n_points; k++) { 00078 i1 = k + 1; 00079 A[i1] = 1.; 00080 } 00081 /* 00082 C 00083 C OTHER COLUMNS 00084 C 00085 */ 00086 RO = -params->rsm; 00087 /* fprintf (stderr,"sm[%d]=%f,ro=%f\n",1,points[1].smooth,RO); */ 00088 for (k = 1; k <= n_points; k++) { 00089 k1 = k * n1 + 1; 00090 k2 = k + 1; 00091 i1 = k1 + k; 00092 if (params->rsm < 0.) { /*indicates variable smoothing */ 00093 A[i1] = -points[k - 1].sm; /* added by Mitasova nov. 96 */ 00094 /* fprintf (stderr,"sm[%d]=%f,a=%f\n",k,points[k-1].sm,A[i1]); */ 00095 } 00096 else { 00097 A[i1] = RO; /* constant smoothing */ 00098 } 00099 /* if (i1 == 100) fprintf (stderr,"A[%d]=%f\n",i1,A[i1]); */ 00100 00101 /* A[i1] = RO; */ 00102 for (l = k2; l <= n_points; l++) { 00103 xx = points[k - 1].x - points[l - 1].x; 00104 yy = points[k - 1].y - points[l - 1].y; 00105 00106 if ((params->theta) && (params->scalex)) { 00107 /* re run anisotropy */ 00108 xxr = xx * rcos + yy * rsin; 00109 yyr = yy * rcos - xx * rsin; 00110 xx = xxr; 00111 yy = yyr; 00112 r = scale * xx * xx + yy * yy; 00113 rfsta2 = fstar2 * (scale * xx * xx + yy * yy); 00114 } 00115 else { 00116 r = xx * xx + yy * yy; 00117 rfsta2 = fstar2 * (xx * xx + yy * yy); 00118 } 00119 00120 if (rfsta2 == 0.) { 00121 fprintf(stderr, "ident. points in segm. \n"); 00122 fprintf(stderr, "x[%d]=%f,x[%d]=%f,y[%d]=%f,y[%d]=%f\n", 00123 k - 1, points[k - 1].x, l - 1, points[l - 1].x, k - 1, 00124 points[k - 1].y, l - 1, points[l - 1].y); 00125 return -1; 00126 } 00127 i1 = k1 + l; 00128 A[i1] = params->interp(r, params->fi); 00129 } 00130 } 00131 /* 00132 C 00133 C SYMMETRISATION 00134 C 00135 */ 00136 amaxa = 1.; 00137 for (k = 1; k <= n1; k++) { 00138 k1 = (k - 1) * n1; 00139 k2 = k + 1; 00140 for (l = k2; l <= n1; l++) { 00141 m = (l - 1) * n1 + k; 00142 A[m] = A[k1 + l]; 00143 amaxa = amax1(A[m], amaxa); 00144 } 00145 } 00146 m = 0; 00147 for (i = 0; i <= n_points; i++) { 00148 for (j = 0; j <= n_points; j++) { 00149 m++; 00150 matrix[i][j] = A[m]; 00151 } 00152 } 00153 00154 if (G_ludcmp(matrix, n_points + 1, indx, &d) <= 0) { /* find the inverse of the mat 00155 rix */ 00156 fprintf(stderr, "G_ludcmp() failed! n=%d\n", n_points); 00157 return -1; 00158 } 00159 /* 00160 G_free_vector(A); 00161 */ 00162 return 1; 00163 }