NMPB08 Shared Libraries
SousCalculs/GroundEffect.cpp
Go to the documentation of this file.
00001 
00008 #include "GroundEffect.h"
00009 #include <math.h>
00010 #include <stdio.h>
00011 
00012 namespace GroundEffectNMPB
00013 {       
00014         // ------------------------------ MeanPlane Class -----------------------------------------
00015 
00021         MeanPlane::MeanPlane(vector<ProfilePointNMPB*> terrainItems)
00022         {
00023                 if (terrainItems.size() < 2)
00024                 {
00025                         return;
00026                 }
00027 
00028                 _source = terrainItems[0];
00029                 _receiver = terrainItems[terrainItems.size() - 1];
00030                 _terrainItems = terrainItems;
00031 
00032                 // Fills the points list with the terrain items positions
00033                 for (int i = 0 ; i < (int) terrainItems.size() ; i++)
00034                 {
00035                         ProfilePointNMPB* terrainItem = terrainItems[i];
00036                         _pointsList.push_back(*(terrainItem->position2D));
00037                 }
00038 
00039                 CalculateData();
00040 
00041         }
00042 
00049         void MeanPlane::FillLineCoefficients()
00050         {
00051                 int pointsNumber = (int) _pointsList.size();
00052                 if (_pointsList[0].d != _pointsList[pointsNumber-1].d)
00053                 {
00054                         // mean plane coefficients
00055                         double valA1 = 0;
00056                         double valA2 = 0;
00057                         double valB1 = 0;
00058                         double valB2 = 0;
00059                         for (int i = 0 ; i < pointsNumber - 1 ; i++)
00060                         {
00061                                 Position2D point_i = _pointsList[i];
00062                                 Position2D point_i1 = _pointsList[i+1];
00063                                 double diffdi = point_i1.d - point_i.d;
00064                                 if (diffdi != 0)
00065                                 {
00066                                         // zi = ai * di + bi
00067                                         double ai = (point_i1.z - point_i.z) / diffdi;
00068                                         double bi = point_i.z - ai * point_i.d;
00069                                         double vald2 = (pow(point_i1.d, 2) - pow(point_i.d, 2));
00070                                         valA1 += ai * (pow(point_i1.d, 3) - pow(point_i.d, 3));
00071                                         valA2 += bi * vald2;
00072                                         valB1 += ai * vald2;
00073                                         valB2 += bi * diffdi;
00074                                 }                       
00075                         }
00076 
00077                         double valA = 2 * valA1 / 3 + valA2;
00078                         double valB = valB1 + 2 * valB2;
00079 
00080                         double dividende = pow(_pointsList[pointsNumber-1].d - _pointsList[0].d, 3);
00081                         _aCoeff =  3 * (2 * valA - valB * (_pointsList[pointsNumber-1].d + _pointsList[0].d)) / dividende;
00082                         _bCoeff = 2 * valB * (pow(_pointsList[pointsNumber-1].d, 3) - pow(_pointsList[0].d, 3))/ ((_pointsList[pointsNumber-1].d - _pointsList[0].d) * dividende) - 3 * (_pointsList[pointsNumber-1].d + _pointsList[0].d) * valA / dividende; 
00083                 }       
00084         }
00085 
00093         Position2D MeanPlane::MeanPlaneProjection(Position2D const* point)
00094         {
00095                 Position2D eqPoint;
00096                 double dividende = 1 + _aCoeff * _aCoeff;
00097                 eqPoint.d = (point->z * _aCoeff + point->d - _aCoeff * _bCoeff) / dividende;
00098                 eqPoint.z = (point->z * _aCoeff * _aCoeff + point->d * _aCoeff + _bCoeff) / dividende;
00099 
00100                 return eqPoint;
00101         }
00102 
00109         void MeanPlane::FillData()
00110         {
00111                 Position2D _sourceWithHeight;
00112                 _sourceWithHeight.d = _source->position2D->d;
00113                 _sourceWithHeight.z = _source->position2D->z + _source->height;
00114                 Position2D _receiverWithHeight;
00115                 _receiverWithHeight.d = _receiver->position2D->d;
00116                 _receiverWithHeight.z = _receiver->position2D->z + _receiver->height;
00117 
00118                 Position2D eq1 = MeanPlaneProjection(&_sourceWithHeight);
00119                 Position2D eq2 = MeanPlaneProjection(&_receiverWithHeight);
00120 
00121                 if (eq1.z <= _sourceWithHeight.z)
00122                 {
00123                         _zEqS = distance2D(&eq1, &_sourceWithHeight);
00124                         _imageS.d = 2 * eq1.d - _sourceWithHeight.d;
00125                         _imageS.z = 2 * eq1.z - _sourceWithHeight.z;
00126                 }
00127                 else
00128                 {
00129                         _zEqS = 0;
00130                         _imageS = eq1;
00131                 }
00132 
00133                 if (eq2.z <= _receiverWithHeight.z)
00134                 {
00135                         _zEqR = distance2D(&eq2, &_receiverWithHeight);
00136                         _imageR.d = 2 * eq2.d - _receiverWithHeight.d;
00137                         _imageR.z = 2 * eq2.z - _receiverWithHeight.z;
00138                 }
00139                 else
00140                 {
00141                         _zEqR = 0;
00142                         _imageR = eq2;
00143                 }
00144 
00145                 _dp = distance2D(&eq1, &eq2);
00146 
00147                 // ground coefficient calculation : 
00148                 _Gpath = 0;
00149                 int terrainsNb = (int) _terrainItems.size();
00150                 if (terrainsNb != 0 )
00151                 {               
00152                         double sumG = 0;
00153                         double sumD = 0;
00154                         for (int i = 0; i < terrainsNb ; i++ )
00155                         {
00156                                 ProfilePointNMPB* terrainItem = _terrainItems[i];
00157                                 if(i == 0)
00158                                 {
00159                                         sumD += terrainItem->position2D->d - _source->position2D->d;
00160                                         sumG += terrainItem->impedance * (terrainItem->position2D->d - _source->position2D->d);
00161                                 }
00162                                 else
00163                                 {
00164                                         sumD += terrainItem->position2D->d - _terrainItems[i-1]->position2D->d;
00165                                         sumG += terrainItem->impedance * (terrainItem->position2D->d - _terrainItems[i-1]->position2D->d);
00166                                 }
00167                         }
00168 
00169                         sumD += _receiver->position2D->d - _terrainItems[terrainsNb-1]->position2D->d;
00170                         sumG += _receiver->impedance * (_receiver->position2D->d - _terrainItems[terrainsNb-1]->position2D->d);
00171                         
00172                         if(sumD != 0)
00173                         {
00174                                 _Gpath = sumG / sumD;
00175                         }               
00176                 }
00177         }
00178 
00179         // ------------------------------ GroundEffect Class -----------------------------------------
00180 
00188         void GroundEffect::CorrectedGroundCoeffCalculation()
00189         {
00190                 double maxDp = 30 * (_zEqS + _zEqR);
00191                 if (_dp <= maxDp && maxDp != 0)
00192                 {
00193                         _correctedGpath = _dp * _Gpath / maxDp + (1 - (_dp / maxDp)) * _Gsource ;
00194                 }
00195                 else
00196                 {
00197                         _correctedGpath = _Gpath;
00198                 }
00199         }
00200     
00211         double GroundEffect::WparamCalculation(double Gw)
00212         {
00213                 double dividende = pow(_freq,1.5) * pow(Gw,2.6) + 1.3 * 1000 * pow(_freq,0.75) * pow(Gw,1.3) + 1.16 * 1000000;
00214                 double wParam;
00215                 if (dividende != 0)
00216                 {
00217                         wParam = 0.0185 * (pow(_freq,2.5) * pow(Gw,2.6))/ dividende;
00218                 }
00219                 else
00220                 {
00221                         printf("Erreur : division par 0 dans le calcul de l'effet de sol \n");
00222                         throw new int(ERRDivZero);
00223                 }
00224 
00225                 return wParam;
00226         }
00227 
00236         void GroundEffect::CfCalculation()
00237         {
00238                 double wdpH, wdpF;
00239                 switch(_groundCalculationType)
00240                 {
00241                 case Asol:
00242                         wdpH = WparamCalculation(_correctedGpath) * _dp;
00243                         wdpF = WparamCalculation(_Gpath) * _dp;
00244                         break;
00245                 case DeltaSol_SO:
00246                         wdpH = WparamCalculation(_correctedGpath) * _dp;
00247                         wdpF = WparamCalculation(_Gpath) * _dp;
00248                         break;
00249                 case DeltaSol_OR:
00250                         wdpH = WparamCalculation(_Gpath) * _dp;
00251                         wdpF = wdpH;
00252                         break;
00253                 default:
00254                         return;
00255                 }
00256  
00257                 if ((1 + wdpH) != 0)
00258                 {
00259                         _cfH = _dp * (1 + 3 * wdpH * exp(-1 * sqrt(wdpH))) / (1 + wdpH);
00260                 }
00261                 else
00262                 {
00263                         printf("Erreur : division par 0 dans le calcul de l'effet de sol \n");
00264                         throw new int(ERRDivZero);
00265                 }
00266 
00267                 if ((1 + wdpF) != 0)
00268                 {
00269                         _cfF = _dp * (1 + 3 * wdpF * exp(-1 * sqrt(wdpF))) / (1 + wdpF);
00270                 }
00271                 else
00272                 {
00273                         printf("Erreur : division par 0 dans le calcul de l'effet de sol \n");
00274                         throw new int(ERRDivZero);
00275                 }
00276         }
00277 
00285         void GroundEffect::KfreqCalculation()
00286         {
00287                 _kFreq = 2 * PI * _freq / cSound;
00288         }
00289 
00306         double GroundEffect::AttenuationCalculation(double z1, double z2, double cf)
00307         {
00308                 double attSol = 0;
00309                 if (_kFreq != 0 && _dp != 0)
00310                 {
00311                         double CfOnK = cf / _kFreq;
00312                         if (CfOnK >= 0)
00313                         {
00314                                 double sqrt2CfOnK = sqrt(2 * CfOnK);
00315                                 double mult1 = 4 * _kFreq * _kFreq / (_dp * _dp);
00316                                 double mult2 = (z1 * z1 - sqrt2CfOnK * z1 + CfOnK);
00317                                 double mult3 = (z2 * z2 - sqrt2CfOnK * z2 + CfOnK);
00318 
00319                                 attSol = -10 * log10(mult1 * mult2 * mult3);
00320                         }
00321                         else
00322                         {
00323                                 printf("Erreur : racine carrée d'un nombre négatif dans le calcul de l'effet de sol \n");
00324                                 throw new int(ERRSqrtNegative);
00325                         }
00326                 }
00327                 else
00328                 {
00329                         printf("Erreur : division par 0 dans le calcul de l'effet de sol \n");
00330                         throw new int(ERRDivZero);
00331                 }
00332 
00333                 return attSol;
00334         }
00335 
00344         double GroundEffect::AttenuationCalculationH()
00345         {
00346                 double attSolH = -3; // default value if Gpath = 0
00347 
00348                 if (_Gpath != 0)
00349                 {
00350                         attSolH = AttenuationCalculation(_zEqS, _zEqR, _cfH);
00351 
00352                         // Checks minimal value :
00353                         double attSolHmin; // formula (26)
00354                         switch(_groundCalculationType)
00355                         {
00356                         case Asol:
00357                                 attSolHmin = -3 * (1 - _correctedGpath);
00358                                 break;
00359                         case DeltaSol_SO:
00360                                 attSolHmin = -3 * (1 - _correctedGpath);
00361                                 break;
00362                         case DeltaSol_OR:
00363                                 attSolHmin = -3 * (1 - _Gpath);
00364                                 break;
00365                         }
00366 
00367                         if (attSolH < attSolHmin)
00368                         {
00369                                 attSolH = attSolHmin;
00370                         }
00371                 }
00372 
00373                 _AsolH = attSolH; 
00374                 return attSolH;
00375         }
00376 
00387         double GroundEffect::AttenuationCalculationF()
00388         {
00389          // modification dvm 20.04.2011
00390          // in case Gpath=0, we have AF = AFmin as decided at the meeting LRS/CSTD in Paris
00391                 double zEqSum = _zEqS + _zEqR;
00392 
00393                 double attSolFmin = 0;
00394                 switch(_groundCalculationType)
00395                 {
00396                 case Asol:
00397                         attSolFmin = -3 * (1 - _correctedGpath);
00398                         break;
00399                 case DeltaSol_SO:
00400                         attSolFmin = -3 * (1 - _correctedGpath);
00401                         break;
00402                 case DeltaSol_OR:
00403                         attSolFmin = -3 * (1 - _Gpath);
00404                         break;
00405                 }
00406 
00407                 double maxDp = 30 * zEqSum;
00408                 if (_dp > maxDp && _dp != 0)
00409                 {
00410                         attSolFmin *= (1 + 2 *(1 - maxDp / _dp));
00411                 }
00412 
00413                 double attSolF = attSolFmin ; // default value if Gpath = 0
00414 
00415                 if (_Gpath != 0)
00416                 {
00417                         if ((zEqSum) != 0)
00418                         {
00419                                 double a0 = 0.0002;
00420                                 double dp2by2 = _dp * _dp / 2; 
00421                                 double delta_zS = a0 * pow((_zEqS/(zEqSum)), 2) * dp2by2;
00422                                 double delta_zR = a0 * pow((_zEqR/(zEqSum)), 2) * dp2by2;
00423                                 double delta_zT = 0.006 * _dp / (zEqSum);
00424                                 _zEqS_f = _zEqS + delta_zS + delta_zT;
00425                                 _zEqR_f = _zEqR + delta_zR + delta_zT;
00426                                 attSolF = AttenuationCalculation(_zEqS_f, _zEqR_f, _cfF);
00427 
00428                                 // Checks minimal value :
00429 
00430                                 if (attSolF < attSolFmin)
00431                                 {
00432                                         attSolF = attSolFmin;
00433                                 }
00434 
00435                         }
00436                         else
00437                         {
00438                                 printf("Erreur : division par 0 dans le calcul de l'effet de sol \n");
00439                                 throw new int(ERRDivZero);
00440                         }
00441                 }
00442 
00443                 _AsolF = attSolF;
00444                 return attSolF;
00445         }
00446 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Defines