NMPB08 Shared Libraries
|
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 }