Introduction
Introduction Statistics Contact Development Disclaimer Help
twerner.c - werner - cellular automata simulation of wind-driven sand transport
git clone git://src.adamsgaard.dk/werner
Log
Files
Refs
LICENSE
---
twerner.c (9111B)
---
1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <gsl/gsl_matrix.h>
4
5 // Return a random number in [0;1[
6 double rnd(void)
7 {
8 return (double)rand()/RAND_MAX;
9 }
10
11 // Initialize matrix M with random ints from low to high
12 void init_rnd_matrix(gsl_matrix* M, double low, double high)
13 {
14 int row, col;
15 for (row = 0; row < M->size1; row++)
16 for (col = 0; col < M->size2; col++)
17 gsl_matrix_set(M, row, col,
18 (int)(rnd() * (high-low) + low));
19 }
20
21 // Write matrix to stdout
22 void print_matrix(gsl_matrix* M)
23 {
24 int row, col;
25 double val;
26 for (row = 0; row < M->size1; row++) {
27 for (col = 0; col < M->size2; col++) {
28 val = gsl_matrix_get(M, row, col);
29 printf("%f\t", val);
30 }
31 puts(""); // newline
32 }
33 }
34
35 // Add a slab of sand from the target cell
36 void deposit(gsl_matrix* Z, int row, int col)
37 {
38 gsl_matrix_set(Z, row, col, gsl_matrix_get(Z, row, col) + 1.0);
39 }
40
41 // Remove a slab of sand from the target cell
42 void erode(gsl_matrix* Z, int row, int col)
43 {
44 gsl_matrix_set(Z, row, col, gsl_matrix_get(Z, row, col) - 1.0);
45 }
46
47 // Find out if cell is in shade (1) or not (0)
48 int inshade(
49 gsl_matrix* Z, // matrix with sand slabs
50 int row, // row idx of interest
51 int col) // col idx of interest
52 {
53 int js;
54 int i = 1;
55 int shade = 0;
56 int nx = Z->size2;
57
58 while ((i <= nx/4) && (shade == 0)) {
59 js = col - i;
60 if (js < 0) // periodic boundary
61 js += nx;
62 if (gsl_matrix_get(Z, row, js) >= gsl_matrix_get(Z, row, col) + …
63 shade = 1;
64 i++;
65 }
66 return shade;
67 }
68
69 void find_max_slope_neighbor_ero(
70 gsl_matrix* Z, // sand slab values
71 int row, // target row
72 int col, // target col
73 int* row_max, // max slope neighbor row
74 int* col_max, // max slope neighbor col
75 double* dh_max) // max slope
76 {
77
78 int n; // 1d neighbor idx
79 int i, j; // 2d neighbor idx
80 double dh; // slope
81
82 // relative neighbor coordinates and distances
83 int drow[] = {1, 1, 1, 0, -1, -1, -1, 0};
84 int dcol[] = {-1, 0, 1, 1, 1, 0, -1, -1};
85 double dx[] = {1.4142, 1.0, 1.4142, 1.0, 1.4142, 1.0, 1.4142, 1.0};
86
87 // matrix dimensions
88 int nrows = Z->size1;
89 int ncols = Z->size2;
90
91 // check the 8 surrounding neighbor cells
92 for (n = 0; n<8; n++) {
93
94 // absolute neighbor cell idx
95 i = row + drow[n];
96 j = col + dcol[n];
97
98 // correct for periodic boundaries
99 if (i < 0)
100 i += nrows;
101 if (i >= nrows)
102 i -= nrows;
103 if (j < 0)
104 j += ncols;
105 if (j >= ncols)
106 j -= ncols;
107
108 // find slope
109 dh = (gsl_matrix_get(Z, i, j) - gsl_matrix_get(Z, row, col)) / d…
110
111 // save position if it is the highest slope so far
112 if (dh > *dh_max) {
113 *dh_max = dh;
114 *row_max = i;
115 *col_max = j;
116 }
117 }
118 }
119
120 void find_max_slope_neighbor_depo(
121 gsl_matrix* Z, // sand slab values
122 int row, // target row
123 int col, // target col
124 int* row_max, // max slope neighbor row
125 int* col_max, // max slope neighbor col
126 double* dh_max) // max slope
127 {
128
129 int n; // 1d neighbor idx
130 int i, j; // 2d neighbor idx
131 double dh; // slope
132
133 // relative neighbor coordinates and distances
134 int drow[] = {1, 1, 1, 0, -1, -1, -1, 0};
135 int dcol[] = {-1, 0, 1, 1, 1, 0, -1, -1};
136 double dx[] = {1.4142, 1.0, 1.4142, 1.0, 1.4142, 1.0, 1.4142, 1.0};
137
138 // matrix dimensions
139 int nrows = Z->size1;
140 int ncols = Z->size2;
141
142 // check the 8 surrounding neighbor cells
143 for (n = 0; n<8; n++) {
144
145 // absolute neighbor cell idx
146 i = row + drow[n];
147 j = col + dcol[n];
148
149 // correct for periodic boundaries
150 if (i < 0)
151 i += nrows;
152 if (i >= nrows)
153 i -= nrows;
154 if (j < 0)
155 j += ncols;
156 if (j >= ncols)
157 j -= ncols;
158
159 // find slope
160 dh = (gsl_matrix_get(Z, row, col) - gsl_matrix_get(Z, i, j)) / d…
161
162 // save position if it is the highest slope so far
163 if (dh > *dh_max) {
164 *dh_max = dh;
165 *row_max = i;
166 *col_max = j;
167 }
168 }
169 }
170
171 // Check and perform avalanche into cell if slope exceeds limit
172 void avalanche_erosion(
173 gsl_matrix* Z, // sand slab values
174 int row, // target row
175 int col) // target col
176 {
177 // find the neighbor with the max. height difference and the slope v…
178 double dh_max = 0.0; // max slope
179 int row_max, col_max;
180 find_max_slope_neighbor_ero(Z, row, col, &row_max, &col_max, &dh_max…
181
182 // Perform avalanche along max slope neighbor
183 double threshold = 2.0; // avalanche threshold
184 if (dh_max > threshold) {
185 deposit(Z, row, col); // put sand in target cell
186 erode(Z, row_max, col_max); // remove sand from max slope neighb…
187
188 // Recursion
189 avalanche_erosion(Z, row, col);
190 }
191 }
192
193 // Check and perform avalanche away from cell if slope exceeds limit
194 void avalanche_deposition(
195 gsl_matrix* Z, // sand slab values
196 int row, // target row
197 int col) // target col
198 {
199 // find the neighbor with the max. height difference and the slope v…
200 double dh_max = 0.0; // max slope
201 int row_max, col_max;
202 find_max_slope_neighbor_depo(Z, row, col, &row_max, &col_max, &dh_ma…
203
204 // Perform avalanche along max slope neighbor
205 double threshold = 2.0; // avalanche threshold
206 if (dh_max > threshold) {
207 erode(Z, row, col);
208 deposit(Z, row_max, col_max);
209
210 // Recursion
211 avalanche_deposition(Z, row, col);
212 }
213 }
214
215 // Move a sand unit and deposit it where it fits
216 void move_and_deposit(
217 gsl_matrix* Z, // matrix with sand slabs
218 double p_ns, // likelihood of deposition in no-sand cell
219 double p_s, // likelihood of deposition in sand cell
220 double d_l, // transport distance in no. of cells
221 int* row, // current cell row
222 int* col, // current cell col
223 int* deposited) // deposited? 1=yes, 0=no
224 {
225 // move sand slab in wind direction
226 *col += d_l;
227 int ncols = Z->size1;
228 if (*col >= ncols)
229 *col -= ncols;
230
231 // is the target in the shade?
232 // 1=yes, 0=no
233 int shade = 0;
234 shade = inshade(Z, *row, *col);
235 if (shade > 0) {
236 deposit(Z, *row, *col);
237 avalanche_deposition(Z, *row, *col); // check for avalanches
238 *deposited = 1; // sand deposited, stop moving it
239 }
240
241 // if not in the shade, check against deposition probabilities
242 else {
243
244 // sand in cell
245 if (gsl_matrix_get(Z, *row, *col) > 0.) {
246 if (rnd() <= p_s) { // deposit in cell with sand
247 deposit(Z, *row, *col);
248 avalanche_deposition(Z, *row, *col); // check for avalan…
249 *deposited = 1; // sand deposited, stop moving it
250 }
251 } else { // no sand in cell
252 if (rnd() <= p_ns) { // deposit in cell without sand
253 deposit(Z, *row, *col);
254 avalanche_deposition(Z, *row, *col); // check for avalan…
255 *deposited = 1; // sand deposited, stop moving it
256 }
257 }
258 }
259 }
260
261 // Perform a single Werner iteration
262 void werner_iter(
263 gsl_matrix* Z, // matrix with sand slabs
264 int d_l, // wind transport distance
265 double p_ns, // likelihood of deposition in sand-free cell
266 double p_s) // likelihood of deposition in sand-covered cell
267 {
268 // Evaluate N=nx*ny cells in random order
269 int row, col, deposited;
270 int nrows = Z->size1; // row major
271 int ncols = Z->size2;
272 long int n;
273 long int N = nrows*ncols;
274 double cellval;
275
276
277 // Perform cycle
278 for (n=0; n<N; n++) {
279
280 // random cell idx
281 row = rand()%nrows;
282 col = rand()%ncols;
283
284 // If the cell has sand in it (val>0), and is not in the shadow …
285 cellval = gsl_matrix_get(Z, row, col);
286 if ((cellval > 0.) && (inshade(Z, row, col) == 0)) {
287
288 // erode
289 erode(Z, row, col);
290
291 // check for avalanche
292 avalanche_erosion(Z, row, col);
293
294 // move eroded sand and deposit where it fits
295 deposited = 0;
296 while (deposited == 0)
297 move_and_deposit(Z, p_ns, p_s, d_l, &row, &col, &deposit…
298
299 }
300 }
301 }
302
303 // Loop system through time
304 void werner_loop(
305 gsl_matrix* Z, // matrix with sand slabs
306 long int N_c, // number of iterations
307 int d_l, // wind transport distance
308 double p_ns, // likelihood of deposition in sand-free cell
309 double p_s) // likelihood of deposition in sand-covered cell
310 {
311 long int t;
312 for (t = 0; t<N_c; t++)
313 werner_iter(Z, d_l, p_ns, p_s);
314 }
You are viewing proxied material from mx1.adamsgaard.dk. The copyright of proxied material belongs to its original authors. Any comments or complaints in relation to proxied material should be directed to the original authors of the content concerned. Please see the disclaimer for more details.