NMPB08 Shared Libraries
SousCalculs/Diffraction.cpp
Go to the documentation of this file.
00001 
00008 #include "Diffraction.h"
00009 #include <math.h>
00010 #include <stdio.h>
00011 #include <assert.h>
00012 
00013 namespace DiffractionNMPB
00014 {
00015 
00016         // ------------------------------ Path difference -------------------------------------------------
00017 
00027         double CurveRayLength(double distMN, double curvatureRadius)
00028         {
00029                 double curveLength = 0;
00030                 if (curvatureRadius != 0)
00031                 {
00032                         curveLength = 2 * curvatureRadius * asin(distMN / (2 * curvatureRadius));
00033                 }
00034 
00035                 return curveLength;
00036         }
00037 
00053         double PathDifference(Position2D const* source2D, Position2D const* receiver2D, vector<ProfilePointNMPB*> screenItems, bool favourableConditions)
00054         {
00055                 double pathDifference;
00056 
00057                 // gets all the kept screens :
00058                 int usedScreensNumber = 0;
00059                 vector<ProfilePointNMPB*> usedScreens;
00060                 for (int i = 0; i < (int) screenItems.size(); i++)
00061                 {
00062                         ProfilePointNMPB* screenToCheck = screenItems[i];
00063                         if (screenToCheck != 0 && screenToCheck->isDiff == 1)
00064                         {
00065                                 usedScreens.push_back(screenToCheck);
00066                                 usedScreensNumber++;
00067                         }
00068                 }
00069 
00070                 if (usedScreensNumber == 0)
00071                 {
00072                         return 0;
00073                 }
00074 
00075                 double distSR = distance2D(source2D, receiver2D);
00076                 if (distSR == 0)
00077                 {
00078                         return 0;
00079                 }
00080 
00081                 if(favourableConditions)
00082                 {
00083                         // favorable conditions (p.47-49 - § 9.4.3.2)
00084                         double curvatureRadius = 8 * distSR; // p.47 - § 9.4.3.2 - Formula (33)
00085                         if (curvatureRadius < 1000)
00086                         {
00087                                 curvatureRadius = 1000;
00088                         }
00089 
00090                         if (usedScreensNumber == 1)
00091                         {
00092                                 // only one screen 
00093                                 ProfilePointNMPB* screen = usedScreens[0];
00094                                 Position2D *posScreen = screen->position2D_withHeight();
00095                                 if (screen->h_ray <= 0)
00096                                 {
00097                                         
00098                                         pathDifference = CurveRayLength(distance2D(source2D, posScreen), curvatureRadius) + CurveRayLength(distance2D(posScreen, receiver2D), curvatureRadius) - CurveRayLength(distSR, curvatureRadius);
00099                                 }
00100                                 else
00101                                 {
00102                                         Position2D pointA;
00103                                         pointA.d = posScreen->d;
00104                                         pointA.z = posScreen->z + screen->h_ray;
00105                                         pathDifference = 2 * CurveRayLength(distance2D(source2D, &pointA), curvatureRadius) + 2 * CurveRayLength(distance2D(&pointA, receiver2D), curvatureRadius) - CurveRayLength(distance2D(source2D, posScreen), curvatureRadius) - CurveRayLength(distance2D(posScreen, receiver2D), curvatureRadius) - CurveRayLength(distSR, curvatureRadius);
00106                                 }
00107                                 delete posScreen;
00108                         }
00109                         else
00110                         {
00111                                 // multiple screens 
00112                                 double distScreens = 0;
00113                                 Position2D *posScreen1, *posScreen2;
00114                                 for (int i = 1; i < usedScreensNumber; i++)
00115                                 {
00116                                         posScreen1 = usedScreens[i-1]->position2D_withHeight();
00117                                         posScreen2 = usedScreens[i]->position2D_withHeight();
00118                                         distScreens += CurveRayLength(distance2D(posScreen1,posScreen2), curvatureRadius);
00119                                         
00120                                         delete posScreen1;
00121                                         delete posScreen2;
00122                                 }
00123                                 posScreen1 = usedScreens[0]->position2D_withHeight();
00124                                 posScreen2 = usedScreens[usedScreensNumber-1]->position2D_withHeight();
00125                                 pathDifference = CurveRayLength(distance2D(source2D, posScreen1), curvatureRadius) + distScreens + CurveRayLength(distance2D(posScreen2, receiver2D), curvatureRadius) - CurveRayLength(distSR, curvatureRadius);
00126                                 delete posScreen1;
00127                                 delete posScreen2;
00128                         }               
00129 
00130                 }
00131                 else
00132                 {
00133                         // homogeneous conditions (p.47 - § 9.4.3.1)
00134                         if (usedScreensNumber == 1)
00135                         {
00136                                 // only one screen 
00137                                 ProfilePointNMPB* screen = usedScreens[0];
00138                                 Position2D *posScreen = screen->position2D_withHeight();
00139                                 pathDifference = distance2D(source2D, posScreen) + distance2D(posScreen, receiver2D) - distSR;
00140                                 
00141                                 /*
00142                                  * BUG CORRECTED
00143                                  *
00144                                  * DVM 28/11/2011
00145                                  *
00146                                  * screen->h_ray is valid only for real source to real receiver
00147                                  *
00148                                  * as we also use this calculation for image sources and/or receivers (e.g. when 
00149                                  * estimating ground effects left and right of a diffracting edge), we must recalculate
00150                                  * the height of the ray-path for every situation...
00151                                  */
00152 
00153                                 double d1 = source2D->d ;
00154                                 double z1 = source2D->z ;
00155                                 double d2 = receiver2D->d ;
00156                                 double z2 = receiver2D->z ;
00157                                 double dd = posScreen->d ;
00158                                 double zd = posScreen->z ;
00159                                 double zr = z1 + (z2 - z1) * (dd - d1) / (d2 - d1) ;
00160 
00161                                 delete posScreen;
00162 
00163                                 //if (screen->h_ray > 0)
00164                                 if (zr > zd)
00165                                 {
00166                                         pathDifference = -1 * pathDifference;                                   
00167                                 }
00168                         }
00169                         else
00170                         {
00171                                 // multiple screens 
00172                                 double distScreens = 0;
00173                                 Position2D *posScreen1, *posScreen2;
00174                                 for (int i = 1; i < usedScreensNumber; i++)
00175                                 {
00176                                         posScreen1 = usedScreens[i-1]->position2D_withHeight();
00177                                         posScreen2 = usedScreens[i]->position2D_withHeight();
00178                                         distScreens += distance2D(posScreen1,posScreen2);
00179                                         delete posScreen1;
00180                                         delete posScreen2;
00181                                 }
00182                                 posScreen1 = usedScreens[0]->position2D_withHeight();
00183                                 posScreen2 = usedScreens[usedScreensNumber-1]->position2D_withHeight();
00184                                 pathDifference = distance2D(source2D, posScreen1) + distScreens + distance2D(posScreen2, receiver2D) - distSR;
00185                                 delete posScreen1;
00186                                 delete posScreen2;
00187                         }               
00188                 }
00189 
00190                 return pathDifference;
00191         }
00192 
00206         double PathDifference(Position2D const* source2D, Position2D const* receiver2D, ProfilePointNMPB* reflectionItem)
00207         {
00208                 double pathDifference;
00209 
00210                 double distSR = distance2D(source2D, receiver2D);
00211                 if (distSR == 0)
00212                 {
00213                         return 0;
00214                 }
00215 
00216                 reflectionItem->height = reflectionItem->ext->height;
00217                 Position2D *reflexPos = reflectionItem->position2D_withHeight();
00218                 pathDifference = -1 * (distance2D(source2D, reflexPos) + distance2D(reflexPos, receiver2D) - distSR);
00219                 delete reflexPos;
00220 
00221                 return pathDifference;
00222         }
00223 
00239         double SidePathDifference(Position3D const* source3D, Position3D const* receiver3D, vector<ProfilePointNMPB*> screenItems, double &totalDiffDist)
00240         {
00241                 double pathDifference = 0;
00242 
00243                 // gets all the kept screens : 
00244                 // 1 or 2 edges and the last reflection screen before edges and the first reflection screen after edges
00245                 int usedEdgesNumber = 0;
00246                 int usedRefNumber = 0;
00247                 vector<ProfilePointNMPB*> usedEdges;
00248                 vector<ProfilePointNMPB*> usedReflections;
00249                 int posFirstEdge = 0;
00250                 int posLastEdge = 0;
00251                 int posLastRefBefore = 0;
00252                 int posFirstRefAfter = 0;
00253                 for (int i = 0; i < (int) screenItems.size(); i++)
00254                 {
00255                         ProfilePointNMPB* screenToCheck = screenItems[i];
00256                         if (screenToCheck != 0 && screenToCheck->ext != 0)
00257                         {
00258                                 if(screenToCheck->ext->type == ETSideDiffraction_NMPB)
00259                                 {
00260                                         if(posFirstEdge == 0)
00261                                         {
00262                                                 posFirstEdge = i;
00263                                         }
00264 
00265                                         posLastEdge = i;
00266 
00267                                         usedEdges.push_back(screenToCheck);
00268                                         usedEdgesNumber++;
00269                                 }
00270                                 else if(screenToCheck->ext->type == ETReflection_NMPB)
00271                                 {
00272                                         if(posFirstEdge == 0)
00273                                         {
00274                                                 posLastRefBefore = i;
00275                                                 usedReflections.push_back(screenToCheck);
00276                                                 usedRefNumber++;
00277                                         }
00278                                         else if(posLastEdge != 0 && posFirstRefAfter == 0)
00279                                         {
00280                                                 posFirstRefAfter = i;
00281                                                 usedReflections.push_back(screenToCheck);
00282                                                 usedRefNumber++;
00283                                         }
00284                                 }
00285                         }
00286                 }
00287 
00288                 if (usedEdgesNumber == 0)
00289                 {
00290                         return 0;
00291                 }
00292                 
00293                 // SR or S'R' calculation (in reflection case) 
00294                 // All the calculations are made in horizontal plane, so z = 0
00295                 double distSR = 0;
00296                 Position3D sourceImage;
00297                 Position3D receiverImage;
00298                 sourceImage.x = source3D->x;
00299                 sourceImage.y = source3D->y;
00300                 sourceImage.z = 0;
00301                 receiverImage.x = receiver3D->x;
00302                 receiverImage.y = receiver3D->y;
00303                 receiverImage.z = 0;
00304                 if(usedRefNumber > 0)
00305                 {
00306                         if (posLastRefBefore != 0 && posLastRefBefore < posFirstEdge)
00307                         {
00308                                 ProfilePointNMPB* firstEdge = usedEdges[0];
00309                                 ProfilePointNMPB* refBefore = usedReflections[0];
00310                                 double mult = (firstEdge->position2D->d - screenItems[0]->position2D->d) / (firstEdge->position2D->d - refBefore->position2D->d);
00311                                 sourceImage.x = firstEdge->position3D->x  + (refBefore->position3D->x - firstEdge->position3D->x) * mult;
00312                                 sourceImage.y = firstEdge->position3D->y  + (refBefore->position3D->y - firstEdge->position3D->y) * mult;               
00313                         }
00314 
00315                         if (posFirstRefAfter != 0 && posFirstRefAfter > posLastEdge)
00316                         {
00317                                 ProfilePointNMPB* lastEdge = usedEdges[usedEdgesNumber - 1];
00318                                 ProfilePointNMPB* refAfter = usedReflections[usedRefNumber - 1];
00319                                 double mult = (lastEdge->position2D->d - screenItems[screenItems.size() - 1]->position2D->d) / (lastEdge->position2D->d - refAfter->position2D->d);
00320                                 receiverImage.x = lastEdge->position3D->x  + (refAfter->position3D->x - lastEdge->position3D->x) * mult;
00321                                 receiverImage.y = lastEdge->position3D->y  + (refAfter->position3D->y - lastEdge->position3D->y) * mult;                                
00322                         }
00323                 }
00324 
00325                 distSR = GroundDistance(&sourceImage, &receiverImage);
00326 
00327                 if (distSR == 0)
00328                 {
00329                         return 0;
00330                 }
00331 
00332                 // Path difference in horizontal plane (only homogeneous conditions)
00333                 if (usedEdgesNumber == 1)
00334                 {
00335                         // only one edge 
00336                         pathDifference = GroundDistance(&sourceImage, usedEdges[0]->position3D) + GroundDistance(usedEdges[0]->position3D, &receiverImage) - distSR;
00337                         totalDiffDist = 0;
00338                 }
00339                 else
00340                 {
00341                         // multiple screens 
00342                         totalDiffDist = usedEdges[usedEdgesNumber-1]->position2D->d - usedEdges[0]->position2D->d;
00343                         pathDifference = GroundDistance(&sourceImage, usedEdges[0]->position3D) + totalDiffDist + GroundDistance(usedEdges[usedEdgesNumber-1]->position3D, &receiverImage) - distSR;
00344                 }                
00345 
00346                 return pathDifference;
00347         }
00348 
00349 
00350         // ------------------------------ Diffraction Class -----------------------------------------
00351 
00369         Diffraction::Diffraction(vector<ProfilePointNMPB*> terrainItems, double freq, int freqPos, bool withEmbankment)
00370         {
00371                 if (terrainItems.size() < 2)
00372                 {
00373                         return;
00374                 }
00375 
00376                 _source = terrainItems[0];
00377                 _receiver = terrainItems[terrainItems.size() - 1];
00378                 _terrainItems = terrainItems;
00379                 _freq = freq;
00380 
00381                 if (freq != 0)
00382                 {
00383                         _wavelength =  cSound / freq;
00384                 }
00385                 else
00386                 {
00387                         printf (".ERROR: frequency = 0 \n") ;
00388                         throw new int( ERRFrequency );
00389                 }
00390 
00391                 ClearData();
00392 
00393                 // reflections attenuation
00394                 ReflectionAttenuation(terrainItems, freqPos);
00395 
00396                 // path differences
00397                 Position2D *posSource = _source->position2D_withHeight();
00398                 Position2D *posRecep = _receiver->position2D_withHeight();
00399                 _pathDifferenceSR_h = PathDifference(posSource, posRecep, _terrainItems, false);
00400                 _pathDifferenceSR_f = PathDifference(posSource, posRecep, _terrainItems, true);
00401                 delete posSource;
00402                 delete posRecep;
00403 
00404                 MeanPlanesDataCalculation(withEmbankment);
00405         }
00406 
00425         Diffraction::Diffraction(vector<ProfilePointNMPB*> terrainItems, double freq, int freqPos, int nbSideDiffractions)
00426         {
00427                 if (terrainItems.size() < 2)
00428                 {
00429                         return;
00430                 }
00431 
00432                 _source = terrainItems[0];
00433                 _receiver = terrainItems[terrainItems.size() - 1];
00434                 _terrainItems = terrainItems;
00435                 _freq = freq;
00436 
00437                 if (freq != 0)
00438                 {
00439                         _wavelength =  cSound / freq;
00440                 }
00441                 else
00442                 {
00443                         printf (".ERROR: frequency = 0 \n") ;
00444                         throw new int( ERRFrequency );
00445                 }
00446 
00447                 ClearData();
00448 
00449                 // reflections attenuation
00450                 ReflectionAttenuation(terrainItems, freqPos);
00451 
00452                 // side diffraction
00453                 Position3D *terrainPos1, *terrainPos2;
00454                 Position2D *sourcePos, *recepPos;
00455                 terrainPos1 = terrainItems[0]->position3D_withHeight();
00456                 terrainPos2 = terrainItems[(int)terrainItems.size()-1]->position3D_withHeight();
00457                 sourcePos = _source->position2D_withHeight();
00458                 recepPos = _receiver->position2D_withHeight();
00459                 _pathDifferenceSR_h = SidePathDifference(terrainPos1, terrainPos2, terrainItems, _totalDiffDist);
00460                 _DeltaDifSR_h = CalculDiffractionPure(sourcePos , recepPos, _terrainItems, false, true, false);
00461                 delete terrainPos1;
00462                 delete terrainPos2;
00463                 delete sourcePos;
00464                 delete recepPos;
00465         }
00466 
00474         void Diffraction::ReflectionAttenuation(vector<ProfilePointNMPB*> &terrainItems, int freqPos) 
00475         {
00476                 _absorptionAttenuation = 0;
00477                 for (int i = 0 ; i < (int)terrainItems.size() ; i++)
00478                 {
00479                         ProfilePointNMPB* pointNMPB = terrainItems[i];
00480                         if(pointNMPB->ext->type == ETReflection_NMPB)
00481                         {
00482                                 _absorptionAttenuation += SoundAbsorption(pointNMPB, i, freqPos);
00483                         }
00484                 }
00485                 // * (-1) because value taken as attenuation 
00486                 _absorptionAttenuation = -1. * _absorptionAttenuation;
00487         }
00488 
00502         double Diffraction::SoundAbsorption(ProfilePointNMPB* pointNMPB, int posRef, int freqPos)
00503         {
00504                 // absorption p.51 - § 9.5.1 Formula (43)
00505                 double* alphaArray = pointNMPB->ext->alphaArray;
00506                 if(freqPos >= 0 && alphaArray != NULL)
00507                 {
00508                         double absorpCoeff = alphaArray[freqPos];
00509 
00510                         double screenAbsorption = 0;
00511                         if (absorpCoeff < 1)
00512                         {
00513                                 screenAbsorption =  10 * log10(1-absorpCoeff);
00514                         }
00515                         else
00516                         {
00517                                 printf (".ERROR: screen absorption %4.2f > 1 \n", screenAbsorption) ;
00518                                 throw new int( ERRScreenAbsorption );
00519                         }
00520 
00521                         // retro diffraction p.52-53-54 - § 9.5.2
00522                         double retroDiff = 0;
00523                         if (pointNMPB->ext->height != 0)
00524                         {
00525                                 pointNMPB->height = pointNMPB->ext->height;
00526 
00527                                 // retro diffraction between S and R, 
00528                                 //      or between E1 and E2 if E1  last diffraction screen between S and the reflection screen
00529                                 //                                                      E2 first diffraction screen between the reflection screen and R
00530                                 ProfilePointNMPB* diffPointBefore = _source;
00531                                 for (int i = 1; i < posRef; i++)
00532                                 {
00533                                         ProfilePointNMPB* screenToCheck = _terrainItems[i];
00534                                         if (screenToCheck != 0 && screenToCheck->isDiff == 1 && screenToCheck->h_ray <= 0)
00535                                         {
00536                                                 diffPointBefore = screenToCheck;
00537                                         }
00538                                 }
00539 
00540                                 ProfilePointNMPB* diffPointAfter = _receiver;
00541                                 for (int i = posRef + 1; i < (int) _terrainItems.size(); i++)
00542                                 {
00543                                         ProfilePointNMPB* screenToCheck = _terrainItems[i];
00544                                         if (screenToCheck != 0 && screenToCheck->isDiff == 1 && screenToCheck->h_ray <= 0)
00545                                         {
00546                                                 diffPointAfter = screenToCheck;
00547                                                 break;
00548                                         }
00549                                 }
00550 
00551                                 Position2D *posBefore, *posAfter;
00552                                 posBefore = diffPointBefore->position2D_withHeight();
00553                                 posAfter = diffPointAfter->position2D_withHeight();
00554                                 retroDiff = CalculDiffractionPure(posBefore, posAfter, pointNMPB);
00555                                 delete posBefore;
00556                                 delete posAfter;
00557                         }
00558 
00559                         //p.53 - § 9.5.2 Formula (46)
00560                         return screenAbsorption - retroDiff;
00561                 }
00562 
00563                 return 0;
00564         }
00565 
00569         void Diffraction::MeanPlanesDataCalculation(bool withEmbankment)
00570         {
00571                 // search first and last used screens :
00572                 _firstScreen = 0;
00573                 _lastScreen = 0;
00574                 int screenNumber = 0;
00575                 int firstScreenPos = 0;
00576                 int lastScreenPos = 0;
00577                 for (int i = 0; i < (int) _terrainItems.size(); i++)
00578                 {
00579                         ProfilePointNMPB* screenToCheck = _terrainItems[i];
00580                         if (screenToCheck != 0 && screenToCheck->isDiff == 1)
00581                         {
00582                                 screenNumber++;
00583                                 if (_firstScreen == 0)
00584                                 {
00585                                         _firstScreen = screenToCheck;
00586                                         firstScreenPos = i;
00587                                 }
00588 
00589                                 _lastScreen = screenToCheck;
00590                                 lastScreenPos = i;
00591                         }
00592                 }
00593                 
00594                 if (screenNumber == 0)
00595                 {
00596                         _Adif_h = 0;
00597                         _Adif_f = 0;
00598                         return;
00599                 }
00600 
00601                 _totalDiffDist = 0;
00602                 if (screenNumber > 1)
00603                 {
00604                         Position2D *firstPos, *lastPos;
00605                         firstPos = _firstScreen->position2D_withHeight();
00606                         lastPos = _lastScreen->position2D_withHeight();
00607                         _totalDiffDist = distance2D(firstPos , lastPos);
00608                         delete firstPos;
00609                         delete lastPos;
00610                 }
00611 
00612                 // Mean planes (p.45 figure 21)
00613                 vector<ProfilePointNMPB*> terrainsS;
00614                 vector<ProfilePointNMPB*> terrainsR;
00615                 for (int i = 0 ; i < (int) _terrainItems.size(); i++ )
00616                 {
00617                         ProfilePointNMPB* terrain = _terrainItems[i];
00618                         if (i <= firstScreenPos)
00619                         {
00620                                 terrainsS.push_back(terrain);
00621                         }
00622 
00623                         if (i >= lastScreenPos)
00624                         {
00625                                 terrainsR.push_back(terrain);
00626                         }
00627                 }
00628 
00629                 MeanPlane meanPlaneSO(terrainsS);
00630                 _meanPlaneSO = meanPlaneSO;
00631                 MeanPlane meanPlaneOR(terrainsR);
00632                 _meanPlaneOR = meanPlaneOR;
00633 
00634                 _h0 = max(meanPlaneSO.Get_zEqReceiver(), meanPlaneOR.Get_zEqSource());  
00635 
00636                 // Embankment calculations
00637                 _AttTalusSO = 0;
00638                 if(withEmbankment)
00639                 {
00640                         Embankment embankment(terrainsS, _freq);
00641                         _AttTalusSO = embankment.GetEmbankmentAttenuation();
00642                         _embankmentSourceImage.d = embankment.GetSourceImage()->d;
00643                         _embankmentSourceImage.z = embankment.GetSourceImage()->z;
00644                 }
00645         }
00646 
00666         double Diffraction::CalculDiffractionPure (Position2D const* source2D, 
00667                                                                                            Position2D const* receiver2D, 
00668                                                                                            vector<ProfilePointNMPB*> screenItems, 
00669                                                                                            bool favourableConditions, 
00670                                                                                            bool totalPath, 
00671                                                                                            bool withCh)
00672         {
00673                 // calculation of coefficient for multiple diffractions (C'')
00674                 // Formula (32)
00675                 double multiDiffCoef = 1;
00676                 if (_totalDiffDist >= 0.3)
00677                 {
00678                         double val = pow(5. * _wavelength / _totalDiffDist, 2);
00679                         multiDiffCoef = (1. + val) / (1./3. + val);
00680                 }
00681 
00682                 // calculation of Ch coefficient
00683                 double Ch = 1;
00684                 if(withCh)
00685                 {
00686                         double Chtemp = _freq * _h0 / 250.0;
00687                         if (Chtemp < 1)
00688                         {
00689                                 Ch = Chtemp;
00690                         }
00691                 }
00692 
00693                 // path difference calculation
00694                 double pathDifference;
00695                 if (totalPath)
00696                 {
00697                         // SR path already calculated
00698                         if (favourableConditions)
00699                         {
00700                                 pathDifference = _pathDifferenceSR_f;
00701                         }
00702                         else
00703                         {
00704                                 pathDifference = _pathDifferenceSR_h;
00705                         }
00706                 }
00707                 else
00708                 {
00709                         pathDifference = PathDifference(source2D, receiver2D, screenItems, favourableConditions);
00710                 }
00711 
00712 
00713                 // diffraction calculation 
00714                 double deltaDif = 0;
00715                 double comparisonVal = 40 * multiDiffCoef * pathDifference / _wavelength;
00716                 if (comparisonVal >= -2)
00717                 {
00718                         deltaDif = 10 * Ch * log10(3. + comparisonVal);
00719                         if (deltaDif < 0) deltaDif = 0;
00720                 }
00721         
00722                 return deltaDif;
00723         }
00724 
00738         double Diffraction::CalculDiffractionPure (Position2D const* source2D, 
00739                                                                                            Position2D const* receiver2D, 
00740                                                                                            ProfilePointNMPB* reflectionItem)
00741         {
00742                 double pathDifference =  PathDifference(source2D, receiver2D, reflectionItem);
00743 
00744                 // diffraction calculation 
00745                 double deltaDif = 0;
00746                 if (_wavelength != 0)
00747                 {
00748                         double comparisonVal = 40 * pathDifference / _wavelength;
00749                         if (comparisonVal >= -2)
00750                         {
00751                                 deltaDif = 10 * log10(3. + comparisonVal);
00752                         }
00753                 }
00754 
00755                 return deltaDif;
00756         }
00757 
00769         double Diffraction::CalculDeltaSol(double aSol, double deltaDif1, double deltaDif2)
00770         {
00771                 /*
00772                  * BUG CORRECTION
00773                  *
00774                  * DVM 28/11/2011
00775                  *
00776                  * weighting of ground effect by means of diffraction may become undefined if the 
00777                  * image source (or receiver) is less attenuated than the real source (or receiver) ;
00778                  * this could happen if the source or receiver was found below the mean plane for
00779                  * the ground segments on either side of a diffracting edge
00780                  *
00781                  * assert added in case this happens... no solution proposed
00782                  */
00783                 assert (deltaDif1 >= deltaDif2) ;
00784 
00785                 double deltaSol = -20 * log10(1 + (pow(10, -1 * aSol / 20.0) - 1) * pow(10, -1 * (deltaDif1 - deltaDif2) / 20.0));
00786                 return deltaSol;
00787         }
00788 
00796         void Diffraction::CalculAttenuationDiffraction (bool withCh)
00797         {
00798                 // Checks path difference :
00799                 double minPathDifference = -1 * _wavelength / 20;
00800                 bool doHomogeneous = _pathDifferenceSR_h >= minPathDifference;
00801                 bool doFavorable = _pathDifferenceSR_f >= minPathDifference;
00802 
00803                 if (!doHomogeneous && !doFavorable)
00804                 {
00805                         _Adif_h = 0;
00806                         _Adif_f = 0;
00807                         return;
00808                 }
00809 
00810                 // intermediate calculations for ground
00811                 GroundEffect grEffSO(_meanPlaneSO.Get_dp(),_meanPlaneSO.Get_zEqSource(), _meanPlaneSO.Get_zEqReceiver(), _meanPlaneSO.Get_Gpath(), _meanPlaneSO.Get_Gsource(), _freq, DeltaSol_SO);
00812                 GroundEffect grEffOR(_meanPlaneOR.Get_dp(),_meanPlaneOR.Get_zEqSource(), _meanPlaneOR.Get_zEqReceiver(), _meanPlaneOR.Get_Gpath(), _meanPlaneOR.Get_Gsource(), _freq, DeltaSol_OR);
00813 
00814                 Position2D *posSource, *posRecep;
00815                 posSource = _source->position2D_withHeight();
00816                 posRecep = _receiver->position2D_withHeight();
00817                 // homogeneous case
00818                 if(doHomogeneous)
00819                 {
00820                         // calculation of pure diffraction between source and receiver :
00821                         _DeltaDifSR_h = CalculDiffractionPure(posSource, posRecep, _terrainItems, false, true, withCh);
00822 
00823                         // calculation of ground attenuation with diffraction, source side: formula (38)
00824                         _aSolSO_h = grEffSO.AttenuationCalculationH();
00825                         Position2D sourceImage = _meanPlaneSO.Get_sourceImage();
00826                         double deltaDifSO_h = CalculDiffractionPure(&sourceImage, posRecep, _terrainItems, false, false, withCh);
00827                         _DeltaSolSO_h = CalculDeltaSol(_aSolSO_h, deltaDifSO_h, _DeltaDifSR_h);
00828                         
00829                         // calculation of ground attenuation with diffraction, receiver side: formula (40)
00830                         double aSolOR_h = grEffOR.AttenuationCalculationH();
00831                         Position2D recepImage = _meanPlaneOR.Get_receiverImage();
00832                         double deltaDifOR_h = CalculDiffractionPure(posSource, &recepImage, _terrainItems, false, false, withCh);
00833                         _DeltaSolOR_h = CalculDeltaSol(aSolOR_h, deltaDifOR_h, _DeltaDifSR_h);
00834                         
00835                         // calculation of embankment source side attenuation : formula (39)
00836                         _DeltaTalusSO_h = 0;
00837                         if (_AttTalusSO != 0)
00838                         {
00839                                 double deltaDifTalusSO_h = CalculDiffractionPure (&_embankmentSourceImage, posRecep, _terrainItems, false, false, withCh);
00840                                 _DeltaTalusSO_h = CalculDeltaSol(_AttTalusSO, deltaDifTalusSO_h, _DeltaDifSR_h);
00841                         }
00842 
00843                         // calculation of the total diffraction attenuation : formula (37)
00844                         if (_DeltaDifSR_h > 25)
00845                         {
00846                                 // p.46 - § 9.4.2 
00847                                 _DeltaDifSR_h = 25;
00848                         }
00849 
00850                         _Adif_h = _DeltaDifSR_h + _DeltaSolSO_h + _DeltaTalusSO_h + _DeltaSolOR_h;
00851                 }
00852                 else
00853                 {
00854                         _Adif_h = 0;
00855                 }
00856 
00857                 // favorable case
00858                 if(doFavorable)
00859                 {
00860                         // calculation of pure diffraction between source and receiver :
00861                         _DeltaDifSR_f = CalculDiffractionPure(posSource, posRecep, _terrainItems, true, true, withCh);
00862 
00863                         // calculation of ground attenuation with diffraction, source side : formula (38)
00864                         _aSolSO_f = grEffSO.AttenuationCalculationF();
00865                         Position2D sourceImage = _meanPlaneSO.Get_sourceImage();
00866                         double deltaDifSO_f = CalculDiffractionPure(&sourceImage, posRecep, _terrainItems, true, false, withCh);
00867                         _DeltaSolSO_f = CalculDeltaSol(_aSolSO_f, deltaDifSO_f, _DeltaDifSR_f);
00868 
00869                         // calculation of ground attenuation with diffraction, receiver side: formula (40)
00870                         double aSolOR_f = grEffOR.AttenuationCalculationF();
00871                         Position2D recepImage = _meanPlaneOR.Get_receiverImage();
00872                         double deltaDifOR_f = CalculDiffractionPure(posSource, &recepImage, _terrainItems, true, false, withCh);
00873                         _DeltaSolOR_f = CalculDeltaSol(aSolOR_f, deltaDifOR_f, _DeltaDifSR_f);
00874 
00875                         // calculation of embankment source side attenuation : formula (39)
00876                         _DeltaTalusSO_f = 0;
00877                         if (_AttTalusSO != 0)
00878                         {
00879                                 double deltaDifTalusSO_f = CalculDiffractionPure(&_embankmentSourceImage, posRecep, _terrainItems, true, false, withCh);
00880                                 _DeltaTalusSO_f = CalculDeltaSol(_AttTalusSO, deltaDifTalusSO_f, _DeltaDifSR_f);
00881                         }
00882 
00883                         // calculation of the total diffraction attenuation : formula (37)
00884                         if (_DeltaDifSR_f > 25)
00885                         {
00886                                 // p.46 - § 9.4.2 
00887                                 _DeltaDifSR_f = 25;
00888                         }
00889 
00890                         _Adif_f = _DeltaDifSR_f + _DeltaSolSO_f + _DeltaTalusSO_f + _DeltaSolOR_f;
00891                 }
00892                 else
00893                 {
00894                         _Adif_f = 0;
00895                 }
00896                 
00897                 delete posRecep;
00898                 delete posSource;
00899         }
00900 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Defines