NMPB08 Shared Libraries
|
00001 00008 #include "CalculPropagation.h" 00009 #include "SousCalculs/Diffraction.h" 00010 #include <math.h> 00011 #include <stdio.h> 00012 #include <assert.h> 00013 00014 #include "../test_mem/safe_new.h" 00015 00016 using namespace std; 00017 using namespace DiffractionNMPB; 00018 00019 namespace CalculPropagationNMPB 00020 { 00021 00022 // ------------------------------ utility function ------------------------------------ 00027 static int modified_floor (double index) 00028 { 00029 double x = floor (index) ; 00030 int i = (int) x ; 00031 if (x != index) i++ ; 00032 return i ; 00033 } 00034 00035 // ------------------------------ Sound levels ----------------------------------------- 00047 double SumLevels(int n, double const* levels) 00048 { 00049 double soundsSum = 0; 00050 for(int i = 0 ; i < n ; i++) 00051 { 00052 soundsSum += pow(10, 0.1 * levels[i]); 00053 } 00054 00055 double soundsSummed = 10 * log10(soundsSum); 00056 00057 return soundsSummed; 00058 } 00059 00075 double SoundLevelForPath(double soundLevel_h, double soundLevel_f, double favourableProbability) 00076 { 00077 if(favourableProbability > 1) 00078 { 00079 printf (".ERROR: probability %4.2f > 1 \n", favourableProbability) ; 00080 throw new int( ERRProbability ); 00081 } 00082 double soundLevel = 10 * log10(favourableProbability * pow(10.0, (double)(soundLevel_f/10.0)) + (1 - favourableProbability) * pow(10., (double)soundLevel_h/10.0)); 00083 return soundLevel; 00084 } 00085 00103 double GetFavorableConditionProbability(Position3D const* sourcePos, Position3D const* receiverPos, int nbAngles, double const* fcpAngles, double angleNorth) 00104 { 00105 // angle of source-receiver line relative to X-axis, positive = counter clockwise 00106 double alphaS = 0; 00107 if (sourcePos->x == receiverPos->x && sourcePos->y == receiverPos->y) 00108 { 00109 alphaS = 0; 00110 } 00111 else 00112 { 00113 alphaS = (180 / PI) * atan2(sourcePos->y - receiverPos->y, sourcePos->x - receiverPos->x); 00114 } 00115 00116 // angle relative to north, positive = clockwise 00117 double alphaNS = angleNorth - alphaS; 00118 while (alphaNS <= 10) 00119 { 00120 alphaNS += 360; 00121 } 00122 while (alphaNS > 370) 00123 { 00124 alphaNS -= 360; 00125 } 00126 00127 // index into array of frequencies of occurrence 00128 00129 assert (alphaNS > 10) ; 00130 assert (alphaNS <= 370) ; 00131 double indexR = (alphaNS - 10) / 20 ; 00132 assert (indexR > 0) ; 00133 assert (indexR <= 18) ; 00134 int indexA = modified_floor(indexR) ; 00135 00136 if (indexA > 0 && indexA <= nbAngles) 00137 { 00138 // convert probability ([0,100] --> [0,1]) : 00139 // note that the first value in the table has index 1, i.e. angular sector ]10,30°] 00140 00141 double fp = fcpAngles[indexA - 1] / 100.0; 00142 return fp; 00143 } 00144 else 00145 { 00146 printf (".ERROR: angle %4.1f not in array containing favorable conditions probabilities for the angles \n", alphaNS) ; 00147 throw new int( ERRAngle ); 00148 } 00149 } 00150 00170 void CalculateLeqLT(int nbFreq, double const* Lw, double const* attH, double const* attF, double fcp, double* LeqH, double* LeqF, double* LeqLT) 00171 { 00172 for (int i = 0 ; i < nbFreq ; i++) 00173 { 00174 //LAi,F p.16 § 5.2.1 formula (3) 00175 LeqH[i] = Lw[i] - attH[i]; 00176 //LAi,H p.17 § 5.2.2 formula (5) 00177 LeqF[i] = Lw[i] - attF[i]; 00178 00179 //Li,LT 00180 LeqLT[i] = SoundLevelForPath(LeqH[i], LeqF[i], fcp); 00181 } 00182 } 00183 00184 00185 // ------------------------------ Attenuation Class ----------------------------------------- 00186 00194 Attenuation::Attenuation(PropagationPath* path, double frequency) 00195 { 00196 _path = path; 00197 _frequency = frequency; 00198 _favourableProbability = 0; 00199 _displayConsoleStep = path->GetOption(TRACE_DETAILS); 00200 _withCh = !path->GetOption (FORCE_CH_EQUAL_ONE) ; 00201 ClearResults(); 00202 00203 // fill frequency map with default values : 00204 FillFrequencyMap(); 00205 } 00206 00216 Attenuation::Attenuation(PropagationPath* path, double frequency, double favourableProbability) 00217 { 00218 _path = path; 00219 _frequency = frequency; 00220 _favourableProbability = favourableProbability; 00221 _displayConsoleStep = path->GetOption(TRACE_DETAILS); 00222 _withCh = !path->GetOption (FORCE_CH_EQUAL_ONE) ; 00223 00224 ClearResults(); 00225 00226 // fill frequency map with default values : 00227 FillFrequencyMap(); 00228 } 00229 00240 double Attenuation::DivergenceAttenuationCalculation(double dist) 00241 { 00242 double aDiv = 20 * log10(dist) + 11; 00243 00244 _divergenceAttenuation = aDiv; 00245 return aDiv; 00246 } 00247 00253 void Attenuation::FillFrequencyMap() 00254 { 00255 _attenuationCoeffByFrequencyMap.clear(); 00256 int frequencyNumber = 24; 00257 // Nominal median frequencies array 00258 int NominalMedianFrequency[] = {50, 63, 80, 100, 125, 160, 200, 250, 315, 400, 500, 630, 800, 1000, 1250, 1600, 2000, 2500, 3150, 4000, 5000, 6300, 8000, 10000}; 00259 // Atmospheric attenuation coefficients array 00260 double AtmosphericAttenuationCoeff[] = {0.07, 0.1, 0.17, 0.25, 0.38, 0.57, 0.82, 1.13, 1.51, 1.92, 2.36, 2.84, 3.38, 4.08, 5.05, 6.51, 8.75, 12.2, 17.7, 26.4, 39.9, 61.1, 93.7, 143.5}; 00261 for(int i = 0; i < frequencyNumber ; i++) 00262 { 00263 _attenuationCoeffByFrequencyMap[NominalMedianFrequency[i]] = AtmosphericAttenuationCoeff[i]; 00264 } 00265 } 00266 00272 void Attenuation::FillFrequencyMap(map<int,double> newMap) 00273 { 00274 _attenuationCoeffByFrequencyMap.clear(); 00275 for(map<int,double>::iterator it = newMap.begin(); it != newMap.end(); ++it) 00276 { 00277 _attenuationCoeffByFrequencyMap[it->first] = it->second; 00278 } 00279 00280 00281 00282 } 00283 00284 00285 00297 double Attenuation::AttenuationCoeffByFrequency(double freq) 00298 { 00299 double alpha = 0; 00300 if (_attenuationCoeffByFrequencyMap.find(static_cast<int>(freq)) != _attenuationCoeffByFrequencyMap.end()) 00301 { 00302 alpha = _attenuationCoeffByFrequencyMap[static_cast<int>(freq)]; 00303 } 00304 else 00305 { 00306 printf (".ERROR: frequency %4.1f not found to get attenuation coefficient \n", freq) ; 00307 throw new int( ERRAttCoeffFrequency ); 00308 } 00309 00310 return alpha; 00311 } 00312 00325 double Attenuation::AtmosphericAbsorptionCalculation(double dist, double freq) 00326 { 00327 // alpha : atmospheric attenuation coefficient 00328 double alpha = AttenuationCoeffByFrequency(freq); 00329 00330 double aAtm = alpha * dist / 1000.0; 00331 _atmosphericAbsorption = aAtm; 00332 return aAtm; 00333 } 00334 00335 00343 void Attenuation::GroundAttenuationCalculation() 00344 { 00345 int nbPoints = (int) _path->pathPoints.size(); 00346 if (nbPoints >= 2) 00347 { 00348 MeanPlane meanPlane(_path->pathPoints); 00349 00350 GroundEffect groundEffect(meanPlane.Get_dp(), meanPlane.Get_zEqSource(), meanPlane.Get_zEqReceiver(), meanPlane.Get_Gpath(), meanPlane.Get_Gsource(), _frequency, Asol); 00351 _groundAttenuation_h = groundEffect.AttenuationCalculationH(); 00352 _groundAttenuation_f = groundEffect.AttenuationCalculationF(); 00353 00354 if (_displayConsoleStep == 1) 00355 { 00356 printf(" - AttsolH = %5.1f dB, AttsolF = %5.1f dB \n", _groundAttenuation_h, _groundAttenuation_f); 00357 printf(" Mean plane : z = %.3f d + %.3f ; dp = %6.2f \n", meanPlane.Get_aCoeff(), meanPlane.Get_bCoeff(), meanPlane.Get_dp()); 00358 printf(" h : zEqS = %6.2f, zEqR = %6.2f \n", groundEffect.Get_zEqSource(), groundEffect.Get_zEqReceiver()); 00359 printf(" f : zEqS_f = %6.2f, zEqR_f = %6.2f \n", groundEffect.Get_zEqSource_f(), groundEffect.Get_zEqReceiver_f()); 00360 printf(" Gpath = %6.2f, G'path = %6.2f \n", meanPlane.Get_Gpath(), groundEffect.get_EquivalentGpath()); 00361 } 00362 } 00363 } 00364 00371 void Attenuation::EmbankmentAttenuationCalculation() 00372 { 00373 _embankmentAttenuation_h = 0; 00374 _embankmentAttenuation_f = 0; 00375 if(_path->GetOption(CHECK_EMBANKMENT)) 00376 { 00377 int nbPoints = (int) _path->pathPoints.size(); 00378 if (nbPoints >= 3) 00379 { 00380 Embankment embankment(_path->pathPoints, _frequency); 00381 double embankmentAttenuation = embankment.GetEmbankmentAttenuation(); 00382 _embankmentAttenuation_h = embankmentAttenuation; 00383 _embankmentAttenuation_f = embankmentAttenuation; 00384 00385 if (_displayConsoleStep == 1) 00386 { 00387 printf(" - Attenuation talus = %5.1f dB \n", embankmentAttenuation); 00388 } 00389 } 00390 } 00391 } 00392 00404 void Attenuation::BoundaryAttenuationCalculation(int nbSideDiffractions) 00405 { 00406 int nbPoints = (int) _path->pathPoints.size(); 00407 if (nbPoints < 1) 00408 { 00409 printf (".ERROR: path empty \n") ; 00410 throw new int( ERRNoPoint ); 00411 } 00412 else if (nbPoints < 2) 00413 { 00414 printf (".ERROR: only one point in the path \n") ; 00415 throw new int( ERROnePoint ); 00416 } 00417 00418 // Calculates diffraction or reflection if necessary: 00419 if (nbPoints >= 3) 00420 { 00421 int fPos = _path->GetFrequencyPosition(_frequency); 00422 if(fPos < 0) 00423 { 00424 printf (".ERROR: frequency %4.1f not found in path spectrum \n", _frequency) ; 00425 throw new int( ERRAttCoeffFrequency ); 00426 } 00427 00428 if(nbSideDiffractions > 0) 00429 { 00430 // Side diffraction case 00431 // no embankment 00432 _embankmentAttenuation_h = 0; 00433 _embankmentAttenuation_f = 0; 00434 // ground attenuation is calculated here 00435 GroundAttenuationCalculation(); 00436 00437 // Side diffractions and reflections : 00438 Diffraction diffractionAttenuation(_path->pathPoints, _frequency, fPos, nbSideDiffractions); 00439 _diffractionAttenuation_h = diffractionAttenuation.Get_DeltaDifSR_h() + diffractionAttenuation.Get_absorptionAttenuation(); 00440 _diffractionAttenuation_f = diffractionAttenuation.Get_DeltaDifSR_h() + diffractionAttenuation.Get_absorptionAttenuation(); 00441 } 00442 else 00443 { 00444 // No side diffraction 00445 // ground attenuation calculated in diffraction attenuation, so equals 0 00446 _groundAttenuation_h = 0; 00447 _groundAttenuation_f = 0; 00448 00449 // embankment attenuation calculated in diffraction attenuation, so equals 0 00450 _embankmentAttenuation_h = 0; 00451 _embankmentAttenuation_f = 0; 00452 00453 // diffraction attenuation : 00454 Diffraction diffractionAttenuation(_path->pathPoints, _frequency, fPos, _path->GetOption(CHECK_EMBANKMENT)); 00455 if (diffractionAttenuation.Get_pathDifferenceSR_h() == 0 && diffractionAttenuation.Get_pathDifferenceSR_f() == 0) 00456 { 00457 // if diffraction is negligible, ground attenuation is calculated : 00458 _diffractionAttenuation_h = 0; 00459 _diffractionAttenuation_f = 0; 00460 GroundAttenuationCalculation(); 00461 EmbankmentAttenuationCalculation(); 00462 } 00463 else 00464 { 00465 diffractionAttenuation.CalculAttenuationDiffraction (_withCh); 00466 _diffractionAttenuation_h = diffractionAttenuation.Get_diffractionAttenuation_h(); 00467 _diffractionAttenuation_f = diffractionAttenuation.Get_diffractionAttenuation_f(); 00468 00469 if (_displayConsoleStep == 1) 00470 { 00471 printf(" - low height correction Ch : %s \n", _withCh ? "ENABLED" : "DISABLED (Ch = 1)") ; 00472 printf(" - AttDiffH = %5.1f dB, avec pathDifference = %7.3f, DeltaDiff = %5.1f, ", _diffractionAttenuation_h, diffractionAttenuation.Get_pathDifferenceSR_h(), diffractionAttenuation.Get_DeltaDifSR_h()); 00473 printf("DeltaSol(S,O) = %5.1f, DeltaSol(O,R) = %5.1f, Asol(S,O) = %5.1f", diffractionAttenuation.Get_DeltaSolSO_h(), diffractionAttenuation.Get_DeltaSolOR_h(), diffractionAttenuation.Get_AttsolSO_h()); 00474 if(_path->GetOption(CHECK_EMBANKMENT)) 00475 { 00476 printf(", DeltaTalus(S,O) = %5.1f ", diffractionAttenuation.Get_DeltaTalusSO_h()); 00477 } 00478 printf("\n"); 00479 printf(" - AttDiffF = %5.1f dB, avec pathDifference = %7.3f, DeltaDiff = %5.1f, ", _diffractionAttenuation_f, diffractionAttenuation.Get_pathDifferenceSR_f(), diffractionAttenuation.Get_DeltaDifSR_f()); 00480 printf("DeltaSol(S,O) = %5.1f, DeltaSol(O,R) = %5.1f, Asol(S,O) = %5.1f", diffractionAttenuation.Get_DeltaSolSO_f(), diffractionAttenuation.Get_DeltaSolOR_f(), diffractionAttenuation.Get_AttsolSO_f()); 00481 if(_path->GetOption(CHECK_EMBANKMENT)) 00482 { 00483 printf(", DeltaTalus(S,O) = %5.1f ", diffractionAttenuation.Get_DeltaTalusSO_f()); 00484 } 00485 printf("\n"); 00486 00487 } 00488 00489 // if diffraction is negligible for homogeneous or favorable path, ground attenuation is calculated : 00490 if (_diffractionAttenuation_h == 0 || _diffractionAttenuation_f == 0) 00491 { 00492 GroundAttenuationCalculation(); 00493 EmbankmentAttenuationCalculation(); 00494 if (_diffractionAttenuation_h != 0) 00495 { 00496 _groundAttenuation_h = 0; 00497 _embankmentAttenuation_h = 0; 00498 } 00499 00500 if (_diffractionAttenuation_f != 0) 00501 { 00502 _groundAttenuation_f = 0; 00503 _embankmentAttenuation_f = 0; 00504 } 00505 } 00506 } 00507 00508 // Add absorption attenuation due to reflections 00509 _diffractionAttenuation_h += diffractionAttenuation.Get_absorptionAttenuation(); 00510 _diffractionAttenuation_f += diffractionAttenuation.Get_absorptionAttenuation(); 00511 } 00512 } 00513 else 00514 { 00515 // No diffractions or reflections 00516 _diffractionAttenuation_h = 0; 00517 _diffractionAttenuation_f = 0; 00518 00519 // ground attenuation calculation : 00520 GroundAttenuationCalculation(); 00521 00522 // embankment attenuation : 00523 EmbankmentAttenuationCalculation(); 00524 } 00525 00526 // border attenuation : 00527 _boundaryAttenuation_h = _diffractionAttenuation_h + _groundAttenuation_h + _embankmentAttenuation_h; 00528 _boundaryAttenuation_f = _diffractionAttenuation_f + _groundAttenuation_f + _embankmentAttenuation_f; 00529 00530 } 00531 00540 void Attenuation::AttenuationCalculation() 00541 { 00542 if (_displayConsoleStep == 1) 00543 { 00544 printf("Attenuation dans la bande d'octave %.1f : \n", _frequency); 00545 } 00546 00547 // Check SR < 2000 : 00548 if (_path->distSR > 2000) 00549 { 00550 printf ("Erreur: la distance de propagation (%.0fm) dépasse la limite de 2000m\n", _path->distSR) ; 00551 throw new int( ERRMaxDistance ); 00552 } 00553 00554 // Check side diffractions : 00555 int nbSideDiffractions = 0; 00556 for ( int i = 0 ; i < (int) _path->pathPoints.size() ; i++) 00557 { 00558 ProfilePointNMPB* profilePoint = _path->pathPoints[i]; 00559 if (profilePoint->ext != NULE && profilePoint->ext->type == ETSideDiffraction_NMPB) 00560 { 00561 nbSideDiffractions ++; 00562 } 00563 } 00564 00565 if(nbSideDiffractions > 2) 00566 { 00567 printf("Erreur : %d diffractions laterales detectees ; ce nombre ne doit pas depasser 2.\n", nbSideDiffractions); 00568 throw new int( ERRSideDiff ); 00569 } 00570 00571 double divergenceAttenuation = 0; 00572 if(!_path->GetOption(EXCLUDE_ADIV)) 00573 { 00574 if(nbSideDiffractions != 0) 00575 { 00576 // SR distance (or S'R' if reflections), go through reflections but no side diffractions 00577 double accuGroundDistance = 0; 00578 ProfilePointNMPB* precItem = _path->GetSource(); 00579 int nbPoints = (int) _path->pathPoints.size(); 00580 for (int i = 1 ; i < nbPoints ; i++) 00581 { 00582 ProfilePointNMPB* curItem = _path->pathPoints[i]; 00583 if (curItem->ext == NULE || curItem->ext->type != ETSideDiffraction_NMPB) 00584 { 00585 accuGroundDistance += GroundDistance(precItem->position3D, curItem->position3D); 00586 precItem = _path->pathPoints[i]; 00587 } 00588 } 00589 00590 Position2D receivPos; 00591 Position2D *recepPos = _path->GetReceiver()->position2D_withHeight(); 00592 receivPos.d = accuGroundDistance; 00593 receivPos.z = recepPos->z; 00594 delete recepPos; 00595 00596 Position2D *sourcePos = _path->GetSource()->position2D_withHeight(); 00597 double pathDirectLength = distance2D(sourcePos, &receivPos); 00598 delete sourcePos; 00599 00600 divergenceAttenuation = DivergenceAttenuationCalculation(pathDirectLength); 00601 } 00602 else 00603 { 00604 divergenceAttenuation = DivergenceAttenuationCalculation(_path->distSR); 00605 } 00606 } 00607 00608 double atmosphericAbsorption = 0; 00609 if(!_path->GetOption(EXCLUDE_AATM)) 00610 { 00611 atmosphericAbsorption = AtmosphericAbsorptionCalculation(_path->distSR, _frequency); 00612 } 00613 00614 BoundaryAttenuationCalculation(nbSideDiffractions); 00615 _attenuation_h = divergenceAttenuation + atmosphericAbsorption + _boundaryAttenuation_h; 00616 _attenuation_f = divergenceAttenuation + atmosphericAbsorption + _boundaryAttenuation_f; 00617 00618 if (_displayConsoleStep == 1) 00619 { 00620 printf(" - Aatm = %5.1f dB, Adiv = %5.1f dB \n", atmosphericAbsorption, divergenceAttenuation); 00621 printf(" - Attenuations totales : AttH = %5.1f dB, AttF = %5.1f dB \n", _attenuation_h, _attenuation_f); 00622 printf(" \n"); 00623 } 00624 } 00625 00629 void Attenuation::ClearResults() 00630 { 00631 _attenuation_h = 0; 00632 _attenuation_f = 0; 00633 _divergenceAttenuation = 0; 00634 _atmosphericAbsorption = 0; 00635 _boundaryAttenuation_h = 0; 00636 _boundaryAttenuation_f = 0; 00637 _groundAttenuation_h = 0; 00638 _groundAttenuation_f = 0; 00639 _diffractionAttenuation_h = 0; 00640 _diffractionAttenuation_f = 0; 00641 _embankmentAttenuation_h = 0; 00642 _embankmentAttenuation_f = 0; 00643 } 00644 00645 } 00646