GRASS Programmer's Manual  6.4.2(2012)
secpar2d.c
Go to the documentation of this file.
00001 
00002 /*-
00003  * Written by H. Mitasova, L. Mitas, I. Kosinovsky, D. Gerdes Fall 1994
00004  * University of Illinois
00005  * US Army Construction Engineering Research Lab  
00006  * Copyright 1994, H. Mitasova (University of Illinois),
00007  * L. Mitas (University of Illinois),
00008  * I. Kosinovsky, (USA-CERL), and D.Gerdes (USA-CERL)   
00009  *
00010  * modified by McCauley in August 1995
00011  * modified by Mitasova in August 1995  
00012  *
00013  */
00014 
00015 #include <stdio.h>
00016 #include <math.h>
00017 #include <unistd.h>
00018 #include <grass/gis.h>
00019 #include <grass/bitmap.h>
00020 #include <grass/interpf.h>
00021 
00022 int IL_secpar_loop_2d(struct interp_params *params, int ngstc,  /* starting column */
00023                       int nszc, /* ending column */
00024                       int k,    /* current row */
00025                       struct BM *bitmask, double *gmin, double *gmax, double *c1min, double *c1max, double *c2min, double *c2max,       /* min,max interp.
00026                                                                                                                                          * values */
00027                       int cond1, int cond2      /* determine if particular values need to
00028                                                  * be computed */
00029     )
00030 
00031 /*
00032  * Computes slope, aspect and curvatures (depending on cond1, cond2) for
00033  * derivative arrays adx,...,adxy between columns ngstc and nszc.
00034  */
00035 {
00036     double dnorm1, ro,          /* rad to deg conv */
00037       dx2 = 0, dy2 = 0, grad2 = 0,      /* gradient squared */
00038         slp = 0, grad,          /* gradient */
00039         oor = 0,                /* aspect  (orientation) */
00040         curn = 0,               /* profile curvature */
00041         curh = 0,               /* tangential curvature */
00042         curm = 0,               /* mean curvature */
00043         temp,                   /* temp  variable */
00044         dxy2;                   /* temp variable   square of part diriv. */
00045 
00046     double gradmin;
00047     int i, got, bmask = 1;
00048     static int first_time_g = 1;
00049 
00050     ro = 57.295779;
00051     gradmin = 0.001;
00052 
00053 
00054     for (i = ngstc; i <= nszc; i++) {
00055         if (bitmask != NULL) {
00056             bmask = BM_get(bitmask, i, k);
00057         }
00058         got = 0;
00059         if (bmask == 1) {
00060             while ((got == 0) && (cond1)) {
00061                 dx2 = (double)(params->adx[i] * params->adx[i]);
00062                 dy2 = (double)(params->ady[i] * params->ady[i]);
00063                 grad2 = dx2 + dy2;
00064                 grad = sqrt(grad2);
00065                 /* slope in %        slp = 100. * grad; */
00066                 /* slope in degrees */
00067                 slp = ro * atan(grad);
00068                 if (grad <= gradmin) {
00069                     oor = 0.;
00070                     got = 3;
00071                     if (cond2) {
00072                         curn = 0.;
00073                         curh = 0.;
00074                         got = 3;
00075                         break;
00076                     }
00077                 }
00078                 if (got == 3)
00079                     break;
00080 
00081         /***********aspect from r.slope.aspect, with adx, ady computed
00082                     from interpol. function RST **************************/
00083 
00084                 if (params->adx[i] == 0.) {
00085                     if (params->ady[i] > 0.)
00086                         oor = 90;
00087                     else
00088                         oor = 270;
00089                 }
00090                 else {
00091                     oor = ro * atan2(params->ady[i], params->adx[i]);
00092                     if (oor <= 0.)
00093                         oor = 360. + oor;
00094                 }
00095 
00096                 got = 1;
00097             }                   /* while */
00098             if ((got != 3) && (cond2)) {
00099 
00100                 dnorm1 = sqrt(grad2 + 1.);
00101                 dxy2 =
00102                     2. * (double)(params->adxy[i] * params->adx[i] *
00103                                   params->ady[i]);
00104 
00105 
00106                 curn =
00107                     (double)(params->adxx[i] * dx2 + dxy2 +
00108                              params->adyy[i] * dy2) / (grad2 * dnorm1 *
00109                                                        dnorm1 * dnorm1);
00110 
00111                 curh =
00112                     (double)(params->adxx[i] * dy2 - dxy2 +
00113                              params->adyy[i] * dx2) / (grad2 * dnorm1);
00114 
00115                 temp = grad2 + 1.;
00116                 curm =
00117                     .5 * ((1. + dy2) * params->adxx[i] - dxy2 +
00118                           (1. + dx2) * params->adyy[i]) / (temp * dnorm1);
00119             }
00120             if (first_time_g) {
00121                 first_time_g = 0;
00122                 *gmin = *gmax = slp;
00123                 *c1min = *c1max = curn;
00124                 *c2min = *c2max = curh;
00125             }
00126             *gmin = amin1(*gmin, slp);
00127             *gmax = amax1(*gmax, slp);
00128             *c1min = amin1(*c1min, curn);
00129             *c1max = amax1(*c1max, curn);
00130             *c2min = amin1(*c2min, curh);
00131             *c2max = amax1(*c2max, curh);
00132             if (cond1) {
00133                 params->adx[i] = (FCELL) slp;
00134                 params->ady[i] = (FCELL) oor;
00135                 if (cond2) {
00136                     params->adxx[i] = (FCELL) curn;
00137                     params->adyy[i] = (FCELL) curh;
00138                     params->adxy[i] = (FCELL) curm;
00139                 }
00140             }
00141         }                       /* bmask == 1 */
00142     }
00143     return 1;
00144 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines