TYCHO  1.3.0
 All Data Structures Files Functions Variables Enumerations Enumerator
/home/kapf/tycho_docu/noise_generator.c
Go to the documentation of this file.
1 /*
2  * noise_generator.c
3  *
4  * Author: Wolfgang Kapferer
5  */
6 
7 #include <stdio.h>
8 #include <stdlib.h>
9 #include <math.h>
10 
11 #include "variables_global.h"
12 #include "prototypes.h"
13 
14 #ifdef _OPENMP
15 #include <omp.h>
16 #endif
17 
22 int noise_generator(int x, int y, int z) {
23 
24  double local_pressure_left, local_pressure_right;
25  int one_boundary;
26  int left, right;
27  int counter_mean_pressure;
28  int i, j, k;
29  double tmp;
30  double perfect_sinus_emitter;
31 
32  perfect_sinus_emitter = 1.0 / sqrt(2);
33 
34  mean_pressure = 0.0;
35  counter_mean_pressure = 0;
36  //==========================================================Here the mean pressure in the domain===============================================
37 
38 #ifdef _OPENMP
39 #pragma omp parallel default(none) \
40  private(j, i, k, local_pressure_left, \
41  local_pressure_right, one_boundary, \
42  left, right) \
43  shared(x, y, z, pre, dom, \
44  pressure_integrated, mean_pressure, counter_mean_pressure)
45  {
46 #endif
47 
48  for (i = 0; i < x; i++) {
49 #ifdef _OPENMP
50 #pragma omp for schedule(static)
51 #endif
52  for (j = 0; j < y; j++) {
53  for (k = 0; k < z; k++) {
54 
55  if (dom[i][j][k] == 0) {
56 
57 #ifdef _OPENMP
58 #pragma omp critical
59 #endif
60  mean_pressure += pre[i][j][k];
61 #ifdef _OPENMP
62 #pragma omp end critical
63 #endif
64 
65 #ifdef _OPENMP
66 #pragma omp critical
67 #endif
68  counter_mean_pressure++;
69 #ifdef _OPENMP
70 #pragma omp end critical
71 #endif
72 
73  }
74  }
75  }
76  }
77 
78 #ifdef _OPENMP
79  }
80 #endif
81 
82  mean_pressure = mean_pressure / counter_mean_pressure;
83  if (with_one_pulse == 0) pressure = mean_pressure + sound_pressure_level * sin(2 * M_PI * sound_frequency * time_sim);
84 
85  //just one pulse
86  if ((with_one_pulse == 1) && (counter == 0)) {
87  pressure = fabs(mean_pressure + sound_pressure_level * sin(0.5 * M_PI));
88  }
89 
90  if (with_one_pulse == 2) {
91  pressure = fabs(mean_pressure + sound_pressure_level * sin(0.5 * M_PI));
92  }
93 
94  if (dimension == 3) {
95 
96 #ifdef _OPENMP
97 #pragma omp parallel default(none) \
98  private(j, i, k) \
99  shared(x, y, z, pre_old, pre, \
100  soundemitter, dom, pressure, \
101  pressure_on_solid, with_one_pulse, counter)
102  {
103 #endif
104  for (i = 1; i < x - 1; i++) {
105 
106 #ifdef _OPENMP
107 #pragma omp for schedule(static)
108 #endif
109 
110  for (j = 1; j < y - 1; j++) {
111  for (k = 1; k < z - 1; k++) {
112 
113  if ((with_one_pulse == 1) && (counter == 0)) {
114  if (soundemitter[i - 1][j - 1][k - 1] != 1) pre[i - 1][j - 1][k - 1] = pressure;
115  if (soundemitter[i - 1][j ][k - 1] != 1) pre[i - 1][j ][k - 1] = pressure;
116  if (soundemitter[i - 1][j + 1][k - 1] != 1) pre[i - 1][j + 1][k - 1] = pressure;
117  if (soundemitter[i - 1][j - 1][k] != 1) pre[i - 1][j - 1][k] = pressure;
118  if (soundemitter[i - 1][j][k] != 1) pre[i - 1][j][k] = pressure;
119  if (soundemitter[i - 1][j + 1][k] != 1) pre[i - 1][j + 1][k] = pressure;
120  if (soundemitter[i - 1][j - 1][k + 1] != 1) pre[i - 1][j - 1][k + 1] = pressure;
121  if (soundemitter[i - 1][j ][k + 1] != 1) pre[i - 1][j][k + 1] = pressure;
122  if (soundemitter[i - 1][j + 1][k + 1] != 1) pre[i - 1][j + 1][k + 1] = pressure;
123 
124  if (soundemitter[i][j - 1][k - 1] != 1) pre[i][j - 1][k - 1] = pressure;
125  if (soundemitter[i][j][k - 1] != 1) pre[i][j][k - 1] = pressure;
126  if (soundemitter[i][j + 1][k - 1] != 1) pre[i][j + 1][k - 1] = pressure;
127  if (soundemitter[i][j - 1][k] != 1) pre[i][j - 1][k] = pressure;
128  if (soundemitter[i][j + 1][k] != 1) pre[i][j + 1][k - 1] = pressure;
129  if (soundemitter[i][j - 1][k + 1] != 1) pre[i][j - 1][k + 1] = pressure;
130  if (soundemitter[i][j][k + 1] != 1) pre[i][j][k + 1] = pressure;
131  if (soundemitter[i][j + 1][k + 1] != 1) pre[i][j + 1][k + 1] = pressure;
132 
133  if (soundemitter[i + 1][j - 1][k - 1] != 1) pre[i + 1][j - 1][k - 1] = pressure;
134  if (soundemitter[i + 1][j][k - 1] != 1) pre[i + 1][j][k - 1] = pressure;
135  if (soundemitter[i + 1][j + 1][k - 1] != 1) pre[i + 1][j + 1][k - 1] = pressure;
136  if (soundemitter[i + 1][j - 1][k] != 1) pre[i + 1][j - 1][k] = pressure;
137  if (soundemitter[i + 1][j][k] != 1) pre[i + 1][j][k] = pressure;
138  if (soundemitter[i + 1][j + 1][k] != 1) pre[i + 1][j + 1][k] = pressure;
139  if (soundemitter[i + 1][j - 1][k + 1] != 1) pre[i + 1][j - 1][k + 1] = pressure;
140  if (soundemitter[i + 1][j][k + 1] != 1) pre[i + 1][j][k + 1] = pressure;
141  if (soundemitter[i + 1][j + 1][k + 1] != 1) pre[i + 1][j + 1][k + 1] = pressure;
142  }
143  if (with_one_pulse == 0) {
144  if (soundemitter[i - 1][j - 1][k - 1] != 1) pre[i - 1][j - 1][k - 1] = pressure;
145  if (soundemitter[i - 1][j ][k - 1] != 1) pre[i - 1][j ][k - 1] = pressure;
146  if (soundemitter[i - 1][j + 1][k - 1] != 1) pre[i - 1][j + 1][k - 1] = pressure;
147  if (soundemitter[i - 1][j - 1][k] != 1) pre[i - 1][j - 1][k] = pressure;
148  if (soundemitter[i - 1][j][k] != 1) pre[i - 1][j][k] = pressure;
149  if (soundemitter[i - 1][j + 1][k] != 1) pre[i - 1][j + 1][k] = pressure;
150  if (soundemitter[i - 1][j - 1][k + 1] != 1) pre[i - 1][j - 1][k + 1] = pressure;
151  if (soundemitter[i - 1][j ][k + 1] != 1) pre[i - 1][j][k + 1] = pressure;
152  if (soundemitter[i - 1][j + 1][k + 1] != 1) pre[i - 1][j + 1][k + 1] = pressure;
153 
154  if (soundemitter[i][j - 1][k - 1] != 1) pre[i][j - 1][k - 1] = pressure;
155  if (soundemitter[i][j][k - 1] != 1) pre[i][j][k - 1] = pressure;
156  if (soundemitter[i][j + 1][k - 1] != 1) pre[i][j + 1][k - 1] = pressure;
157  if (soundemitter[i][j - 1][k] != 1) pre[i][j - 1][k] = pressure;
158  if (soundemitter[i][j + 1][k] != 1) pre[i][j + 1][k - 1] = pressure;
159  if (soundemitter[i][j - 1][k + 1] != 1) pre[i][j - 1][k + 1] = pressure;
160  if (soundemitter[i][j][k + 1] != 1) pre[i][j][k + 1] = pressure;
161  if (soundemitter[i][j + 1][k + 1] != 1) pre[i][j + 1][k + 1] = pressure;
162 
163  if (soundemitter[i + 1][j - 1][k - 1] != 1) pre[i + 1][j - 1][k - 1] = pressure;
164  if (soundemitter[i + 1][j][k - 1] != 1) pre[i + 1][j][k - 1] = pressure;
165  if (soundemitter[i + 1][j + 1][k - 1] != 1) pre[i + 1][j + 1][k - 1] = pressure;
166  if (soundemitter[i + 1][j - 1][k] != 1) pre[i + 1][j - 1][k] = pressure;
167  if (soundemitter[i + 1][j][k] != 1) pre[i + 1][j][k] = pressure;
168  if (soundemitter[i + 1][j + 1][k] != 1) pre[i + 1][j + 1][k] = pressure;
169  if (soundemitter[i + 1][j - 1][k + 1] != 1) pre[i + 1][j - 1][k + 1] = pressure;
170  if (soundemitter[i + 1][j][k + 1] != 1) pre[i + 1][j][k + 1] = pressure;
171  if (soundemitter[i + 1][j + 1][k + 1] != 1) pre[i + 1][j + 1][k + 1] = pressure;
172  }
173 
174  if ((dom[i][j][k] == 1) && (soundemitter[i][j][k] == 0)) {
175  if (dom[i - 1][j - 1][k - 1] != 1) pressure_on_solid[i - 1][j - 1][k - 1] = fabs(pre[i - 1][j - 1][k - 1] - pre_old[i - 1][j - 1][k - 1]);
176  if (dom[i - 1][j ][k - 1] != 1) pressure_on_solid[i - 1][j ][k - 1] = fabs(pre[i - 1][j ][k - 1] - pre_old[i - 1][j ][k - 1]);
177  if (dom[i - 1][j + 1][k - 1] != 1) pressure_on_solid[i - 1][j + 1][k - 1] = fabs(pre[i - 1][j + 1][k - 1] - pre_old[i - 1][j + 1][k - 1]);
178  if (dom[i - 1][j - 1][k] != 1) pressure_on_solid[i - 1][j - 1][k] = fabs(pre[i - 1][j - 1][k] - pre_old[i - 1][j - 1][k]);
179  if (dom[i - 1][j][k] != 1) pressure_on_solid[i - 1][j][k] = fabs(pre[i - 1][j][k] - pre_old[i - 1][j][k]);
180  if (dom[i - 1][j + 1][k] != 1) pressure_on_solid[i - 1][j + 1][k] = fabs(pre[i - 1][j + 1][k] - pre_old[i - 1][j + 1][k]);
181  if (dom[i - 1][j - 1][k + 1] != 1) pressure_on_solid[i - 1][j - 1][k + 1] = fabs(pre[i - 1][j - 1][k + 1] - pre_old[i - 1][j - 1][k + 1]);
182  if (dom[i - 1][j ][k + 1] != 1) pressure_on_solid[i - 1][j][k + 1] = fabs(pre[i - 1][j][k + 1] - pre_old[i - 1][j][k + 1]);
183  if (dom[i - 1][j + 1][k + 1] != 1) pressure_on_solid[i - 1][j + 1][k + 1] = fabs(pre[i - 1][j + 1][k + 1] - pre_old[i - 1][j + 1][k + 1]);
184 
185  if (dom[i][j - 1][k - 1] != 1) pressure_on_solid[i][j - 1][k - 1] = fabs(pre[i][j - 1][k - 1] - pre_old[i][j - 1][k - 1]);
186  if (dom[i][j][k - 1] != 1) pressure_on_solid[i][j][k - 1] = fabs(pre[i][j][k - 1] - pre_old[i][j][k - 1]);
187  if (dom[i][j + 1][k - 1] != 1) pressure_on_solid[i][j + 1][k - 1] = fabs(pre[i][j + 1][k - 1] - pre_old[i][j + 1][k - 1]);
188  if (dom[i][j - 1][k] != 1) pressure_on_solid[i][j - 1][k] = fabs(pre[i][j - 1][k] - pre_old[i][j - 1][k]);
189  if (dom[i][j + 1][k] != 1) pressure_on_solid[i][j + 1][k - 1] = fabs(pre[i][j + 1][k - 1] - pre_old[i][j + 1][k - 1]);
190  if (dom[i][j - 1][k + 1] != 1) pressure_on_solid[i][j - 1][k + 1] = fabs(pre[i][j - 1][k + 1] - pre_old[i][j - 1][k + 1]);
191  if (dom[i][j][k + 1] != 1) pressure_on_solid[i][j][k + 1] = fabs(pre[i][j][k + 1] - pre_old[i][j][k + 1]);
192  if (dom[i][j + 1][k + 1] != 1) pressure_on_solid[i][j + 1][k + 1] = fabs(pre[i][j + 1][k + 1] - pre_old[i][j + 1][k + 1]);
193 
194  if (dom[i + 1][j - 1][k - 1] != 1) pressure_on_solid[i + 1][j - 1][k - 1] = fabs(pre[i + 1][j - 1][k - 1] - pre_old[i + 1][j - 1][k - 1]);
195  if (dom[i + 1][j][k - 1] != 1) pressure_on_solid[i + 1][j][k - 1] = fabs(pre[i + 1][j][k - 1] - pre_old[i + 1][j][k - 1]);
196  if (dom[i + 1][j + 1][k - 1] != 1) pressure_on_solid[i + 1][j + 1][k - 1] = fabs(pre[i + 1][j + 1][k - 1] - pre_old[i + 1][j + 1][k - 1]);
197  if (dom[i + 1][j - 1][k] != 1) pressure_on_solid[i + 1][j - 1][k] = fabs(pre[i + 1][j - 1][k] - pre_old[i + 1][j - 1][k]);
198  if (dom[i + 1][j][k] != 1) pressure_on_solid[i + 1][j][k] = fabs(pre[i + 1][j][k] - pre_old[i + 1][j][k]);
199  if (dom[i + 1][j + 1][k] != 1) pressure_on_solid[i + 1][j + 1][k] = fabs(pre[i + 1][j + 1][k] - pre_old[i + 1][j + 1][k]);
200  if (dom[i + 1][j - 1][k + 1] != 1) pressure_on_solid[i + 1][j - 1][k + 1] = fabs(pre[i + 1][j - 1][k + 1] - pre_old[i + 1][j - 1][k + 1]);
201  if (dom[i + 1][j][k + 1] != 1) pressure_on_solid[i + 1][j][k + 1] = fabs(pre[i + 1][j][k + 1] - pre_old[i + 1][j][k + 1]);
202  if (dom[i + 1][j + 1][k + 1] != 1) pressure_on_solid[i + 1][j + 1][k + 1] = fabs(pre[i + 1][j + 1][k + 1] - pre_old[i + 1][j + 1][k + 1]);
203  }
204  }
205  }
206  }
207 
208 #ifdef _OPENMP
209  }
210 #endif
211 
212 
213 
214  //==========================================================Here the acoustic energy exchange===========================================================
215  //==========================================================Here the acoustic energy exchange===========================================================
216 
217 
218  local_pressure_right = 0.0;
219  local_pressure_left = 0.0;
220 
221 #ifdef _OPENMP
222 #pragma omp parallel default(none) \
223  private(j, i, k, local_pressure_left, \
224  local_pressure_right, one_boundary, \
225  left, right) \
226  shared(x, y, z, pre_old, pre, \
227  soundemitter, dom, pressure, \
228  pressure_on_solid, obstalce_absorption_coefficient, \
229  pressure_integrated, mean_pressure, counter, \
230  sound_reflexion)
231  {
232 #endif
233 
234  for (i = 1; i < x - 1; i++) {
235 #ifdef _OPENMP
236 #pragma omp for schedule(static)
237 #endif
238  for (j = 1; j < y - 1; j++) {
239  for (k = 1; k < z - 1; k++) {
240 
241  left = -1;
242  //now we simulate the obstacle acoustic-energy exchange
243  if ((dom[i][j][k] == 1) && (soundemitter[i][j][k] == 0)) {
244  if (dom[i - 1][j][k] != 1) {
245  local_pressure_left = pressure_on_solid[i - 1][j][k];
246  left = i - 1;
247  }
248  if ((dom[i + 1][j][k] != 1) && (left > 0)) {
249  local_pressure_right = pressure_on_solid[i + 1][j][k];
250  right = i + 1;
251  if (fabs(local_pressure_left) > fabs(local_pressure_right)) {
252  pre[right][j][k] += obstalce_absorption_coefficient*local_pressure_left;
253  pre[left][j][k] = pre_old[left][j][k] + (pressure_on_solid[left][j][k] * sound_reflexion);
254  }
255  if (fabs(local_pressure_left) < fabs(local_pressure_right)) {
256  pre[left][j][k] += obstalce_absorption_coefficient*local_pressure_right;
257  pre[right][j][k] = pre_old[right][j][k] + (pressure_on_solid[right][j][k] * sound_reflexion);
258  }
259  }
260  }
261  }
262  }
263  }
264 
265  for (i = 1; i < x - 1; i++) {
266 #ifdef _OPENMP
267 #pragma omp for schedule(static)
268 #endif
269  for (j = 1; j < y - 1; j++) {
270  for (k = 1; k < z - 1; k++) {
271  left = -1;
272  //now we simulate the obstacle acoustic-energy exchange
273  if ((dom[i][j][k] == 1) && (soundemitter[i][j][k] == 0)) {
274  if (dom[i][j - 1][k] != 1) {
275  local_pressure_left = pressure_on_solid[i][j - 1][k];
276  left = j - 1;
277  }
278 
279  if ((dom[i][j + 1][k] != 1) && (left > 0)) {
280  local_pressure_right = pressure_on_solid[i][j + 1][k];
281  right = j + 1;
282  if (fabs(local_pressure_left) > fabs(local_pressure_right)) {
283  pre[i][right][k] += obstalce_absorption_coefficient*local_pressure_left;
284  pre[i][left][k] = pre_old[i][left][k] + (pressure_on_solid[i][left][k] * sound_reflexion);
285 
286  }
287  if (fabs(local_pressure_left) < fabs(local_pressure_right)) {
288  pre[i][left][k] += obstalce_absorption_coefficient*local_pressure_right;
289  pre[i][right][k] = pre_old[i][right][k] + (pressure_on_solid[i][right][k] * sound_reflexion);
290  }
291  }
292  }
293  }
294  }
295  }
296 
297  for (i = 1; i < x - 1; i++) {
298 #ifdef _OPENMP
299 #pragma omp for schedule(static)
300 #endif
301  for (j = 1; j < y - 1; j++) {
302  for (k = 1; k < z - 1; k++) {
303  left = -1;
304  //now we simulate the obstacle acoustic-energy exchange
305  if ((dom[i][j][k] == 1) && (soundemitter[i][j][k] == 0)) {
306  if (dom[i][j][k - 1] != 1) {
307  local_pressure_left = pressure_on_solid[i][j][k - 1];
308  left = k - 1;
309  }
310  if ((dom[i][j][k + 1] != 1) && (left > 0)) {
311  local_pressure_right = pressure_on_solid[i][j][k + 1];
312  right = k + 1;
313  if (fabs(local_pressure_left) > fabs(local_pressure_right)) {
314  pre[i][j][right] += obstalce_absorption_coefficient*local_pressure_left;
315  pre[i][j][left] = pre_old[i][j][left] + (pressure_on_solid[i][j][left] * sound_reflexion);
316  }
317  if (fabs(local_pressure_left) < fabs(local_pressure_right)) {
318  pre[i][j][left] += obstalce_absorption_coefficient*local_pressure_right;
319  pre[i][j][right] = pre_old[i][j][right] + (pressure_on_solid[i][j][right] * sound_reflexion);
320  }
321  }
322  }
323  }
324  }
325  }
326 #ifdef _OPENMP
327  }
328 #endif
329 
330  //==========================================================Here the acoustic energy exchange END========================================================
331 
332  }
333 
334  if (dimension == 2) {
335 
336  k = 0;
337 
338 #ifdef _OPENMP
339 #pragma omp parallel default(none) \
340  private(j, i, k) \
341  shared(x, y, z, pre_old, pre, \
342  soundemitter, dom, pressure, \
343  pressure_on_solid, mean_pressure, \
344  with_one_pulse, counter)
345  {
346 #endif
347 
348 #ifdef _OPENMP
349 #pragma omp for schedule(static)
350 #endif
351  for (i = 1; i < x - 1; i++) {
352  for (j = 1; j < y - 1; j++) {
353 
354  if (soundemitter[i][j][k] == 1) {
355 
356  if ((with_one_pulse == 1) && (counter == 0)) {
357  if (soundemitter[i - 1][j - 1][k] != 1) pre[i - 1][j - 1][k] = pressure;
358  if (soundemitter[i - 1][j][k] != 1) pre[i - 1][j][k] = pressure;
359  if (soundemitter[i - 1][j + 1][k] != 1) pre[i - 1][j + 1][k] = pressure;
360 
361  if (soundemitter[i][j - 1][k] != 1) pre[i][j - 1][k] = pressure;
362  if (soundemitter[i][j + 1][k] != 1) pre[i][j + 1][k] = pressure;
363 
364 
365  if (soundemitter[i + 1][j - 1][k] != 1) pre[i + 1][j - 1][k] = pressure;
366  if (soundemitter[i + 1][j][k] != 1) pre[i + 1][j][k] = pressure;
367  if (soundemitter[i + 1][j + 1][k] != 1) pre[i + 1][j + 1][k] = pressure;
368  }
369  if (with_one_pulse == 0) {
370  if (soundemitter[i - 1][j - 1][k] != 1) pre[i - 1][j - 1][k] = pressure;
371  if (soundemitter[i - 1][j][k] != 1) pre[i - 1][j][k] = pressure;
372  if (soundemitter[i - 1][j + 1][k] != 1) pre[i - 1][j + 1][k] = pressure;
373 
374  if (soundemitter[i][j - 1][k] != 1) pre[i][j - 1][k] = pressure;
375  if (soundemitter[i][j + 1][k] != 1) pre[i][j + 1][k] = pressure;
376 
377 
378  if (soundemitter[i + 1][j - 1][k] != 1) pre[i + 1][j - 1][k] = pressure;
379  if (soundemitter[i + 1][j][k] != 1) pre[i + 1][j][k] = pressure;
380  if (soundemitter[i + 1][j + 1][k] != 1) pre[i + 1][j + 1][k] = pressure;
381  }
382  }
383 
384  if ((dom[i][j][k] == 1) && (soundemitter[i][j][k] == 0)) {
385 
386  if (dom[i - 1][j - 1][k] != 1) pressure_on_solid[i - 1][j - 1][k] = pre[i - 1][j - 1][k] - pre_old[i - 1][j - 1][k];
387  if (dom[i - 1][j][k] != 1) pressure_on_solid[i - 1][j][k] = pre[i - 1][j][k] - pre_old[i - 1][j][k];
388  if (dom[i - 1][j + 1][k] != 1) pressure_on_solid[i - 1][j + 1][k] = pre[i - 1][j + 1][k] - pre_old[i - 1][j + 1][k];
389 
390  if (dom[i][j - 1][k] != 1) pressure_on_solid[i][j - 1][k] = pre[i][j - 1][k] - pre_old[i][j - 1][k];
391  if (dom[i][j + 1][k] != 1) pressure_on_solid[i][j + 1][k] = pre[i][j + 1][k] - pre_old[i][j + 1][k];
392 
393 
394  if (dom[i + 1][j - 1][k] != 1) pressure_on_solid[i + 1][j - 1][k] = pre[i + 1][j - 1][k] - pre_old[i + 1][j - 1][k];
395  if (dom[i + 1][j][k] != 1) pressure_on_solid[i + 1][j][k] = pre[i + 1][j][k] - pre_old[i + 1][j][k];
396  if (dom[i + 1][j + 1][k] != 1) pressure_on_solid[i + 1][j + 1][k] = pre[i + 1][j + 1][k] - pre_old[i + 1][j + 1][k];
397  }
398 
399  }
400  }
401 #ifdef _OPENMP
402  }
403 #endif
404 
405 
406 
407  //==========================================================Here the acoustic energy exchange===========================================================
408  //==========================================================Here the acoustic energy exchange===========================================================
409 
410 
411  local_pressure_right = 0.0;
412  local_pressure_left = 0.0;
413  k = 0;
414 
415 
416 #ifdef _OPENMP
417 #pragma omp parallel default(none) \
418  private(j, i, k, local_pressure_left, \
419  local_pressure_right, \
420  left, right) \
421  shared(x, y, z, pre_old, pre, \
422  soundemitter, dom, pressure_on_solid, \
423  obstalce_absorption_coefficient, \
424  sound_reflexion)
425  {
426 #endif
427 
428 #ifdef _OPENMP
429 #pragma omp for schedule(static)
430 #endif
431  for (i = 1; i < x - 1; i++) {
432  for (j = 1; j < y - 1; j++) {
433  left = -1;
434  //now we simulate the obstacle acoustic-energy exchange
435  if ((dom[i][j][k] == 1) && (soundemitter[i][j][k] == 0)) {
436  if (dom[i - 1][j][k] != 1) {
437  local_pressure_left = pressure_on_solid[i - 1][j][k];
438  left = i - 1;
439 
440  }
441  if ((dom[i + 1][j][k] != 1) && (left > 0)) {
442  local_pressure_right = pressure_on_solid[i + 1][j][k];
443  right = i + 1;
444 
445  if (fabs(local_pressure_left) > fabs(local_pressure_right)) {
446  pre[right][j][k] += obstalce_absorption_coefficient*local_pressure_left;
447  pre[left][j][k] = pre_old[left][j][k] + (pressure_on_solid[left][j][k] * sound_reflexion);
448  }
449  if (fabs(local_pressure_left) < fabs(local_pressure_right)) {
450  pre[left][j][k] += obstalce_absorption_coefficient*local_pressure_right;
451  pre[right][j][k] = pre_old[right][j][k] + (pressure_on_solid[right][j][k] * sound_reflexion);
452  }
453  }
454  }
455  }
456  }
457 #ifdef _OPENMP
458  }
459 #endif
460 
461  local_pressure_right = 0.0;
462  local_pressure_left = 0.0;
463  k = 0;
464 
465 #ifdef _OPENMP
466 #pragma omp parallel default(none) \
467  private(j, i, k, local_pressure_left, \
468  local_pressure_right, \
469  left, right) \
470  shared(x, y, z, pre_old, pre, \
471  soundemitter, dom, pressure_on_solid, \
472  obstalce_absorption_coefficient, \
473  sound_reflexion)
474  {
475 #endif
476 
477 #ifdef _OPENMP
478 #pragma omp for schedule(static)
479 #endif
480  for (i = 1; i < x - 1; i++) {
481  for (j = 1; j < y - 1; j++) {
482 
483  left = -1;
484  //now we simulate the obstacle acoustic-energy exchange
485  if ((dom[i][j][k] == 1) && (soundemitter[i][j][k] == 0)) {
486  if (dom[i][j - 1][k] != 1) {
487  local_pressure_left = pressure_on_solid[i][j - 1][k];
488  left = j - 1;
489  }
490  if ((dom[i][j + 1][k] != 1) && (left > 0)) {
491  local_pressure_right = pressure_on_solid[i][j + 1][k];
492  right = j + 1;
493  if (fabs(local_pressure_left) > fabs(local_pressure_right)) {
494  pre[i][right][k] += obstalce_absorption_coefficient*local_pressure_left;
495  pre[i][left][k] = pre_old[i][left][k] + (pressure_on_solid[i][left][k] * sound_reflexion);
496  }
497  if (fabs(local_pressure_left) < fabs(local_pressure_right)) {
498  pre[i][left][k] += obstalce_absorption_coefficient*local_pressure_right;
499  pre[i][right][k] = pre_old[i][right][k] + (pressure_on_solid[i][right][k] * sound_reflexion);
500  }
501  }
502  }
503  }
504  }
505 #ifdef _OPENMP
506  }
507 #endif
508 
509 
510  //==========================================================Here the acoustic energy exchange END========================================================
511 
512 
513  }
514 
515 
516  local_pressure_right = 0.0;
517  local_pressure_left = 0.0;
518  k = 0;
519  j = 0;
520 
521  if (dimension == 1) {
522  for (i = 1; i < x - 1; i++) {
523  if ((with_one_pulse == 1) && (counter == 0)) {
524  if (soundemitter[i - 1][j][k] != 1) pre[i - 1][j][k] = pressure;
525  if (soundemitter[i + 1][j][k] != 1) pre[i + 1][j][k] = pressure;
526  }
527  if (with_one_pulse == 0) {
528  if (soundemitter[i - 1][j][k] != 1) pre[i - 1][j][k] = pressure;
529  if (soundemitter[i + 1][j][k] != 1) pre[i + 1][j][k] = pressure;
530  }
531  if ((dom[i][j][k] == 1) && (soundemitter[i][j][k] == 0)) {
532  if (dom[i - 1][j][k] != 1) pressure_on_solid[i - 1][j][k] = fabs(pre[i - 1][j][k] - pre_old[i - 1][j][k]);
533  if (dom[i + 1][j][k] != 1) pressure_on_solid[i + 1][j][k] = fabs(pre[i + 1][j][k] - pre_old[i + 1][j][k]);
534  }
535  }
536 
537  for (i = 1; i < x - 1; i++) {
538  left = -1;
539  //now we simulate the obstacle acoustic-energy exchange
540  if ((dom[i][j][k] == 1) && (soundemitter[i][j][k] == 0)) {
541  if (dom[i - 1][j][k] != 1) {
542  local_pressure_left = pressure_on_solid[i - 1][j][k];
543  left = i - 1;
544  }
545  if ((dom[i + 1][j][k] != 1) && (left > 0)) {
546  local_pressure_right = pressure_on_solid[i + 1][j][k];
547  right = i + 1;
548  if (fabs(local_pressure_left) > fabs(local_pressure_right)) {
549  pre[right][j][k] += obstalce_absorption_coefficient*local_pressure_left;
550  pre[left][j][k] = pre_old[left][j][k] + (pressure_on_solid[left][j][k] * sound_reflexion);
551  }
552  if (fabs(local_pressure_left) < fabs(local_pressure_right)) {
553  pre[left][j][k] += obstalce_absorption_coefficient*local_pressure_right;
554  pre[right][j][k] = pre_old[right][j][k] + (pressure_on_solid[right][j][k] * sound_reflexion);
555  }
556  }
557  }
558  }
559  }
560 
561 
562  //==========================================================Here parts of the dB Map are calculated============================================================
563 
564  //here the preparation for the periodic dB - Map is done
565  if (with_one_pulse == 0) {
566 
567 #ifdef _OPENMP
568 #pragma omp parallel default(none) \
569  private(j, i, k, local_pressure_left, \
570  local_pressure_right, one_boundary, \
571  left, right) \
572  shared(x, y, z, pre, dom, \
573  pressure_integrated, mean_pressure, dt)
574  {
575 #endif
576 
577  for (i = 0; i < x; i++) {
578 #ifdef _OPENMP
579 #pragma omp for schedule(static)
580 #endif
581  for (j = 0; j < y; j++) {
582  for (k = 0; k < z; k++) {
583 
584  if (dom[i][j][k] == 0) {
585  pressure_integrated[i][j][k] += pow((mean_pressure - pre[i][j][k]), 2) * dt;
586  }
587  }
588  }
589  }
590 
591 #ifdef _OPENMP
592  }
593 #endif
594  }
595 
596  //here the preparation for the pulse dB - Map is done
597  if (with_one_pulse == 1) {
598 
599 #ifdef _OPENMP
600 #pragma omp parallel default(none) \
601  private(j, i, k, local_pressure_left, \
602  local_pressure_right, one_boundary, \
603  left, right, tmp) \
604  shared(x, y, z, pre, pre_old, dom, \
605  pressure_integrated, mean_pressure, dt, \
606  perfect_sinus_emitter)
607  {
608 #endif
609 
610  for (i = 0; i < x; i++) {
611 #ifdef _OPENMP
612 #pragma omp for schedule(static)
613 #endif
614  for (j = 0; j < y; j++) {
615  for (k = 0; k < z; k++) {
616  if (dom[i][j][k] == 0) {
617 
618  tmp = fabs(pre_old[i][j][k] - pre[i][j][k]) * perfect_sinus_emitter;
619  if (tmp < 2.0e-5) tmp = 0.0;
620  if (tmp > pressure_integrated[i][j][k]) pressure_integrated[i][j][k] = tmp;
621  }
622  }
623  }
624  }
625 #ifdef _OPENMP
626  }
627 #endif
628  }
629 
630  return 0;
631 
632 }
633 
634 int reset_pressure_integrated(int x, int y, int z) {
635 
636  int i, j, k;
637 
638 #ifdef _OPENMP
639 #pragma omp parallel default(none) \
640  private(j, i, k) \
641  shared(x, y, z, pressure_integrated)
642  {
643 #endif
644 
645  for (i = 0; i < x; i++) {
646 #ifdef _OPENMP
647 #pragma omp for schedule(static)
648 #endif
649 
650  for (j = 0; j < y; j++) {
651  for (k = 0; k < z; k++) {
652 
653  pressure_integrated[i][j][k] = 0.0;
654  }
655  }
656  }
657 
658 #ifdef _OPENMP
659  }
660 #endif
661 
662  return 0;
663 
664 }
665 
666 int prepare_the_dB_map(int x, int y, int z) {
667 
668  int i, j, k;
669 
670 #ifdef _OPENMP
671 #pragma omp parallel default(none) \
672  private(j, i, k) \
673  shared(x, y, z, pressure_integrated, \
674  dB_map, mean_pressure, dt_integrated, \
675  with_one_pulse, dB_Map_ready)
676  {
677 #endif
678 
679  //just to make sure 0.0 dB-Maps are not written out
680  dB_Map_ready = 0;
681 
682 #ifdef _OPENMP
683 #pragma omp for schedule(static)
684 #endif
685  for (i = 0; i < x; i++) {
686  for (j = 0; j < y; j++) {
687  for (k = 0; k < z; k++) {
688  if (pressure_integrated[i][j][k] > 0.0) {
689  if (with_one_pulse == 0) dB_map[i][j][k] = 20 * log10((sqrt((1.0 / dt_integrated) * pressure_integrated[i][j][k])) / 2.0E-5);
690  //Assumption perfect sinus sound emitter
691  if (with_one_pulse == 1) dB_map[i][j][k] = 20 * log10(pressure_integrated[i][j][k] / 2.0E-5);
692 
693  dB_Map_ready = 1;
694 
695  } else {
696  dB_map[i][j][k] = 0.0;
697  }
698  }
699  }
700  }
701 
702 #ifdef _OPENMP
703  }
704 #endif
705 
706  return 0;
707 }
Generated on Thu Oct 10 2013 17:15:52 for TYCHO by   doxygen 1.8.1.1