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