NMPB08 Shared Libraries
|
00001 00008 #include "Embankment.h" 00009 #include <math.h> 00010 #include <stdio.h> 00011 00012 namespace EmbankmentNMPB 00013 { 00023 Embankment::Embankment(vector<ProfilePointNMPB*> terrainItems, double freq) 00024 { 00025 _embankmentAttenuation = 0; 00026 _source = new Position2D; 00027 _receiver = new Position2D; 00028 _platformEnd = new Position2D; 00029 _O1 = new Position2D; 00030 _O2 = new Position2D; 00031 _sourceImage = new Position2D; 00032 00033 // Fill and check data 00034 if(freq <= 0) 00035 { 00036 printf (".ERROR: frequency = 0 \n") ; 00037 throw new int( ERRFrequency ); 00038 } 00039 00040 if (!DataInPerpendicularPlane(terrainItems)) 00041 { 00042 return; 00043 } 00044 00045 if (!CheckData(terrainItems)) 00046 { 00047 return; 00048 } 00049 00050 // source image S'' 00051 double aSlope = 0; 00052 double bSlope = 0; 00053 bool verticalSlope = LineCoeff(_O1, _O2, aSlope, bSlope); 00054 SourceImage(aSlope, bSlope, verticalSlope); 00055 if(_receiver->d == _sourceImage->d && _receiver->z == _sourceImage->z) 00056 { 00057 return; 00058 } 00059 00060 // Intersection between (S''R) and (O1O2) 00061 double aSimageR = 0; 00062 double bSimageR = 0; 00063 bool verticalSimageR = LineCoeff(_sourceImage, _receiver, aSimageR, bSimageR); 00064 Position2D intersect; 00065 if(!LinesIntersection(aSlope, bSlope, verticalSlope, aSimageR, bSimageR, verticalSimageR, &intersect)) 00066 { 00067 return; 00068 } 00069 00070 // Attenuation calculation 00071 EmbankmentCalculation(&intersect, freq); 00072 } 00073 00081 void Embankment::EmbankmentCalculation(Position2D const* intersect, double freq) 00082 { 00083 // Half-length e 00084 // § 9.3.5 formula (28) p.43 00085 double angleTheta = ThetaAngle(); 00086 double cosAngleTheta = cos(angleTheta); 00087 if(cosAngleTheta == 0) 00088 { 00089 return; 00090 } 00091 double SimageI = distance2D(_sourceImage, intersect); 00092 double IR = distance2D(intersect, _receiver); 00093 if(SimageI == 0 && IR == 0) 00094 { 00095 return; 00096 } 00097 00098 double halfLength = HalfLength(cosAngleTheta, freq, SimageI, IR); 00099 00100 // I x-coordinate on (O1O2) : 00101 double xI = distance2D(_O1, intersect); 00102 double prodScal = (_O2->d - _O1->d) * (intersect->d - _O1->d) + (_O2->z - _O1->z) * (intersect->z - _O1->z); 00103 if (prodScal < 0) 00104 { 00105 xI = -1 * xI; 00106 } 00107 00108 // Epsilon (formula (29) p. 44) 00109 double lT = distance2D(_O1, _O2); 00110 double epsilon = 0; 00111 if(xI > (-1 * halfLength) && xI < (lT + halfLength)) 00112 { 00113 double val1 = lT / (2 * halfLength); 00114 double val2 = (xI - halfLength) / (2 * halfLength); 00115 double val3 = (xI + halfLength) / (2 * halfLength); 00116 if (val1 <= val3) //(val1 <= val2) //// 00117 { 00118 epsilon = val1; 00119 } 00120 else 00121 { 00122 epsilon = val3; 00123 } 00124 00125 if(val2 > 0) 00126 { 00127 epsilon -= val2; 00128 } 00129 } 00130 00131 // Attenuation (formula (30) p.44) 00132 _embankmentAttenuation = -1.5 * epsilon * (2 - _Gembankment); 00133 } 00134 00144 void Embankment::FillPosInPerpendicularPlane(Position2D* posToFill, Position2D const* posToRead, double cosTheta) 00145 { 00146 if(posToFill == NULL || posToRead == NULL || cosTheta > 1 || cosTheta < -1) 00147 { 00148 return; 00149 } 00150 00151 posToFill->d = posToRead->d * cosTheta; 00152 posToFill->z = posToRead->z; 00153 } 00154 00165 bool Embankment::DataInPerpendicularPlane(vector<ProfilePointNMPB*> terrainItems) 00166 { 00167 int nbItems = (int) terrainItems.size(); 00168 if (nbItems < 3) 00169 { 00170 return false; 00171 } 00172 00173 ProfilePointNMPB* source = terrainItems[0]; 00174 if(source->ext == NULL || source->ext->type != ETRoadSource_NMPB) 00175 { 00176 return false; 00177 } 00178 00179 double cosTheta = source->ext->cosTheta; 00180 00181 if(_source) delete _source; 00182 _source = new Position2D; 00183 Position2D *posSource = source->position2D_withHeight(); 00184 FillPosInPerpendicularPlane(_source, posSource, cosTheta); 00185 delete posSource; 00186 bool embankmentFound = false; 00187 bool platformFound = false; 00188 for (int i = 1 ; i < nbItems ; i ++) 00189 { 00190 ProfilePointNMPB* terrainItem = terrainItems[i]; 00191 00192 if(terrainItem->ext != NULL) 00193 { 00194 if(terrainItem->ext->type == ETPlatform_NMPB) 00195 { 00196 // platform only before embankment 00197 if(!embankmentFound) 00198 { 00199 if(_platformEnd) delete _platformEnd; 00200 _platformEnd = new Position2D; 00201 FillPosInPerpendicularPlane(_platformEnd, terrainItem->position2D, cosTheta); 00202 } 00203 } 00204 else if(terrainItem->ext->type == ETEmbankment_NMPB) 00205 { 00206 if(embankmentFound == true) 00207 { 00208 // only one embankment is accepted 00209 printf (".ERROR: 2 embankments detected \n") ; 00210 throw new int( ERREmbankment ); 00211 } 00212 if(_O1) delete _O1; 00213 _O1 = new Position2D; 00214 FillPosInPerpendicularPlane(_O1, terrainItems[i-1]->position2D, cosTheta); 00215 platformFound = true; 00216 00217 if(_O2) delete _O2; 00218 _O2 = new Position2D; 00219 FillPosInPerpendicularPlane(_O2, terrainItem->position2D, cosTheta); 00220 _Gembankment = terrainItem->impedance; 00221 embankmentFound = true; 00222 } 00223 } 00224 } 00225 00226 if(!embankmentFound || !platformFound) 00227 { 00228 return false; 00229 } 00230 00231 if(_receiver) delete _receiver; 00232 _receiver = new Position2D; 00233 Position2D *lastPos = terrainItems[nbItems - 1]->position2D_withHeight(); 00234 FillPosInPerpendicularPlane(_receiver, lastPos, cosTheta); 00235 delete lastPos; 00236 00237 // Calculation of _angleAlpha 00238 if(_O1->d == _O2->d && _O1->z == _O2->z) 00239 { 00240 return false; 00241 } 00242 _angleAlpha = atan2(_O2->z - _O1->z, _O2->d - _O1->d); 00243 00244 return true; 00245 } 00246 00256 bool Embankment::CheckData(vector<ProfilePointNMPB*> terrainItems) 00257 { 00258 // check slope : 00259 if (_angleAlpha < PI / 12 || _angleAlpha > PI / 4) 00260 { 00261 return false; 00262 } 00263 00264 // the foot of the bank and the receiver are in direct view : 00265 double aCoef = 0; 00266 double bCoef = 0; 00267 bool vert = LineCoeff(_O1, _receiver, aCoef, bCoef); 00268 if(!vert) 00269 { 00270 for (int i = 0 ; i < (int) terrainItems.size() ; i++) 00271 { 00272 ProfilePointNMPB* elem = terrainItems[i]; 00273 Position2D* posElem = elem->position2D_withHeight(); 00274 if(posElem->d > _O1->d && posElem->d < _receiver->d) 00275 { 00276 if(posElem->z >= aCoef * posElem->d + bCoef) 00277 { 00278 delete posElem; 00279 return false; 00280 } 00281 } 00282 delete posElem; 00283 } 00284 } 00285 00286 // the foot of the bank is at a maximum distance of 10 metres from the edge of the closest platform 00287 double distToPlatform = distance2D(_platformEnd, _O1); 00288 if(distToPlatform > 10) 00289 { 00290 return false; 00291 } 00292 00293 // the foot of the bank is lower than or at the same height as the edge of the closest platform 00294 if(_O1->z > _O2->z) 00295 { 00296 return false; 00297 } 00298 00299 return true; 00300 } 00301 00316 bool Embankment::LineCoeff(Position2D const* pos1, Position2D const* pos2, double &aCoef, double &bCoef) 00317 { 00318 double dividende = pos2->d - pos1->d ; 00319 if(dividende != 0) 00320 { 00321 // z = a*d + b 00322 aCoef = (pos2->z - pos1->z) / dividende; 00323 bCoef = pos1->z - aCoef * pos1->d; 00324 return false; 00325 } 00326 else 00327 { 00328 // d = b 00329 aCoef = 0; 00330 bCoef = pos1->d; 00331 return true; 00332 } 00333 } 00334 00344 void Embankment::SourceImage(double aCoef, double bCoef, bool verticalSlope) 00345 { 00346 if(_sourceImage) delete _sourceImage; 00347 _sourceImage = new Position2D; 00348 if(verticalSlope) 00349 { 00350 // d = bCoef 00351 _sourceImage->d = 2 * bCoef - _source->d; 00352 _sourceImage->z = _source->z; 00353 00354 } 00355 else 00356 { 00357 // z = aCoef * d + bCoef 00358 double dividende = aCoef * aCoef + 1; 00359 _sourceImage->d = (-1 * aCoef * aCoef * _source->d + 2 * aCoef * _source->z + _source->d - 2 * aCoef * bCoef) / dividende; 00360 _sourceImage->z = (aCoef * aCoef * _source->z + 2 * aCoef * _source->d - _source->z + 2 * bCoef) / dividende; 00361 } 00362 } 00363 00383 bool Embankment::LinesIntersection(double a1, double b1, bool vert1, double a2, double b2, bool vert2, Position2D* Ipoint) 00384 { 00385 if(Ipoint == NULL) 00386 { 00387 return false; 00388 } 00389 if (vert1) 00390 { 00391 if(vert2) 00392 { 00393 // parallel lines : 0 or infinity intersections 00394 return false; 00395 } 00396 else 00397 { 00398 Ipoint->d = b1; 00399 Ipoint->z = a2 * b1 + b2; 00400 } 00401 } 00402 else 00403 { 00404 if(vert2) 00405 { 00406 Ipoint->d = b2; 00407 Ipoint->z = a1 * b2 + b1; 00408 } 00409 else 00410 { 00411 double dividende = a2 - a1; 00412 if(dividende == 0) 00413 { 00414 // parallel lines : 0 or infinity intersections 00415 return false; 00416 } 00417 00418 Ipoint->d = (b1 - b2) / dividende; 00419 Ipoint->z = a1 * Ipoint->d + b1; 00420 } 00421 } 00422 00423 return true; 00424 } 00425 00431 double Embankment::ThetaAngle() 00432 { 00433 // NB : _angleAlpha already calculated when check data 00434 00435 // the (O1x,S''R) angle 00436 double angleBeta = atan2(_receiver->z - _sourceImage->z, _receiver->d - _sourceImage->d); 00437 00438 return _angleAlpha + angleBeta; 00439 } 00440 00456 double Embankment::HalfLength(double cosAngleTheta, double freq, double SimageI, double IR) 00457 { 00458 // Already checked : freq > 0, cosAngleTheta != 0, SimageI + IR != 0 00459 double wavelength = cSound / freq; 00460 double val1 = (wavelength / 2) * (SimageI * IR)/(SimageI + IR) + (wavelength * wavelength / 16); 00461 double halfLength = sqrt(val1) / cosAngleTheta; 00462 00463 return halfLength; 00464 } 00465 00466 } 00467