GRASS Programmer's Manual  6.5.svn(2014)-r66266
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
secpar2d.c
Go to the documentation of this file.
1 
2 /*-
3  * Written by H. Mitasova, L. Mitas, I. Kosinovsky, D. Gerdes Fall 1994
4  * University of Illinois
5  * US Army Construction Engineering Research Lab
6  * Copyright 1994, H. Mitasova (University of Illinois),
7  * L. Mitas (University of Illinois),
8  * I. Kosinovsky, (USA-CERL), and D.Gerdes (USA-CERL)
9  *
10  * modified by McCauley in August 1995
11  * modified by Mitasova in August 1995
12  *
13  */
14 
15 #include <stdio.h>
16 #include <math.h>
17 #include <unistd.h>
18 #include <grass/gis.h>
19 #include <grass/bitmap.h>
20 #include <grass/interpf.h>
21 
22 int IL_secpar_loop_2d(struct interp_params *params, int ngstc, /* starting column */
23  int nszc, /* ending column */
24  int k, /* current row */
25  struct BM *bitmask, double *gmin, double *gmax, double *c1min, double *c1max, double *c2min, double *c2max, /* min,max interp.
26  * values */
27  int cond1, int cond2 /* determine if particular values need to
28  * be computed */
29  )
30 
31 /*
32  * Computes slope, aspect and curvatures (depending on cond1, cond2) for
33  * derivative arrays adx,...,adxy between columns ngstc and nszc.
34  */
35 {
36  double dnorm1, ro, /* rad to deg conv */
37  dx2 = 0, dy2 = 0, grad2 = 0, /* gradient squared */
38  slp = 0, grad, /* gradient */
39  oor = 0, /* aspect (orientation) */
40  curn = 0, /* profile curvature */
41  curh = 0, /* tangential curvature */
42  curm = 0, /* mean curvature */
43  temp, /* temp variable */
44  dxy2; /* temp variable square of part diriv. */
45 
46  double gradmin;
47  int i, got, bmask = 1;
48  static int first_time_g = 1;
49 
50  ro = 57.295779;
51  gradmin = 0.001;
52 
53 
54  for (i = ngstc; i <= nszc; i++) {
55  if (bitmask != NULL) {
56  bmask = BM_get(bitmask, i, k);
57  }
58  got = 0;
59  if (bmask == 1) {
60  while ((got == 0) && (cond1)) {
61  dx2 = (double)(params->adx[i] * params->adx[i]);
62  dy2 = (double)(params->ady[i] * params->ady[i]);
63  grad2 = dx2 + dy2;
64  grad = sqrt(grad2);
65  /* slope in % slp = 100. * grad; */
66  /* slope in degrees */
67  slp = ro * atan(grad);
68  if (grad <= gradmin) {
69  oor = 0.;
70  got = 3;
71  if (cond2) {
72  curn = 0.;
73  curh = 0.;
74  got = 3;
75  break;
76  }
77  }
78  if (got == 3)
79  break;
80 
81  /***********aspect from r.slope.aspect, with adx, ady computed
82  from interpol. function RST **************************/
83 
84  if (params->adx[i] == 0.) {
85  if (params->ady[i] > 0.)
86  oor = 90;
87  else
88  oor = 270;
89  }
90  else {
91  oor = ro * atan2(params->ady[i], params->adx[i]);
92  if (oor <= 0.)
93  oor = 360. + oor;
94  }
95 
96  got = 1;
97  } /* while */
98  if ((got != 3) && (cond2)) {
99 
100  dnorm1 = sqrt(grad2 + 1.);
101  dxy2 =
102  2. * (double)(params->adxy[i] * params->adx[i] *
103  params->ady[i]);
104 
105 
106  curn =
107  (double)(params->adxx[i] * dx2 + dxy2 +
108  params->adyy[i] * dy2) / (grad2 * dnorm1 *
109  dnorm1 * dnorm1);
110 
111  curh =
112  (double)(params->adxx[i] * dy2 - dxy2 +
113  params->adyy[i] * dx2) / (grad2 * dnorm1);
114 
115  temp = grad2 + 1.;
116  curm =
117  .5 * ((1. + dy2) * params->adxx[i] - dxy2 +
118  (1. + dx2) * params->adyy[i]) / (temp * dnorm1);
119  }
120  if (first_time_g) {
121  first_time_g = 0;
122  *gmin = *gmax = slp;
123  *c1min = *c1max = curn;
124  *c2min = *c2max = curh;
125  }
126  *gmin = amin1(*gmin, slp);
127  *gmax = amax1(*gmax, slp);
128  *c1min = amin1(*c1min, curn);
129  *c1max = amax1(*c1max, curn);
130  *c2min = amin1(*c2min, curh);
131  *c2max = amax1(*c2max, curh);
132  if (cond1) {
133  params->adx[i] = (FCELL) slp;
134  params->ady[i] = (FCELL) oor;
135  if (cond2) {
136  params->adxx[i] = (FCELL) curn;
137  params->adyy[i] = (FCELL) curh;
138  params->adxy[i] = (FCELL) curm;
139  }
140  }
141  } /* bmask == 1 */
142  }
143  return 1;
144 }
int BM_get(struct BM *map, int x, int y)
Gets &#39;val&#39; from the bitmap.
Definition: bitmap.c:220
double amax1(double, double)
Definition: minmax.c:54
DCELL * adxx
Definition: interpf.h:48
DCELL * adx
Definition: interpf.h:48
int IL_secpar_loop_2d(struct interp_params *, int, int, int, struct BM *, double *, double *, double *, double *, double *, double *, int, int)
Definition: secpar2d.c:22
double amin1(double, double)
Definition: minmax.c:67
DCELL * ady
Definition: interpf.h:48
return NULL
Definition: dbfopen.c:1394
DCELL * adxy
Definition: interpf.h:48
DCELL * adyy
Definition: interpf.h:48