GRASS GIS 8 Programmer's Manual  8.2.2dev(2023)-3d2c704037
secpar2d.c
Go to the documentation of this file.
1 
2 /*!
3  * \file secpar2d.c
4  *
5  * \author H. Mitasova, L. Mitas, I. Kosinovsky, D. Gerdes Fall 1994 (original authors)
6  * \author modified by McCauley in August 1995
7  * \author modified by Mitasova in August 1995
8  * \author H. Mitasova (University of Illinois)
9  * \author L. Mitas (University of Illinois)
10  * \author I. Kosinovsky, (USA-CERL)
11  * \author D.Gerdes (USA-CERL)
12  *
13  * \copyright
14  * (C) 1994-1995 by Helena Mitasova and the GRASS Development Team
15  *
16  * \copyright
17  * This program is free software under the
18  * GNU General Public License (>=v2).
19  * Read the file COPYING that comes with GRASS
20  * for details.
21  *
22  */
23 
24 
25 #include <stdio.h>
26 #include <math.h>
27 #include <unistd.h>
28 #include <grass/gis.h>
29 #include <grass/bitmap.h>
30 #include <grass/interpf.h>
31 
32 
33 /*!
34  * Compute slope aspect and curvatures
35  *
36  * Computes slope, aspect and curvatures (depending on cond1, cond2) for
37  * derivative arrays adx,...,adxy between columns ngstc and nszc.
38  */
39 int IL_secpar_loop_2d(struct interp_params *params,
40  int ngstc, /*!< starting column */
41  int nszc, /*!< ending column */
42  int k, /*!< current row */
43  struct BM *bitmask,
44  double *gmin, double *gmax,
45  double *c1min, double *c1max,
46  double *c2min, double *c2max, /*!< min,max interp. values */
47  int cond1,
48  int cond2 /*!< determine if particular values need to be computed */
49  )
50 {
51  double dnorm1, ro, /* rad to deg conv */
52  dx2 = 0, dy2 = 0, grad2 = 0, /* gradient squared */
53  slp = 0, grad, /* gradient */
54  oor = 0, /* aspect (orientation) */
55  curn = 0, /* profile curvature */
56  curh = 0, /* tangential curvature */
57  curm = 0, /* mean curvature */
58  temp, /* temp variable */
59  dxy2; /* temp variable square of part diriv. */
60 
61  double gradmin;
62  int i, got, bmask = 1;
63  static int first_time_g = 1;
64 
65  ro = M_R2D;
66  gradmin = 0.001;
67 
68 
69  for (i = ngstc; i <= nszc; i++) {
70  if (bitmask != NULL) {
71  bmask = BM_get(bitmask, i, k);
72  }
73  got = 0;
74  if (bmask == 1) {
75  while ((got == 0) && (cond1)) {
76  dx2 = (double)(params->adx[i] * params->adx[i]);
77  dy2 = (double)(params->ady[i] * params->ady[i]);
78  grad2 = dx2 + dy2;
79  grad = sqrt(grad2);
80  /* slope in % slp = 100. * grad; */
81  /* slope in degrees */
82  slp = ro * atan(grad);
83  if (grad <= gradmin) {
84  oor = 0.;
85  got = 3;
86  if (cond2) {
87  curn = 0.;
88  curh = 0.;
89  got = 3;
90  break;
91  }
92  }
93  if (got == 3)
94  break;
95 
96  /***********aspect from r.slope.aspect, with adx, ady computed
97  from interpol. function RST **************************/
98 
99  if (params->adx[i] == 0.) {
100  if (params->ady[i] > 0.)
101  oor = 90;
102  else
103  oor = 270;
104  }
105  else {
106  oor = ro * atan2(params->ady[i], params->adx[i]);
107  if (oor <= 0.)
108  oor = 360. + oor;
109  }
110 
111  got = 1;
112  } /* while */
113  if ((got != 3) && (cond2)) {
114 
115  dnorm1 = sqrt(grad2 + 1.);
116  dxy2 =
117  2. * (double)(params->adxy[i] * params->adx[i] *
118  params->ady[i]);
119 
120 
121  curn =
122  (double)(params->adxx[i] * dx2 + dxy2 +
123  params->adyy[i] * dy2) / (grad2 * dnorm1 *
124  dnorm1 * dnorm1);
125 
126  curh =
127  (double)(params->adxx[i] * dy2 - dxy2 +
128  params->adyy[i] * dx2) / (grad2 * dnorm1);
129 
130  temp = grad2 + 1.;
131  curm =
132  .5 * ((1. + dy2) * params->adxx[i] - dxy2 +
133  (1. + dx2) * params->adyy[i]) / (temp * dnorm1);
134  }
135  if (first_time_g) {
136  first_time_g = 0;
137  *gmin = *gmax = slp;
138  *c1min = *c1max = curn;
139  *c2min = *c2max = curh;
140  }
141  *gmin = amin1(*gmin, slp);
142  *gmax = amax1(*gmax, slp);
143  *c1min = amin1(*c1min, curn);
144  *c1max = amax1(*c1max, curn);
145  *c2min = amin1(*c2min, curh);
146  *c2max = amax1(*c2max, curh);
147  if (cond1) {
148  params->adx[i] = (FCELL) slp;
149  params->ady[i] = (FCELL) oor;
150  if (cond2) {
151  params->adxx[i] = (FCELL) curn;
152  params->adyy[i] = (FCELL) curh;
153  params->adxy[i] = (FCELL) curm;
154  }
155  }
156  } /* bmask == 1 */
157  }
158  return 1;
159 }
int IL_secpar_loop_2d(struct interp_params *params, int ngstc, int nszc, int k, struct BM *bitmask, double *gmin, double *gmax, double *c1min, double *c1max, double *c2min, double *c2max, int cond1, int cond2)
Definition: secpar2d.c:39
Definition: bitmap.h:17
double amax1(double, double)
Definition: minmax.c:54
DCELL * adxx
Definition: interpf.h:78
#define M_R2D
Definition: gis.h:153
#define NULL
Definition: ccmath.h:32
DCELL * adx
Definition: interpf.h:78
int BM_get(struct BM *, int, int)
Gets &#39;val&#39; from the bitmap.
Definition: bitmap.c:223
float FCELL
Definition: gis.h:615
double amin1(double, double)
Definition: minmax.c:67
DCELL * ady
Definition: interpf.h:78
DCELL * adxy
Definition: interpf.h:78
DCELL * adyy
Definition: interpf.h:78