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