GRASS Programmer's Manual  6.4.2(2012)
matrix.c
Go to the documentation of this file.
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 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines