NMPB08 Shared Libraries
CalculPropagation.cpp
Go to the documentation of this file.
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 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Defines