Code_TYMPAN  4.2.0
Industrial site acoustic simulation
entities.cpp
Go to the documentation of this file.
1 
9 #include <vector>
10 #include <numeric>
11 #include "entities.hpp"
12 #include "Tympan/models/common/atmospheric_conditions.h"
13 
14 namespace tympan
15 {
16  AcousticMaterialBase::AcousticMaterialBase(
17  const string& name_
18  ) : name(name_) {}
19 
20  // ---------
21  AtmosphericConditions *AcousticGroundMaterial::atmosphere = new AtmosphericConditions(101325., 20., 70.);
22 
23  AcousticBuildingMaterial::AcousticBuildingMaterial(
24  const string& name_,
25  const ComplexSpectrum& spectrum_
26  ) : AcousticMaterialBase(name_), spectrum(spectrum_) {}
27 
28  ComplexSpectrum AcousticBuildingMaterial::asEyring() const
29  {
30  ComplexSpectrum eyring(spectrum);
31  // 1-e^(-alphaS)
32 
33  return eyring.exp(-1.0).mult(-1.0).sum(1.0);
34  }
35 
36  // ---------
37  AcousticGroundMaterial::AcousticGroundMaterial(
38  const string& name_, double resistivity_, double deviation_, double length_) :
39  AcousticMaterialBase(name_), resistivity(resistivity_), thickness(1.0), deviation(deviation_), length(length_)
40  {
41  init();
42  }
43 
44  void AcousticGroundMaterial::init()
45  {
46  computeZc();
47  computeK();
48  }
49 
50  ComplexSpectrum AcousticGroundMaterial::get_absorption (double incidence_angle, double length)
51  {
52  ComplexSpectrum Zs, Rp, W, Fw, Q;
53  computeZs(incidence_angle, Zc, Zs);
54  computeZf(incidence_angle, Zs);
55  computeRp(incidence_angle, Zf, Rp);
56  computeW(incidence_angle, length, Zf, W);
57  computeFw(W, Fw);
58  computeQ(incidence_angle, Rp, Fw, Q);
59 
60  Q = Q.toModulePhase();
61  Q.setType(SPECTRE_TYPE_ABSO);
62 
63  return Q;
64  }
65 
66 
67  double AcousticGroundMaterial::get_ISO9613_G()
68  {
69  return MIN( pow(300.0 / resistivity, 0.57), 1.0 );
70  }
71 
72  void AcousticGroundMaterial::computeZc()
73  {
74  double sigmaSurF;
75  OTabFreq tabFreq = OSpectre::getTabFreqExact(); // On travaille en frequence reelle
76 
77  for (unsigned int i = 0 ; i < Zc.getNbValues() ; i++)
78  {
79  sigmaSurF = resistivity / tabFreq[i];
80 
81  Zc.getTabValReel()[i] = 1.0 + 9.08 * pow(sigmaSurF, 0.75);
82  Zc.getTabValImag()[i] = 11.9 * pow(sigmaSurF, 0.73);
83  }
84 
85  Zc.setType(SPECTRE_TYPE_AUTRE);
86  }
87 
88  void AcousticGroundMaterial::computeK()
89  {
90  double k_value, FSurSigma;
91 
92  OTabFreq tabFreq = OSpectre::getTabFreqExact(); // On travaille en frequence reelle
93 
94  for (unsigned int i = 0 ; i < K.getNbValues() ; i++)
95  {
96  FSurSigma = tabFreq[i] / resistivity;
97 
98  k_value = atmosphere->get_k().getTabValReel()[i];
99 
100  K.getTabValReel()[i] = k_value * (1 + 10.8 * pow(FSurSigma, -0.7));
101  K.getTabValImag()[i] = k_value * (10.3 * pow(FSurSigma, -0.59));
102  }
103 
104  K.setType(SPECTRE_TYPE_AUTRE);
105  }
106 
107  void AcousticGroundMaterial::computeZs(double angle, ComplexSpectrum Z, ComplexSpectrum& spectrum)
108  {
109  double k_value;
110 
111  TYComplex operand, K_value, Z_value;
112  double cosPhi = cos(angle);
113 
114  TYComplex cplxVal;
115 
116  for (unsigned int i = 0 ; i < Z.getNbValues(); i++)
117  {
118  k_value = atmosphere->get_k().getTabValReel()[i];
119  K_value = TYComplex(K.getTabValReel()[i], K.getTabValImag()[i]);
120  Z_value = TYComplex(Z.getTabValReel()[i], Z.getTabValImag()[i]);
121 
122  operand = std::sqrt(CPLX_UN - ((k_value * cosPhi / K_value) * (k_value * cosPhi / K_value)));
123 
124  cplxVal = Z_value * cotanh(CPLX_MUN * CPLX_J * K_value * thickness * operand) / operand;
125  spectrum.getTabValReel()[i] = cplxVal.real();
126  spectrum.getTabValImag()[i] = cplxVal.imag();
127  }
128  }
129 
130  void AcousticGroundMaterial::computeZf(double angle, ComplexSpectrum Zs)
131  {
132  double k0_value, k_value, intv_a, intv_b;
133  size_t size = Zf.getNbValues();
134  ComplexSpectrum Bf;
135 
136  for (unsigned int i=0; i < size; i++)
137  {
138  k0_value = atmosphere->get_k().getTabValReel()[i];
139  k_value = k0_value*sin(angle);
140 
141  // Partie réelle
142  intv_a = sqrt(k0_value)/100;
143  std::vector<double> u_alpha(size),integrande_alpha1(size),integrande_alpha2(size);
144  for (int j=0; j<=100; j++)
145  {
146  u_alpha.push_back(j*intv_a);
147  double u_value = k0_value - u_alpha[j]*u_alpha[j];
148  double u_2value = 2*k0_value - u_alpha[j]*u_alpha[j];
149  double expression1a = 1.0/(k0_value*sqrt(u_2value))\
150  *((k0_value*k0_value+k_value*u_value)*(k0_value*k0_value+k_value*u_value))\
151  *gaussianSpectrum(k_value+u_value, deviation, length);
152  double expression2a = 1.0/(k0_value*sqrt(u_2value))\
153  *((k0_value*k0_value-k_value*u_value)*(k0_value*k0_value-k_value*u_value))\
154  *gaussianSpectrum(k_value-u_value, deviation, length);
155  integrande_alpha1.push_back(expression1a);
156  integrande_alpha2.push_back(expression2a);
157  }
158  double alpha1 = trapz(u_alpha, integrande_alpha1) ;
159  double alpha2 = trapz(u_alpha, integrande_alpha2) ;
160  double alpha = alpha1 + alpha2;
161  Bf.getTabValReel()[i] = alpha;
162  u_alpha.clear();
163  integrande_alpha1.clear();
164  integrande_alpha2.clear();
165 
166  // Partie imaginaire
167  intv_b = (6/length)/100;// Condition sur length
168  std::vector<double> u_beta(size),integrande_beta1(size),integrande_beta2(size);
169  for (int j=0; j<=100; j++)
170  {
171  u_beta.push_back(j*intv_b);
172  double u_value = k0_value*k0_value+u_beta[j]*u_beta[j];
173  double expression1b = 1.0/(k0_value*sqrt(u_value))\
174  *(k0_value*k0_value+k_value*sqrt(u_value))*(k0_value*k0_value+k_value*sqrt(u_value))\
175  *gaussianSpectrum(k_value+sqrt(u_value), deviation, length);
176  double expression2b = 1.0/(k0_value*sqrt(u_value))\
177  *(k0_value*k0_value-k_value*sqrt(u_value))*(k0_value*k0_value-k_value*sqrt(u_value))\
178  *gaussianSpectrum(k_value-sqrt(u_value), deviation, length);
179  integrande_beta1.push_back(expression1b);
180  integrande_beta2.push_back(expression2b);
181  }
182  double beta1 = -trapz(u_beta, integrande_beta1) ;
183  double beta2 = -trapz(u_beta, integrande_beta2) ;
184  double beta = beta1 + beta2;
185  Bf.getTabValImag()[i] = beta;
186  u_beta.clear();
187  integrande_beta1.clear();
188  integrande_beta2.clear();
189 
190  // Récupération de Zs
191  TYComplex cplxVal;
192  TYComplex Zs_value = TYComplex(Zs.getTabValReel()[i], Zs.getTabValImag()[i]);
193  cplxVal = CPLX_UN/Zs_value;
194 
195  // Calcul de l'impedance effective
196  Bf.getTabValReel()[i] += cplxVal.real();
197  Bf.getTabValImag()[i] += cplxVal.imag();
198 
199  TYComplex ZcplxVal;
200  TYComplex Bf_value = TYComplex(Bf.getTabValReel()[i], Bf.getTabValImag()[i]);
201  ZcplxVal = CPLX_UN/Bf_value;
202 
203  Zf.getTabValReel()[i] = ZcplxVal.real();
204  Zf.getTabValImag()[i] = ZcplxVal.imag();
205  }
206  Zf.setType(SPECTRE_TYPE_AUTRE);
207  }
208 
209  double AcousticGroundMaterial::gaussianSpectrum(double k, double sigma, double lc)
210  {
211  return (sigma*sigma)*(lc/(2*sqrt(M_PI)))*exp(-1.0/4*((k*lc)*(k*lc)));
212  }
213 
214  double AcousticGroundMaterial::trapz(std::vector<double> u, std::vector<double> integrande)
215  {
216  unsigned int size = u.size();
217  double sum = std::accumulate(integrande.begin()+1,integrande.end()-1,(integrande.front()+integrande.back())/2);
218  sum *= (u.back()-u.front())/size;
219  return sum;
220  }
221 
222  void AcousticGroundMaterial::computeRp(double angle, const ComplexSpectrum& Zs, ComplexSpectrum& Rp)
223  {
224  // const TYComplex CPLX_UN = TYComplex(1.0, 0.0);
225  double sinPhi = sin(angle);
226  TYComplex Z;
227 
228  for (unsigned int i = 0 ; i < Zs.getNbValues() ; i++)
229  {
230  Z = TYComplex(Zs.getTabValReel()[i], Zs.getTabValImag()[i]) ;
231 
232  TYComplex val = (Z * sinPhi - CPLX_UN) / (Z * sinPhi + CPLX_UN);
233 
234  Rp.getTabValReel()[i] = val.real();
235  Rp.getTabValImag()[i] = val.imag();
236  }
237  }
238 
239  void AcousticGroundMaterial::computeW(double angle, double length, const ComplexSpectrum& Zs, ComplexSpectrum &W)
240  {
241  TYComplex zs_value, invZs_value, cplxTempo, wTemp;
242 
243  double k_value; // nombre d'onde acoustique
244  double sinPhi = sin(angle);
245 
246  for (unsigned int i = 0 ; i < W.getNbValues() ; i++)
247  {
248  k_value = atmosphere->get_k().getTabValReel()[i];
249  zs_value = TYComplex(Zs.getTabValReel()[i], Zs.getTabValImag()[i]);
250 
251  invZs_value = 1.0 / zs_value;
252  cplxTempo = (sinPhi + invZs_value) * (sinPhi + invZs_value);
253 
254  wTemp = std::sqrt(1.0 / 2.0 * CPLX_J * k_value * length * cplxTempo) ;
255 
256  W.getTabValReel()[i] = wTemp.real();
257  W.getTabValImag()[i] = wTemp.imag();
258  }
259  }
260 
261  void AcousticGroundMaterial::computeFw(ComplexSpectrum localW, ComplexSpectrum& Fw)
262  {
263  // localW is intentionnaly a copy of w (values may change)
264 
265  OSpectreComplex G;
266 
267  // Recherche et remplacement des parties reelles et imaginaires < epsilon par epsilon
268  limit_W_values(localW);
269 
270  // Filtrage des differents cas
271  erfc_G_computation(localW, G);
272 
273  // Operations en fonction du signe des parties reelles et imaginaires de w
274  sgn_G_computation(localW, G);
275 
276  // Calcul de la fonction Fw
277  double sqrt_pi = sqrt(M_PI);
278 
279  for (size_t i = 0; i < localW.getNbValues(); i++)
280  {
281  TYComplex v1 = TYComplex(localW.getTabValReel()[i], localW.getTabValImag()[i]);
282  TYComplex v2 = TYComplex(G.getTabValReel()[i], G.getTabValImag()[i]);
283  TYComplex val = CPLX_UN + CPLX_J * sqrt_pi * v1 * v2;
284  Fw.getTabValReel()[i] = val.real();
285  Fw.getTabValImag()[i] = val.imag();
286  }
287  }
288 
289  void AcousticGroundMaterial::limit_W_values(ComplexSpectrum& localW)
290  {
291  const double epsilon = 1.0E-20;
292  double valeur, tampon;
293 
294  for(size_t i = 0; i < localW.getNbValues() ; i++)
295  {
296  // Si partie reelle < epsilon alors partie reelle = epsilon
297  tampon = localW.getTabValReel()[i];
298  valeur = (fabs(tampon) < epsilon ? epsilon : tampon) ;
299  localW.getTabValReel()[i] = valeur;
300 
301  // Si partie imaginaire < epsilon alors partie imaginaire = epsilon
302  tampon = localW.getTabValImag()[i];
303  valeur = (fabs(tampon) < epsilon ? epsilon : tampon) ;
304  localW.getTabValImag()[i] = valeur;
305  }
306  }
307 
308  void AcousticGroundMaterial::erfc_G_computation(const ComplexSpectrum& localW, ComplexSpectrum& G)
309  {
310  double reelW, imagW;
311  TYComplex val;
312 
313  for (size_t i = 0 ; i < localW.getNbValues() ; i++)
314  {
315  reelW = localW.getTabValReel()[i];
316  imagW = localW.getTabValImag()[i];
317 
318  if (((fabs(reelW) < 3.9) && (fabs(imagW) < 3.0))) // Cas 1
319  {
320  val = erfcCas1(TYComplex(reelW, imagW));
321  }
322  else if (((fabs(reelW) >= 3.9) && (fabs(reelW) < 6.0)) || // Cas 2
323  ((fabs(imagW) >= 3.0) && (fabs(imagW) < 6.0)))
324  {
325  val = erfcCas2(TYComplex(fabs(reelW), fabs(imagW)));
326  }
327  else // Cas 3
328  {
329  val = erfcCas3(TYComplex(fabs(reelW), fabs(imagW)));
330  }
331 
332  G.getTabValReel()[i] = val.real();
333  G.getTabValImag()[i] = val.imag();
334  }
335  }
336 
337  void AcousticGroundMaterial::sgn_G_computation(const ComplexSpectrum& localW, ComplexSpectrum& G)
338  {
339  double re, im;
340  for (size_t i = 0 ; i < localW.getNbValues() ; i++)
341  {
342  re = localW.getTabValReel()[i];
343  im = localW.getTabValImag()[i];
344 
345  TYComplex conjG = conj(TYComplex(G.getTabValReel()[i], G.getTabValImag()[i]));
346 
347  if ((re < 0) && (im > 0)) // imp
348  {
349  G.getTabValReel()[i] = conjG.real() ;
350  G.getTabValImag()[i] = conjG.imag() ;
351  }
352  else if ((re > 0) && (im < 0)) //ipm
353  {
354  TYComplex v1 = TYComplex(localW.getTabValReel()[i], localW.getTabValImag()[i]);
355  TYComplex v2 = TYComplex(G.getTabValReel()[i], G.getTabValImag()[i]);
356  TYComplex tampon = sgnReIm(v1, v2);
357 
358  G.getTabValReel()[i] = tampon.real();
359  G.getTabValImag()[i] = tampon.imag(); ;
360  }
361  else if ((re < 0) && (im < 0)) // imm
362  {
363  TYComplex v1 = TYComplex(localW.getTabValReel()[i], localW.getTabValImag()[i]);
364  TYComplex v2 = TYComplex(G.getTabValReel()[i], G.getTabValImag()[i]);
365  TYComplex tampon = sgnReIm(v1, v2);
366 
367  TYComplex val = conj(tampon);
368 
369  G.getTabValReel()[i] = val.real();
370  G.getTabValImag()[i] = val.imag();
371  }
372  else if ((re > 0) && (im > 0)) //ipp
373  {
374  }
375  else
376  {
377  }
378  }
379  }
380 
381  void AcousticGroundMaterial::computeQ(double angle, ComplexSpectrum &Rp, ComplexSpectrum &Fw, ComplexSpectrum &Q)
382  {
383  TYComplex Rp_value, Fw_value ;
384 
385  for (unsigned int i = 0 ; i < Fw.getNbValues(); i++)
386  {
387  Rp_value = TYComplex(Rp.getTabValReel()[i], Rp.getTabValImag()[i]);
388  Fw_value = TYComplex(Fw.getTabValReel()[i], Fw.getTabValImag()[i]);
389 
390  TYComplex val = Rp_value + (CPLX_UN - Rp_value) * Fw_value;
391 
392  Q.getTabValReel()[i] = val.real();
393  Q.getTabValImag()[i] = val.imag();
394  }
395  }
396 
397  TYComplex AcousticGroundMaterial::erfcCas1(const TYComplex& wValue) const
398  {
399  const double pi = M_PI;
400  const int nbIterations = 5; // Valeur normale = 5
401  const double h = 0.8; // valeur proposee = 0.8;
402  const double h2 = h * h;
403 
404  double x = fabs(wValue.real());
405  double y = fabs(wValue.imag());
406  double x2 = x * x;
407  double y2 = y * y;
408  double modW = x2 + y2;
409 
410  double A1 = cos(2.0 * x * y);
411  double B1 = sin(2.0 * x * y);
412  double C1 = exp(-2.0 * pi * y / h) - cos(2.0 * pi * x / h) ;
413  double D1 = sin(2.0 * pi * x / h);
414  double P2 = 2.0 * exp(-(x2 + 2.0 * pi * y / h - y2)) * (A1 * C1 - B1 * D1) / (C1 * C1 + D1 * D1);
415  double Q2 = 2.0 * exp(-(x2 + 2.0 * pi * y / h - y2)) * (A1 * D1 + B1 * C1) / (C1 * C1 + D1 * D1);
416 
417  double H0, H, K0, K;
418 
419  if (y < (pi / h))
420  {
421  H0 = h * y / (pi * modW) + P2;
422  }
423  else if (y == (pi / h))
424  {
425  H0 = h * y / (pi * modW) + (P2 / 2.0);
426  }
427  else
428  {
429  H0 = h * y / (pi * modW);
430  }
431 
432  if (x < (pi / h))
433  {
434  K0 = h * x / (pi * modW) - Q2;
435  }
436  else if (x == (pi / h))
437  {
438  K0 = h * x / (pi * modW) - (Q2 / 2.0) ;
439  }
440  else
441  {
442  K0 = h * x / (pi * modW);
443  }
444 
445  /* Solution MATLAB */
446  H = H0;
447  K = K0;
448  double n2 = 0.0;
449  double coeff_y = 2.0 * y * h / pi;
450  double coeff_x = 2.0 * x * h / pi;
451  double diviseur = 0.0;
452  double expMn2h2 = 0.0;
453  for (int n = 1 ; n <= nbIterations; n++)
454  {
455  n2 = (double)(n * n);
456  diviseur = ((y2 - x2 + n2 * h2) * (y2 - x2 + n2 * h2) + (4.0 * y2 * x2));
457  expMn2h2 = exp(-n2 * h2);
458 
459  H = H + coeff_y * expMn2h2 * (modW + n2 * h2) / diviseur;
460  K = K + coeff_x * expMn2h2 * (modW - n2 * h2) / diviseur;
461  }
462 
463  return TYComplex(H, K);
464  }
465 
466  TYComplex AcousticGroundMaterial::erfcCas2(const TYComplex& wValue) const
467  {
468  TYComplex G ;
469  TYComplex w2 = wValue * wValue;
470 
471  G = CPLX_J * wValue * ((0.4613135 / (w2 - 0.1901635)) +
472  (0.09999216 / (w2 - 1.7844927)) +
473  (0.002883894 / (w2 - 5.5253437))) ;
474 
475  return G;
476  }
477 
478  TYComplex AcousticGroundMaterial::erfcCas3(const TYComplex& wValue) const
479  {
480  TYComplex G ;
481 
482  TYComplex w2 = wValue * wValue;
483 
484  G = CPLX_J * wValue * ((0.5124242 / (w2 - 0.2752551)) +
485  (0.05176536 / (w2 - 2.724745))) ;
486 
487  return G;
488  }
489 
490  TYComplex AcousticGroundMaterial::sgnReIm(const TYComplex& W, const TYComplex& G) const
491  {
492  double x = fabs(W.real());
493  double y = fabs(W.imag());
494 
495  double cos2xy = cos(2 * x * y);
496  double sin2xy = sin(2 * x * y);
497 
498  return 2.0 * exp(y * y - x * x) * (cos2xy + CPLX_J * sin2xy) - conj(G);
499  }
500 
501  // ---------
502 
503  AtmosphericConditions *CommonFaceDirectivity::atmosphere = new AtmosphericConditions(101325., 20., 70.);
504 
505  // ---------
506  /*static*/ const double VolumeFaceDirectivity::_tabRA[] = { 1.0,
507  2.0,
508  3.0,
509  4.0,
510  5.0,
511  10.0,
512  20.0,
513  30.0,
514  40.0,
515  50.0,
516  100.0,
517  200.0,
518  300.0,
519  };
520  /*static*/ const double VolumeFaceDirectivity::_tabCor[] = {0.608,
521  0.817,
522  0.879,
523  0.909,
524  0.928,
525  0.964,
526  0.982,
527  0.988,
528  0.991,
529  0.993,
530  0.996,
531  0.998,
532  1.000
533  };
534  double VolumeFaceDirectivity::calculC(double distance)
535  {
536  double X = distance / support_size;
537 
538  if (X < 1.0) { return _tabCor[0]; }
539  if (X > 300.0) { return 1.0; }
540 
541  int i = 1;
542  while (_tabRA[i] < X && i < 12) { i++; } // Recherche de l'indice le plus proche
543 
544  // Calcul de la valeur de c par interpolation lineaire
545  return (_tabCor[i] - _tabCor[i - 1]) / (_tabRA[i] - _tabRA[i - 1]) * (X - _tabRA[i - 1]) + _tabCor[i - 1];
546  }
547 
548  Spectrum VolumeFaceDirectivity::lwAdjustment(Vector direction, double distance)
549  {
550  Spectrum directivity_spectrum;
551  double q, ka;
552 
553  // Angle angle between support normal and (source to receptor) vector
554  double cosPhi = cos( support_normal.angle( direction ) );
555 
556  double C = calculC(distance); //Coefficient de directivite
557 
558  for (unsigned int i = 0 ; i < directivity_spectrum.getNbValues() ; i++)
559  {
560  ka = atmosphere->get_k().getTabValReel()[i] * support_size ;
561 
562  q = (1 + (ka / (ka + 0.4)) * cosPhi) * C;
563 
564  directivity_spectrum.getTabValReel()[i] = q;
565  }
566 
567  return directivity_spectrum;
568  }
569  //
570  // ----------------------------------------
571  //
572 #ifdef NB_KA
573 # undef NB_KA
574 # define NB_KA 38
575 #else
576 # define NB_KA 38
577 #endif
578 
579 #ifdef NB_THETA
580 # undef NB_THETA
581 # define NB_THETA 21
582 #else
583 # define NB_THETA 21
584 #endif
585  const double ChimneyFaceDirectivity::_tabQ[NB_KA][NB_THETA] =
586  {
587  {
588  /*ka1*/ 3.9213E-02, 3.8299E-02, 3.5626E-02, 3.1399E-02, 2.5933E-02, 1.9633E-02, 1.2955E-02,
589  6.3631E-03, 2.9321E-04, -4.8880E-03, -8.9124E-03, -1.1631E-02, -1.3023E-02, -1.3192E-02,
590  -1.2351E-02, -1.0793E-02, -8.8587E-03, -6.8993E-03, -5.2361E-03, -4.1269E-03, -3.7380E-03
591  },
592  {
593  /*ka2*/ 1.5295E-01, 1.4940E-01, 1.3903E-01, 1.2260E-01, 1.0132E-01, 7.6753E-02, 5.0654E-02,
594  2.4837E-02, 1.0112E-03, -1.9370E-02, -3.5229E-02, -4.5960E-02, -5.1463E-02, -5.2134E-02,
595  -4.8805E-02, -4.2637E-02, -3.4980E-02, -2.7223E-02, -2.0640E-02, -1.6250E-02, -1.4711E-02
596  },
597  {
598  /*ka3*/ 3.3341E-01, 3.2573E-01, 3.0323E-01, 2.6755E-01, 2.2126E-01, 1.6765E-01, 1.1054E-01,
599  5.3866E-02, 1.4017E-03, -4.3595E-02, -7.8677E-02, -1.0242E-01, -1.1455E-01, -1.1590E-01,
600  -1.0831E-01, -9.4383E-02, -7.7136E-02, -5.9688E-02, -4.4893E-02, -3.5032E-02, -3.1576E-02
601  },
602  {
603  /*ka4*/ 5.7143E-01, 5.5835E-01, 5.2003E-01, 4.5914E-01, 3.7992E-01, 2.8790E-01, 1.8948E-01,
604  9.1440E-02, 3.3045E-04, -7.8073E-02, -1.3934E-01, -1.8081E-01, -2.0183E-01, -2.0382E-01,
605  -1.8997E-01, -1.6488E-01, -1.3394E-01, -1.0271E-01, -7.6267E-02, -5.8657E-02, -5.2487E-02
606  },
607  {
608  /*ka5*/ 8.5759E-01, 8.3811E-01, 7.8096E-01, 6.8997E-01, 5.7125E-01, 4.3285E-01, 2.8419E-01,
609  1.3542E-01, -3.4683E-03, -1.2348E-01, -2.1752E-01, -2.8120E-01, -3.1323E-01, -3.1559E-01,
610  -2.9322E-01, -2.5329E-01, -2.0429E-01, -1.5494E-01, -1.1321E-01, -8.5459E-02, -7.5741E-02
611  },
612  {
613  /*ka6*/ 1.1831E+00, 1.1564E+00, 1.0781E+00, 9.5307E-01, 7.8950E-01, 5.9808E-01, 3.9154E-01,
614  1.8377E-01, -1.1197E-02, -1.8049E-01, -3.1367E-01, -4.0394E-01, -4.4902E-01, -4.5143E-01,
615  -4.1811E-01, -3.5943E-01, -2.8771E-01, -2.1567E-01, -1.5487E-01, -1.1446E-01, -1.0033E-01
616  },
617  {
618  /*ka7*/ 1.5403E+00, 1.5058E+00, 1.4044E+00, 1.2424E+00, 1.0296E+00, 7.7974E-01, 5.0881E-01,
619  2.3480E-01, -2.3835E-02, -2.4970E-01, -4.2827E-01, -5.4964E-01, -6.0999E-01, -6.1221E-01,
620  -5.6552E-01, -4.8406E-01, -3.8483E-01, -2.8537E-01, -2.0156E-01, -1.4594E-01, -1.2649E-01
621  },
622  {
623  /*ka8*/ 1.9232E+00, 1.8805E+00, 1.7546E+00, 1.5530E+00, 1.2876E+00, 9.7474E-01, 6.3391E-01,
624  2.8725E-01, -4.2025E-02, -3.3151E-01, -5.6184E-01, -7.1921E-01, -7.9754E-01, -7.9976E-01,
625  -7.3751E-01, -6.2935E-01, -4.9783E-01, -3.6618E-01, -2.5541E-01, -1.8197E-01, -1.5631E-01
626  },
627  {
628  /*ka9*/ 2.3273E+00, 2.2758E+00, 2.1243E+00, 1.8812E+00, 1.5605E+00, 1.1808E+00, 7.6540E-01,
629  3.4042E-01, -6.6025E-02, -4.2610E-01, -7.1493E-01, -9.1390E-01, -1.0138E+00, -1.0170E+00,
630  -9.3769E-01, -7.9931E-01, -6.3089E-01, -4.6239E-01, -3.2076E-01, -2.2695E-01, -1.9418E-01
631  },
632  {
633  /*ka10*/ 2.7492E+00, 2.6887E+00, 2.5106E+00, 2.2244E+00, 1.8459E+00, 1.3965E+00, 9.0236E-01,
634  3.9399E-01, -9.5749E-02, -5.3340E-01, -8.8806E-01, -1.1353E+00, -1.2618E+00, -1.2684E+00,
635  -1.1715E+00, -1.0001E+00, -7.9063E-01, -5.8088E-01, -4.0461E-01, -2.8794E-01, -2.4721E-01
636  },
637  {
638  /*ka11*/ 3.1894E+00, 3.1197E+00, 2.9141E+00, 2.5834E+00, 2.1451E+00, 1.6229E+00, 1.0464E+00,
639  4.4972E-01, -1.2959E-01, -6.5256E-01, -1.0819E+00, -1.3865E+00, -1.5474E+00, -1.5627E+00,
640  -1.4504E+00, -1.2450E+00, -9.9164E-01, -7.3698E-01, -5.2275E-01, -3.8095E-01, -3.3145E-01
641  },
642  {
643  /*ka12*/ 3.6467E+00, 3.5674E+00, 3.3337E+00, 2.9572E+00, 2.4572E+00, 1.8599E+00, 1.1978E+00,
644  5.0838E-01, -1.6650E-01, -7.8278E-01, -1.2969E+00, -1.6703E+00, -1.8769E+00, -1.9100E+00,
645  -1.7875E+00, -1.5497E+00, -1.2509E+00, -9.4848E-01, -6.9332E-01, -5.2427E-01, -4.6524E-01
646  },
647  {
648  /*ka13*/ 4.0966E+00, 4.0077E+00, 3.7455E+00, 3.3227E+00, 2.7605E+00, 2.0874E+00, 1.3388E+00,
649  5.5573E-01, -2.1615E-01, -9.2821E-01, -1.5310E+00, -1.9786E+00, -2.2366E+00, -2.2918E+00,
650  -2.1605E+00, -1.8882E+00, -1.5398E+00, -1.1845E+00, -8.8398E-01, -6.8475E-01, -6.1517E-01
651  },
652  {
653  /*ka14*/ 4.5369E+00, 4.4382E+00, 4.1471E+00, 3.6774E+00, 3.0522E+00, 2.3026E+00, 1.4670E+00,
654  5.8981E-01, -2.7973E-01, -1.0889E+00, -1.7833E+00, -2.3100E+00, -2.6256E+00, -2.7083E+00,
655  -2.5709E+00, -2.2632E+00, -1.8614E+00, -1.4484E+00, -1.0982E+00, -8.6591E-01, -7.8478E-01
656  },
657  {
658  /*ka15*/ 4.9652E+00, 4.8566E+00, 4.5361E+00, 4.0188E+00, 3.3298E+00, 2.5030E+00, 1.5798E+00,
659  6.0833E-01, -3.5883E-01, -1.2655E+00, -2.0530E+00, -2.6626E+00, -3.0418E+00, -3.1586E+00,
660  -3.0192E+00, -2.6764E+00, -2.2180E+00, -1.7428E+00, -1.3386E+00, -1.0703E+00, -9.7666E-01
661  },
662  {
663  /*ka16*/ 5.3802E+00, 5.2615E+00, 4.9111E+00, 4.3454E+00, 3.5916E+00, 2.6865E+00, 1.6752E+00,
664  6.0927E-01, -4.5493E-01, -1.4582E+00, -2.3388E+00, -3.0336E+00, -3.4819E+00, -3.6400E+00,
665  -3.5043E+00, -3.1278E+00, -2.6105E+00, -2.0689E+00, -1.6065E+00, -1.2994E+00, -1.1921E+00
666  },
667  {
668  /*ka17*/ 5.7809E+00, 5.6519E+00, 5.2709E+00, 4.6558E+00, 3.8360E+00, 2.8514E+00, 1.7510E+00,
669  5.9048E-01, -5.6989E-01, -1.6679E+00, -2.6400E+00, -3.4203E+00, -3.9418E+00, -4.1484E+00,
670  -4.0235E+00, -3.6166E+00, -3.0391E+00, -2.4271E+00, -1.9025E+00, -1.5536E+00, -1.4318E+00
671  },
672  {
673  /*ka18*/ 6.1672E+00, 6.0276E+00, 5.6153E+00, 4.9496E+00, 4.0623E+00, 2.9966E+00, 1.8057E+00,
674  5.5010E-01, -7.0558E-01, -1.8960E+00, -2.9562E+00, -3.8199E+00, -4.4165E+00, -4.6786E+00,
675  -4.5731E+00, -4.1408E+00, -3.5029E+00, -2.8170E+00, -2.2262E+00, -1.8328E+00, -1.6954E+00
676  },
677  {
678  /*ka19*/ 6.5395E+00, 6.3890E+00, 5.9446E+00, 5.2269E+00, 4.2703E+00, 3.1214E+00, 1.8382E+00,
679  4.8645E-01, -8.6398E-01, -2.1439E+00, -3.2877E+00, -4.2304E+00, -4.9013E+00, -5.2244E+00,
680  -5.1475E+00, -4.6969E+00, -4.0000E+00, -3.2376E+00, -2.5767E+00, -2.1360E+00, -1.9820E+00
681  },
682  {
683  /*ka20*/ 6.8987E+00, 6.7370E+00, 6.2595E+00, 5.4883E+00, 4.4602E+00, 3.2257E+00, 1.8476E+00,
684  3.9797E-01, -1.0472E+00, -2.4141E+00, -3.6356E+00, -4.6504E+00, -5.3917E+00, -5.7788E+00,
685  -5.7401E+00, -5.2805E+00, -4.5278E+00, -3.6868E+00, -2.9524E+00, -2.4616E+00, -2.2901E+00
686  },
687  {
688  /*ka21*/ 7.2462E+00, 7.0730E+00, 6.5614E+00, 5.7349E+00, 4.6328E+00, 3.3096E+00, 1.8335E+00,
689  2.8322E-01, -1.2576E+00, -2.7091E+00, -4.0024E+00, -5.0802E+00, -5.8841E+00, -6.3348E+00,
690  -6.3433E+00, -5.8859E+00, -5.0828E+00, -4.1626E+00, -3.3514E+00, -2.8078E+00, -2.6178E+00
691  },
692  {
693  /*ka22*/ 7.5836E+00, 7.3985E+00, 6.8516E+00, 5.9679E+00, 4.7890E+00, 3.3735E+00, 1.7953E+00,
694  1.4071E-01, -1.4979E+00, -3.0326E+00, -4.3913E+00, -5.5215E+00, -6.3766E+00, -6.8863E+00,
695  -6.9488E+00, -6.5063E+00, -5.6609E+00, -4.6621E+00, -3.7715E+00, -3.1726E+00, -2.9630E+00
696  },
697  {
698  /*ka23*/ 7.9125E+00, 7.7152E+00, 7.1317E+00, 6.1886E+00, 4.9299E+00, 3.4179E+00, 1.7326E+00,
699  -3.1215E-02, -1.7709E+00, -3.3886E+00, -4.8069E+00, -5.9777E+00, -6.8694E+00, -7.4281E+00,
700  -7.5480E+00, -7.1343E+00, -6.2574E+00, -5.1825E+00, -4.2103E+00, -3.5537E+00, -3.3237E+00
701  },
702  {
703  /*ka24*/ 8.2345E+00, 8.0244E+00, 7.4032E+00, 6.3983E+00, 5.0563E+00, 3.4431E+00, 1.6449E+00,
704  -2.3448E-01, -2.0804E+00, -3.7824E+00, -5.2552E+00, -6.4540E+00, -7.3646E+00, -7.9573E+00,
705  -8.1330E+00, -7.7614E+00, -6.8670E+00, -5.7206E+00, -4.6657E+00, -3.9491E+00, -3.6977E+00
706  },
707  {
708  /*ka25*/ 8.5510E+00, 8.3278E+00, 7.6673E+00, 6.5982E+00, 5.1690E+00, 3.4493E+00, 1.5312E+00,
709  -4.7147E-01, -2.4305E+00, -4.2200E+00, -5.7434E+00, -6.9574E+00, -7.8668E+00, -8.4733E+00,
710  -8.6970E+00, -8.3789E+00, -7.4837E+00, -6.2733E+00, -5.1355E+00, -4.3570E+00, -4.0835E+00
711  },
712  {
713  /*ka26*/ 8.8632E+00, 8.6262E+00, 7.9252E+00, 6.7892E+00, 5.2687E+00, 3.4364E+00, 1.3906E+00,
714  -7.4501E-01, -2.8265E+00, -4.7089E+00, -6.2803E+00, -7.4966E+00, -8.3829E+00, -8.9783E+00,
715  -9.2352E+00, -8.9780E+00, -8.1010E+00, -6.8372E+00, -5.6177E+00, -4.7758E+00, -4.4793E+00
716  },
717  {
718  /*ka27*/ 9.1717E+00, 8.9207E+00, 8.1775E+00, 6.9719E+00, 5.3555E+00, 3.4041E+00, 1.2212E+00,
719  -1.0587E+00, -3.2745E+00, -5.2580E+00, -6.8767E+00, -8.0824E+00, -8.9219E+00, -9.4773E+00,
720  -9.7454E+00, -9.5504E+00, -8.7116E+00, -7.4086E+00, -6.1105E+00, -5.2042E+00, -4.8840E+00
721  },
722  {
723  /*ka28*/ 9.4773E+00, 9.2117E+00, 8.4247E+00, 7.1465E+00, 5.4296E+00, 3.3517E+00, 1.0212E+00,
724  -1.4170E+00, -3.7824E+00, -5.8785E+00, -7.5458E+00, -8.7278E+00, -9.4948E+00, -9.9776E+00,
725  -1.0228E+01, -1.0089E+01, -9.3081E+00, -7.9839E+00, -6.6123E+00, -5.6413E+00, -5.2969E+00
726  },
727  {
728  /*ka29*/ 9.7800E+00, 9.4993E+00, 8.6669E+00, 7.3131E+00, 5.4905E+00, 3.2779E+00, 7.8776E-01,
729  -1.8253E+00, -4.3595E+00, -6.5844E+00, -8.3040E+00, -9.4487E+00, -1.0115E+01, -1.0489E+01,
730  -1.0688E+01, -1.0589E+01, -9.8825E+00, -8.5589E+00, -7.1217E+00, -6.0864E+00, -5.7173E+00
731  },
732  {
733  /*ka30*/ 1.0080E+01, 9.7834E+00, 8.9039E+00, 7.4713E+00, 5.5375E+00, 3.1815E+00, 5.1771E-01,
734  -2.2903E+00, -5.0181E+00, -7.3937E+00, -9.1727E+00, -1.0265E+01, -1.0798E+01, -1.1024E+01,
735  -1.1130E+01, -1.1049E+01, -1.0427E+01, -9.1291E+00, -7.6373E+00, -6.5392E+00, -6.1452E+00
736  },
737  {
738  /*ka31*/ 1.0376E+01, 1.0064E+01, 9.1354E+00, 7.6206E+00, 5.5698E+00, 3.0605E+00, 2.0709E-01,
739  -2.8202E+00, -5.7735E+00, -8.3303E+00, -1.0180E+01, -1.1202E+01, -1.1563E+01, -1.1594E+01,
740  -1.1565E+01, -1.1468E+01, -1.0935E+01, -9.6893E+00, -8.1576E+00, -6.9995E+00, -6.5807E+00
741  },
742  {
743  /*ka32*/ 1.0668E+01, 1.0339E+01, 9.3604E+00, 7.7600E+00, 5.5861E+00, 2.9127E+00, -1.4879E-01,
744  -3.4253E+00, -6.6463E+00, -9.4274E+00, -1.1367E+01, -1.2295E+01, -1.2433E+01, -1.2217E+01,
745  -1.2003E+01, -1.1852E+01, -1.1400E+01, -1.0234E+01, -8.6812E+00, -7.4676E+00, -7.0243E+00
746  },
747  {
748  /*ka33*/ 1.0956E+01, 1.0609E+01, 9.5782E+00, 7.8886E+00, 5.5849E+00, 2.7355E+00, -5.5553E-01,
749  -4.1185E+00, -7.6641E+00, -1.0733E+01, -1.2791E+01, -1.3593E+01, -1.3439E+01, -1.2909E+01,
750  -1.2457E+01, -1.2206E+01, -1.1820E+01, -1.0758E+01, -9.2063E+00, -7.9436E+00, -7.4766E+00
751  },
752  {
753  /*ka34*/ 1.1237E+01, 1.0873E+01, 9.7877E+00, 8.0051E+00, 5.5648E+00, 2.5261E+00, -1.0198E+00,
754  -4.9166E+00, -8.8666E+00, -1.2322E+01, -1.4549E+01, -1.5171E+01, -1.4623E+01, -1.3692E+01,
755  -1.2939E+01, -1.2540E+01, -1.2193E+01, -1.1254E+01, -9.7309E+00, -8.4279E+00, -7.9385E+00
756  },
757  {
758  /*ka35*/ 1.1512E+01, 1.1129E+01, 9.9878E+00, 8.1084E+00, 5.5239E+00, 2.2812E+00, -1.5494E+00,
759  -5.8416E+00, -1.0313E+01, -1.4320E+01, -1.6807E+01, -1.7156E+01, -1.6046E+01, -1.4592E+01,
760  -1.3464E+01, -1.2864E+01, -1.2521E+01, -1.1718E+01, -1.0253E+01, -8.9208E+00, -8.4107E+00
761  },
762  {
763  /*ka36*/ 1.1778E+01, 1.1376E+01, 1.0177E+01, 8.1972E+00, 5.4605E+00, 1.9974E+00, -2.1540E+00,
764  -6.9238E+00, -1.2098E+01, -1.6967E+01, -1.9918E+01, -1.9789E+01, -1.7805E+01, -1.5641E+01,
765  -1.4046E+01, -1.3190E+01, -1.2810E+01, -1.2145E+01, -1.0769E+01, -9.4225E+00, -8.8943E+00
766  },
767  {
768  /*ka37*/ 1.2035E+01, 1.1614E+01, 1.0355E+01, 8.2703E+00, 5.3729E+00, 1.6709E+00, -2.8454E+00,
769  -8.2070E+00, -1.4393E+01, -2.0830E+01, -2.4847E+01, -2.3638E+01, -2.0073E+01, -1.6886E+01,
770  -1.4703E+01, -1.3528E+01, -1.3066E+01, -1.2530E+01, -1.1275E+01, -9.9329E+00, -9.3900E+00
771  },
772  {
773  /*ka38*/ 1.2282E+01, 1.1841E+01, 1.0521E+01, 8.3267E+00, 5.2594E+00, 1.2974E+00, -3.6384E+00,
774  -9.7578E+00, -1.7557E+01, -2.7951E+01, -3.7266E+01, -3.0744E+01, -2.3208E+01, -1.8396E+01,
775  -1.5453E+01, -1.3891E+01, -1.3297E+01, -1.2873E+01, -1.1769E+01, -1.0452E+01, -9.8987E+00
776  }
777  };
778 
779  Spectrum ChimneyFaceDirectivity::lwAdjustment(Vector direction, double distance)
780  {
781  Spectrum directivity_spectrum;
782 
783  Spectrum spectre_ka = atmosphere->get_k().mult(support_size);
784 
785  double theta = direction.angle(support_normal); // Angle du segment par rapport a la normale au support de la source
786 
787  // angle compris entre 0 et pi
788  if (theta < -M_PI) { theta = theta + M_2PI; }
789  if (theta > M_PI) { theta = theta - M_2PI; }
790  theta = ABS(theta); // On prend la valeur absolue de theta
791 
792  int indice_theta = (int)(20 * theta / M_PI); // Indice de l'angle theta dans le tableau
793  for (unsigned int i = 0 ; i < spectre_ka.getNbValues(); i++)
794  {
795  double ka = spectre_ka.getTabValReel()[i];
796  ka = ka > 3.8 ? 3.8 : ka ;
797  int indice_Ka = (int)(10 * ka);
798 
799  directivity_spectrum.getTabValReel()[i] = compute_q(indice_Ka, indice_theta, ka, theta);
800  }
801 
802  return directivity_spectrum;
803  }
804 
805  double ChimneyFaceDirectivity::compute_q(int ka_idx, int theta_idx, double ka, double theta)
806  {
807  /* ===========================================================================================
808  Recherche par interpolation lineaire dans le tableau _tabQ
809  de la valeur de Q la plus proche
810  ===========================================================================================*/
811  int indiceKaTab = ka_idx > 0 ? ka_idx - 1 : 0; // Eviter les depassement de tableau
812  indiceKaTab = indiceKaTab > (NB_KA - 2) ? NB_KA - 2 : indiceKaTab; // Eviter les depassement de tableau
813 
814  double tabQ1_1 = _tabQ[indiceKaTab] [theta_idx]; //_tabQ[indice_Ka] [indice_theta]
815  double tabQ1_2 = _tabQ[indiceKaTab] [theta_idx + 1]; //_tabQ[indice_Ka] [indice_theta+1]
816  double tabQ2_1 = _tabQ[indiceKaTab + 1][theta_idx]; //_tabQ[indice_Ka+1][indice_theta]
817  double tabQ2_2 = _tabQ[indiceKaTab + 1][theta_idx + 1]; //_tabQ[indice_Ka+1][indice_theta+1]
818 
819  double ka1 = ka_idx / 10.0;
820  double ka2 = (ka_idx + 1) / 10.0;
821 
822  double theta1 = theta_idx * M_PI / 20.0;
823  double theta2 = (theta_idx + 1) * M_PI / 20.0;
824 
825  double val1 = tabQ1_1 + (theta - theta1) * (tabQ1_2 - tabQ1_1) / (theta2 - theta1);
826  double val2 = tabQ2_1 + (theta - theta1) * (tabQ2_2 - tabQ2_1) / (theta2 - theta1);
827 
828  double val = val1 + (ka2 - ka) * (val2 - val1) / (ka2 - ka1);
829 
830 
831  return pow(10.0, val / 10.0);
832  }
833  //
834  // --------------------------
835  //
836 #ifdef NB_KA
837 # undef NB_KA
838 # define NB_KA 9
839 #else
840 # define NB_KA 9
841 #endif
842 
843  // nombre de valeurs de theta dans le tableau
844 #ifdef NB_THETA
845 # undef NB_THETA
846 # define NB_THETA 41
847 #else
848 # define NB_THETA 41
849 #endif
850 
851  const double BaffledFaceDirectivity::_tabKa[NB_KA] = { 0.05, 0.1, 0.25, 0.5, 1.0, 2.0, 5.0, 10.0, 20.0 };
852 
853  const double BaffledFaceDirectivity::_tabQ[NB_KA][NB_THETA] =
854  {
855  {
856  /*ka=0,05*/ 2.0382E+00, 2.0382E+00, 1.7752E+00, 1.7752E+00, 2.0382E+00, 2.0382E+00, 1.6953E+00,
857  1.4101E+00, 1.2860E+00, 1.1729E+00, 8.1142E-01, 6.0152E-01, 4.4591E-01, 3.3056E-01,
858  2.4505E-01, 1.8165E-01, 1.3466E-01, 9.9827E-02, 7.4003E-02, 5.4859E-02, 4.0667E-02,
859  5.4859E-02, 7.4003E-02, 9.9827E-02, 1.3466E-01, 1.8165E-01, 2.4505E-01, 3.3056E-01,
860  4.4591E-01, 6.0152E-01, 8.1142E-01, 1.1729E+00, 1.2860E+00, 1.4101E+00, 1.6953E+00,
861  2.0382E+00, 2.0382E+00, 1.7752E+00, 1.7752E+00, 2.0382E+00, 2.0382E+00
862  },
863  {
864  /*ka=0,1*/ 1.0215E+00, 1.0215E+00, 1.0215E+00, 1.1201E+00, 1.2860E+00, 1.2860E+00, 1.2860E+00,
865  1.4765E+00, 1.6190E+00, 1.6190E+00, 1.6190E+00, 1.1201E+00, 7.7490E-01, 5.3610E-01,
866  3.7089E-01, 2.5659E-01, 1.7752E-01, 1.2281E-01, 8.4966E-02, 5.8782E-02, 4.0667E-02,
867  5.8782E-02, 8.4966E-02, 1.2281E-01, 1.7752E-01, 2.5659E-01, 3.7089E-01, 5.3610E-01,
868  7.7490E-01, 1.1201E+00, 1.6190E+00, 1.6190E+00, 1.6190E+00, 1.4765E+00, 1.2860E+00,
869  1.2860E+00, 1.2860E+00, 1.1201E+00, 1.0215E+00, 1.0215E+00, 1.0215E+00
870  },
871  {
872  /*ka=0,25*/ 2.1534E+00, 1.7911E+00, 1.4898E+00, 1.3587E+00, 1.4227E+00, 1.7105E+00, 1.7105E+00,
873  1.4898E+00, 1.3587E+00, 1.3587E+00, 1.3587E+00, 9.6189E-01, 6.8096E-01, 4.8209E-01,
874  3.4129E-01, 2.4162E-01, 1.7105E-01, 1.2109E-01, 8.5728E-02, 6.0691E-02, 4.2966E-02,
875  6.0691E-02, 8.5728E-02, 1.2109E-01, 1.7105E-01, 2.4162E-01, 3.4129E-01, 4.8209E-01,
876  6.8096E-01, 9.6189E-01, 1.3587E+00, 1.3587E+00, 1.3587E+00, 1.4898E+00, 1.7105E+00,
877  1.7105E+00, 1.4227E+00, 1.3587E+00, 1.4898E+00, 1.7911E+00, 2.1534E+00
878  },
879  {
880  /*ka=0,5*/ 1.4385E+00, 1.4385E+00, 1.4385E+00, 1.5773E+00, 1.8109E+00, 1.8109E+00, 2.1772E+00,
881  1.9857E+00, 1.5063E+00, 1.0912E+00, 9.0762E-01, 6.7283E-01, 4.9877E-01, 3.6975E-01,
882  2.7410E-01, 2.0319E-01, 1.5063E-01, 1.1166E-01, 8.2776E-02, 6.1363E-02, 4.5489E-02,
883  6.1363E-02, 8.2776E-02, 1.1166E-01, 1.5063E-01, 2.0319E-01, 2.7410E-01, 3.6975E-01,
884  4.9877E-01, 6.7283E-01, 9.0762E-01, 1.0912E+00, 1.5063E+00, 1.9857E+00, 2.1772E+00,
885  1.8109E+00, 1.8109E+00, 1.5773E+00, 1.4385E+00, 1.4385E+00, 1.4385E+00
886  },
887  {
888  /*ka=1*/ 3.5678E+00, 4.2895E+00, 4.4916E+00, 3.7360E+00, 2.7065E+00, 2.2511E+00, 1.2954E+00,
889  8.5586E-01, 7.8055E-01, 8.5586E-01, 7.1187E-01, 5.4001E-01, 4.0964E-01, 3.1074E-01,
890  2.3572E-01, 1.7881E-01, 1.3564E-01, 1.0290E-01, 7.8055E-02, 5.9211E-02, 4.4916E-02,
891  5.9211E-02, 7.8055E-02, 1.0290E-01, 1.3564E-01, 1.7881E-01, 2.3572E-01, 3.1074E-01,
892  4.0964E-01, 5.4001E-01, 7.1187E-01, 8.5586E-01, 7.8055E-01, 8.5586E-01, 1.2954E+00,
893  2.2511E+00, 2.7065E+00, 3.7360E+00, 4.4916E+00, 4.2895E+00, 3.5678E+00
894  },
895  {
896  /*ka=2*/ 1.1891E+01, 9.8907E+00, 7.1652E+00, 4.5209E+00, 2.6015E+00, 1.4970E+00, 1.2452E+00,
897  7.8565E-01, 5.4353E-01, 4.5209E-01, 3.7603E-01, 3.0565E-01, 2.4844E-01, 2.0194E-01,
898  1.6414E-01, 1.3342E-01, 1.0845E-01, 8.8151E-02, 7.1652E-02, 5.8241E-02, 4.7340E-02,
899  5.8241E-02, 7.1652E-02, 8.8151E-02, 1.0845E-01, 1.3342E-01, 1.6414E-01, 2.0194E-01,
900  2.4844E-01, 3.0565E-01, 3.7603E-01, 4.5209E-01, 5.4353E-01, 7.8565E-01, 1.2452E+00,
901  1.4970E+00, 2.6015E+00, 4.5209E+00, 7.1652E+00, 9.8907E+00, 1.1891E+01
902  },
903  {
904  /*ka=5*/ 4.5941E+00, 4.5941E+00, 4.5941E+00, 4.1899E+00, 3.1784E+00, 1.8290E+00, 1.2653E+00,
905  1.0051E+00, 8.3600E-01, 6.9535E-01, 5.7837E-01, 4.4896E-01, 3.4850E-01, 2.7052E-01,
906  2.0999E-01, 1.6301E-01, 1.2653E-01, 9.8221E-02, 7.6244E-02, 5.9184E-02, 4.5941E-02,
907  5.9184E-02, 7.6244E-02, 9.8221E-02, 1.2653E-01, 1.6301E-01, 2.0999E-01, 2.7052E-01,
908  3.4850E-01, 4.4896E-01, 5.7837E-01, 6.9535E-01, 8.3600E-01, 1.0051E+00, 1.2653E+00,
909  1.8290E+00, 3.1784E+00, 4.1899E+00, 4.5941E+00, 4.5941E+00, 4.5941E+00
910  },
911  {
912  /*ka=10*/ 3.8007E+00, 3.1613E+00, 3.0190E+00, 2.7534E+00, 2.3981E+00, 2.3981E+00, 1.9946E+00,
913  1.4450E+00, 1.0961E+00, 8.7070E-01, 6.0237E-01, 4.6759E-01, 3.6297E-01, 2.8175E-01,
914  2.1871E-01, 1.6977E-01, 1.3179E-01, 1.0230E-01, 7.9408E-02, 6.1641E-02, 4.7848E-02,
915  6.1641E-02, 7.9408E-02, 1.0230E-01, 1.3179E-01, 1.6977E-01, 2.1871E-01, 2.8175E-01,
916  3.6297E-01, 4.6759E-01, 6.0237E-01, 8.7070E-01, 1.0961E+00, 1.4450E+00, 1.9946E+00,
917  2.3981E+00, 2.3981E+00, 2.7534E+00, 3.0190E+00, 3.1613E+00, 3.8007E+00
918  },
919  {
920  /*ka=20*/ 2.3330E+00, 2.3330E+00, 2.3330E+00, 2.3330E+00, 2.3330E+00, 2.3330E+00, 2.3330E+00,
921  1.7698E+00, 1.2244E+00, 8.4708E-01, 5.8603E-01, 4.5491E-01, 3.5312E-01, 2.7411E-01,
922  2.1278E-01, 1.6517E-01, 1.2821E-01, 9.9523E-02, 7.7254E-02, 5.9968E-02, 4.6550E-02,
923  5.9968E-02, 7.7254E-02, 9.9523E-02, 1.2821E-01, 1.6517E-01, 2.1278E-01, 2.7411E-01,
924  3.5312E-01, 4.5491E-01, 5.8603E-01, 8.4708E-01, 1.2244E+00, 1.7698E+00, 2.3330E+00,
925  2.3330E+00, 2.3330E+00, 2.3330E+00, 2.3330E+00, 2.3330E+00, 2.3330E+00
926  }
927  };
928 
929  Spectrum BaffledFaceDirectivity::lwAdjustment(Vector direction, double distance)
930  {
931  Spectrum directivity_spectrum;
932  Spectrum spectre_Ka = atmosphere->get_k().mult(support_size); // 1/2 longueur de la diagonale de la bouche
933 
934  double theta = direction.angle(support_normal); // Angle du segment par rapport a la normale au support de la source
935 
936  // Angle compris entre 0 et 2pi;
937  if (theta < 0) { theta = theta + M_2PI; }
938  if (theta > M_2PI) { theta = theta - M_2PI; }
939 
940  int indice_theta = (int)(20.0 * theta / M_PI);
941 
942  indice_theta = indice_theta > (NB_THETA - 2) ? NB_THETA - 2 : indice_theta; // Eviter les depassement de tableau
943 
944  for (unsigned int i = 0 ; i < directivity_spectrum.getNbValues(); i++)
945  {
946  double ka = spectre_Ka.getTabValReel()[i];
947  ka = ka > 20.0 ? 20.0 : ka ;
948 
949  int indice_Ka = find_Ka_idx(ka);
950 
951  directivity_spectrum.getTabValReel()[i] = compute_q(indice_Ka, indice_theta, ka, theta);
952  }
953 
954  return directivity_spectrum;
955  }
956 
957  double BaffledFaceDirectivity::compute_q(int indice_Ka, int indice_theta, double ka, double theta)
958  {
959  double tabQ1_1 = _tabQ[indice_Ka] [indice_theta];
960  double tabQ1_2 = _tabQ[indice_Ka] [indice_theta + 1];
961  double tabQ2_1 = _tabQ[indice_Ka + 1][indice_theta];
962  double tabQ2_2 = _tabQ[indice_Ka + 1][indice_theta + 1];
963 
964  double ka1 = _tabKa[indice_Ka];
965  double ka2 = _tabKa[indice_Ka + 1];
966 
967  double theta1 = indice_theta * M_PI / 20;
968  double theta2 = (indice_theta + 1) * M_PI / 20;
969 
970  double val1 = tabQ1_1 + (theta - theta1) * (tabQ1_2 - tabQ1_1) / (theta2 - theta1);
971  double val2 = tabQ2_1 + (theta - theta1) * (tabQ2_2 - tabQ2_1) / (theta2 - theta1);
972 
973  return val1 + (ka - ka1) * (val2 - val1) / (ka2 - ka1);
974  }
975 
976  int BaffledFaceDirectivity::find_Ka_idx(double ka)
977  {
978  int indice = 0;
979 
980  while ((_tabKa[indice] < ka) && (indice < (NB_KA - 1))) { indice++; }
981 
982  return indice > (NB_KA - 2) ? NB_KA - 2 : indice - 1; // Eviter les depassement de tableau
983  // return indice > 0 ? indice - 1 : 0;
984  }
985  //
986  // --------------------------
987  //
988  AcousticSource::AcousticSource(
989  const Point& position_,
990  const Spectrum& spectrum_ ,
991  SourceDirectivityInterface* directivity_)
992  : position(position_)
993  , spectrum(spectrum_)
994  , directivity(directivity_)
995  {}
996 
997  // ---------
998 
999  AcousticReceptor::AcousticReceptor(
1000  const Point& position_
1001  ) : position(position_) {}
1002 
1003 
1004  Point ComputeAcousticCentroid(const source_pool_t &tabSources_)
1005  {
1006  Point centroid;
1007  double cumul_x=0., cumul_y=0., cumul_z=0., cumul_level = 0.;
1008  std::deque<double> tabLevels;
1009 
1010  for(unsigned int i=0; i<tabSources_.size(); i++)
1011  {
1012  tabLevels.push_back( ::pow( 10, tabSources_[i].spectrum.valGlobDBA()/10 ) );
1013  }
1014 
1015  for(unsigned int i=0; i<tabSources_.size(); i++)
1016  {
1017  cumul_x = cumul_x + tabSources_[i].position._x * tabLevels[i];
1018  cumul_y = cumul_y + tabSources_[i].position._y * tabLevels[i];
1019  cumul_z = cumul_z + tabSources_[i].position._z * tabLevels[i];
1020  cumul_level = cumul_level + tabLevels[i];
1021  }
1022 
1023  centroid._x = cumul_x / cumul_level;
1024  centroid._y = cumul_y / cumul_level;
1025  centroid._z = cumul_z / cumul_level;
1026 
1027  return centroid;
1028  }
1029 } /* namespace tympan */
tympan::BaffledFaceDirectivity::compute_q
double compute_q(int indice_Ka, int indice_theta, double ka, double theta)
Definition: entities.cpp:957
entities.hpp
This file provides the declaration of the entities of the model, which inherit from BaseEntity...
CPLX_J
#define CPLX_J
Definition: macros.h:30
tympan::AcousticGroundMaterial::erfc_G_computation
void erfc_G_computation(const ComplexSpectrum &localW, ComplexSpectrum &G)
Definition: entities.cpp:308
MIN
#define MIN(A, B)
Definition: 3d.cpp:25
OVector3D
The 3D vector class.
Definition: 3d.h:296
tympan::AcousticGroundMaterial::computeQ
void computeQ(double angle, ComplexSpectrum &Rp, ComplexSpectrum &Fw, ComplexSpectrum &Q)
Compute reflection coefficient.
Definition: entities.cpp:381
tympan::VolumeFaceDirectivity::_tabRA
static const double _tabRA[]
RA form factor.
Definition: entities.hpp:230
name
const char * name
Definition: TYRouteWidget.cpp:33
cotanh
TYComplex cotanh(const TYComplex &valeur)
Definition: macros.h:35
OSpectre::sum
virtual OSpectre sum(const OSpectre &spectre) const
Arithmetic sum of two spectrums in one-third Octave.
Definition: spectre.cpp:314
tympan::AcousticSource::position
Point position
Destructor.
Definition: entities.hpp:325
tympan::AcousticGroundMaterial::limit_W_values
void limit_W_values(ComplexSpectrum &localW)
Definition: entities.cpp:289
AtmosphericConditions::get_k
const OSpectre & get_k() const
Get the wave number spectrum.
Definition: atmospheric_conditions.h:39
tympan::ChimneyFaceDirectivity::compute_q
double compute_q(int ka_idx, int theta_idx, double ka, double theta)
Definition: entities.cpp:805
tympan::VolumeFaceDirectivity::_tabCor
static const double _tabCor[]
Correction factors.
Definition: entities.hpp:231
tympan::AcousticGroundMaterial::computeZs
void computeZs(double angle, ComplexSpectrum Z, ComplexSpectrum &spectrum)
Compute specific impedance.
Definition: entities.cpp:107
tympan::AcousticGroundMaterial::erfcCas1
TYComplex erfcCas1(const TYComplex &wValue) const
: Functions used in Fw computation
Definition: entities.cpp:397
tympan::AcousticBuildingMaterial::AcousticBuildingMaterial
AcousticBuildingMaterial(const string &name_, const ComplexSpectrum &spectrum)
Constructor.
Definition: entities.cpp:23
tympan::AcousticGroundMaterial::sgn_G_computation
void sgn_G_computation(const ComplexSpectrum &localW, ComplexSpectrum &G)
Definition: entities.cpp:337
tympan::AcousticGroundMaterial::K
ComplexSpectrum K
Wave number.
Definition: entities.hpp:135
tympan::AcousticGroundMaterial::erfcCas3
TYComplex erfcCas3(const TYComplex &wValue) const
Definition: entities.cpp:478
tympan::AcousticGroundMaterial::erfcCas2
TYComplex erfcCas2(const TYComplex &wValue) const
Definition: entities.cpp:466
tympan::ChimneyFaceDirectivity::_tabQ
static const double _tabQ[NB_KA][NB_THETA]
Definition: entities.hpp:268
atmospheric_conditions.h
tympan::AcousticGroundMaterial::trapz
double trapz(std::vector< double > u, std::vector< double > integrande)
Definition: entities.cpp:214
tympan::AcousticGroundMaterial::computeRp
void computeRp(double angle, const ComplexSpectrum &Zs, ComplexSpectrum &Rp)
Compute reflection coefficient for plane waves.
Definition: entities.cpp:222
TYComplex
std::complex< double > TYComplex
Definition: macros.h:26
NB_KA
#define NB_KA
Definition: entities.cpp:838
OSpectre::getTabFreqExact
static OTabFreq getTabFreqExact()
Definition: spectre.cpp:644
OSpectreComplex::getTabValReel
virtual double * getTabValReel()
Get an array of the real values of the spectrum.
Definition: spectre.h:378
OSpectre::getNbValues
virtual unsigned int getNbValues() const
Number of values in the spectrum.
Definition: spectre.cpp:789
OSpectre::setType
virtual void setType(TYSpectreType type)
Set the spectrum type.
Definition: spectre.h:101
CPLX_MUN
#define CPLX_MUN
Definition: macros.h:29
tympan::AcousticGroundMaterial::computeFw
void computeFw(ComplexSpectrum localW, ComplexSpectrum &Fw)
Compute function of numeric distance.
Definition: entities.cpp:261
tympan::AcousticGroundMaterial::sgnReIm
TYComplex sgnReIm(const TYComplex &W, const TYComplex &G) const
: function used in G computation
Definition: entities.cpp:490
OSpectreComplex::toModulePhase
OSpectreComplex toModulePhase() const
Conversion to module/phase.
Definition: spectre.cpp:1149
tympan::AcousticBuildingMaterial::asEyring
ComplexSpectrum asEyring() const
Returns a spectrum with Eyring formulae.
Definition: entities.cpp:28
OCoord3D::_y
double _y
y coordinate of OCoord3D
Definition: 3d.h:281
OCoord3D::_x
double _x
x coordinate of OCoord3D
Definition: 3d.h:280
tympan::AcousticReceptor::AcousticReceptor
AcousticReceptor(const Point &position_)
Constructor to build a receptor defined by the its position.
Definition: entities.cpp:999
tympan::AcousticGroundMaterial::computeZc
void computeZc()
Compute characteristic impedance.
Definition: entities.cpp:72
tympan::AcousticMaterialBase::AcousticMaterialBase
AcousticMaterialBase(const string &name_)
Constructor.
Definition: entities.cpp:16
CPLX_UN
#define CPLX_UN
Definition: macros.h:28
tympan::AcousticGroundMaterial::computeK
void computeK()
Compute wave number.
Definition: entities.cpp:88
tympan::AcousticGroundMaterial::atmosphere
static AtmosphericConditions * atmosphere
Pointer to current atmosphere.
Definition: entities.hpp:105
tympan::BaffledFaceDirectivity::_tabQ
static const double _tabQ[NB_KA][NB_THETA]
Definition: entities.hpp:304
tympan::CommonFaceDirectivity::atmosphere
static AtmosphericConditions * atmosphere
Characteristic size of support face.
Definition: entities.hpp:212
tympan::AcousticSource::AcousticSource
AcousticSource(const Point &point_, const Spectrum &spectrum_, SourceDirectivityInterface *directivity)
Constructor to build a source defined by a point, spectrum, directivity.
Definition: entities.cpp:988
tympan::AcousticGroundMaterial::Zc
ComplexSpectrum Zc
Characteristic impedance.
Definition: entities.hpp:134
tympan::VolumeFaceDirectivity::lwAdjustment
virtual Spectrum lwAdjustment(Vector direction, double distance)
Directivity of a volume face.
Definition: entities.cpp:548
OSpectreComplex::getTabValImag
double * getTabValImag()
Get an array of the imaginary values of the spectrum.
Definition: spectre.h:381
tympan::BaffledFaceDirectivity::lwAdjustment
virtual Spectrum lwAdjustment(Vector direction, double distance)
Directivity of a baffled face source.
Definition: entities.cpp:929
tympan::AcousticMaterialBase
Base class for material.
Definition: entities.hpp:22
OSpectre::getTabValReel
virtual double * getTabValReel()
Definition: spectre.h:112
M_2PI
#define M_2PI
2Pi.
Definition: 3d.h:57
OTabFreq
std::vector< double > OTabFreq
Frequencies collection.
Definition: spectre.h:40
OPoint3D
The 3D point class.
Definition: 3d.h:484
tympan::source_pool_t
std::deque< AcousticSource > source_pool_t
Array of sources.
Definition: entities.hpp:334
OSpectre::mult
virtual OSpectre mult(const OSpectre &spectre) const
Multiplication of two spectrums.
Definition: spectre.cpp:383
AtmosphericConditions
Class for the definition of atmospheric conditions.
Definition: atmospheric_conditions.h:12
tympan::AcousticBuildingMaterial::spectrum
ComplexSpectrum spectrum
Spectrum to store acoustic values at different frequencies.
Definition: entities.hpp:60
tympan::ChimneyFaceDirectivity::lwAdjustment
virtual Spectrum lwAdjustment(Vector direction, double distance)
Directivity of a chimney face source.
Definition: entities.cpp:779
OSpectre::exp
virtual OSpectre exp(const double coef=1.0)
Compute e^(coef * spectre) of this spectrum.
Definition: spectre.cpp:537
OCoord3D::_z
double _z
z coordinate of OCoord3D
Definition: 3d.h:282
OVector3D::angle
double angle(const OVector3D &vector) const
Computes the angle between this vector and another vector.
Definition: 3d.cpp:265
OSpectre
Store acoustic power values for different frequencies.
Definition: spectre.h:49
tympan::ComputeAcousticCentroid
Point ComputeAcousticCentroid(const source_pool_t &tabSources_)
Definition: entities.cpp:1004
ABS
#define ABS(x)
Definition: mathlib.h:84
tympan::AcousticGroundMaterial::gaussianSpectrum
double gaussianSpectrum(double const k, double const sigma, double const lc)
Definition: entities.cpp:209
M_PI
#define M_PI
Pi.
Definition: color.cpp:24
tympan::AcousticGroundMaterial::get_absorption
virtual ComplexSpectrum get_absorption(double incidence_angle, double length)
Destructor.
Definition: entities.cpp:50
tympan::AcousticGroundMaterial::computeW
void computeW(double angle, double length, const ComplexSpectrum &Zs, ComplexSpectrum &W)
Compute numeric distance.
Definition: entities.cpp:239
NB_THETA
#define NB_THETA
Definition: entities.cpp:846
tympan::AcousticGroundMaterial::Zf
ComplexSpectrum Zf
Effective impedance.
Definition: entities.hpp:136
tympan::BaffledFaceDirectivity::_tabKa
static const double _tabKa[NB_KA]
Definition: entities.hpp:305
tympan::SourceDirectivityInterface
Interface for source directivity classes (SphericalSourceDirectivity, CommonFaceDirectivity, ...)
Definition: entities.hpp:164
tympan::AcousticGroundMaterial::get_ISO9613_G
double get_ISO9613_G()
Absorption given by ISO9613.
Definition: entities.cpp:67
tympan::AcousticGroundMaterial::computeZf
void computeZf(double angle, ComplexSpectrum Zs)
Compute effective impedance in rough ground.
Definition: entities.cpp:130
tympan::VolumeFaceDirectivity::calculC
double calculC(double distance)
Compute directivity factor.
Definition: entities.cpp:534
tympan::AcousticGroundMaterial::AcousticGroundMaterial
AcousticGroundMaterial(const string &name_, double resistivity_, double deviation_, double length_)
Constructor.
Definition: entities.cpp:37
Generated on Tue May 25 2021 08:06:51 for Code_TYMPAN by   doxygen 1.8.13