TYAcousticModel.cpp

Go to the documentation of this file.

/*
 * Copyright (C) <2012-2014> <EDF-R&D> <FRANCE>
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 * See the GNU General Public License for more details.
 * You should have received a copy of the GNU General Public License along
 * with this program; if not, write to the Free Software Foundation, Inc.,
 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
 */
 
#include <deque>
#include <list>
#include <cmath>
#include <algorithm>
#include "Tympan/core/defines.h"
#include "Tympan/models/solver/config.h"
#include "Tympan/solvers/9613Solver/TYTrajet.h"
#include "Tympan/solvers/9613Solver/TYSolver.h"
#include "Tympan/models/common/mathlib.h"
#include "Tympan/geometric_methods/AcousticRaytracer/Geometry/Scene.h"
#include "TYAcousticModel.h"
#include "TYSolver.h"
#include "TYSolverHelper.h"
#include "TYFaceSelector.h"
 
TYAcousticModel::TYAcousticModel(TYSolver& solver)
    : _useSol(true), _useReflex(false), _propaCond(0), _interference(false), _paramH(10.0), _solver(solver)
{
    _absoNulle = OSpectreOctave(1.0);
    _absoNulle.setType(SPECTRE_TYPE_ABSO); // Spectre d'absorption
}
 
TYAcousticModel::~TYAcousticModel() {}
 
void TYAcousticModel::init()
{
    tympan::LPSolverConfiguration config = tympan::SolverConfiguration::get();
    // Calcul avec sol reel
    _useSol = config->UseRealGround;
    // Calcul avec reflexion sur les parois verticales
    _useReflex = config->UseReflection;
    // Calcul en conditions favorables
    _propaCond = config->PropaConditions;
    // Definit l'atmosphere courante du site
    double pression = config->AtmosPressure;
    double temperature = config->AtmosTemperature;
    double hygrometrie = config->AtmosHygrometry;
 
    pSolverAtmos =
        std::unique_ptr<AtmosphericConditions>(new AtmosphericConditions(pression, temperature, hygrometrie));
    // Calcul avec interference
    _interference = config->ModSummation;
 
    // Compute wave length
    double c_9613_2 = 340.0;
    _lambda = OSpectreOctave::getLambda(c_9613_2);
    // Coefficient multiplicateur pour le calcul des reflexions supplementaires en condition favorable
    _paramH = config->H1parameter;
}
 
void TYAcousticModel::compute(const std::deque<TYSIntersection>& tabIntersect, TYTrajet& trajet,
                              TabPoint3D& ptsTop, TabPoint3D& ptsLeft, TabPoint3D& ptsRight)
{
    bool vertical = true, horizontal = false;
    bool left = true, right = false;
 
    // Construction du rayon SR
    OSegment3D rayon;
    trajet.getPtSetPtRfromOSeg3D(rayon);
    bool conditionFav = false;
 
    // Calcul des conditions de propagation suivant la direction du vent
    tympan::LPSolverConfiguration config = tympan::SolverConfiguration::get();
    assert(config->DSWindDirection >= 0 && config->DSWindDirection <= 360);
 
    double windRadian = DEGTORAD(config->DSWindDirection);
    OVector3D windDirection = OVector3D(-sin(windRadian), -cos(windRadian), 0);
    OVector3D propaDirection = rayon.toVector3D();
    propaDirection._z = 0;
    double angle =
        RADTODEG(acos(windDirection.dot(propaDirection) /
                      (windDirection.norme() * propaDirection.norme()))); // Angle always between 0-180
    assert(180 >= angle >= 0);
    assert(180 >= config->AngleFavorable >= 0);
 
    if (angle <= config->AngleFavorable)
    {
        conditionFav = true;
    }
    else
    {
        conditionFav = false;
    }
 
    // Recuperation de la source
    tympan::AcousticSource& source = trajet.asrc;
 
    // Distance de la source au recepteur
    double distance = trajet.getDistance();
 
    TYTabChemin& tabChemins = trajet.getChemins();
 
    // Calcul du chemin direct
    computeCheminSansEcran(tabIntersect, rayon, source, tabChemins, distance, conditionFav);
 
    if (ptsTop.size() > 1 || ptsLeft.size() > 1 || ptsRight.size() > 1)
    {
        // Calcul des parcours lateraux
        // 1. Vertical
        computeCheminsAvecEcran(rayon, source, ptsTop, vertical, tabChemins, distance, right);
 
        // 2. Horizontal gauche
        computeCheminsAvecEcran(rayon, source, ptsLeft, horizontal, tabChemins, distance, left);
 
        // 3. Horizontal droite
        computeCheminsAvecEcran(rayon, source, ptsRight, horizontal, tabChemins, distance, right);
    }
 
    // Calcul des reflexions si necessaire
    computeCheminReflexion(tabIntersect, rayon, source, tabChemins, distance);
 
    // Calcul la pression cumulee de tous les chemins au point de reception du trajet
    solve(trajet);
 
    // Le calcul est fini pour ce trajet, on peut effacer les tableaux des chemins
    tabChemins.clear();
}
 
void TYAcousticModel::computeCheminAPlat(const OSegment3D& rayon, const tympan::AcousticSource& source,
                                         TYTabChemin& TabChemins, double distance) const
{
    TYTabEtape tabEtapes;
 
    // Calcul de la pente moyenne sur le trajet source-recepteur
    OSegment3D penteMoyenne;
    meanSlope(rayon, penteMoyenne);
 
    // Etape directe Source-Recepteur
    TYEtape etape1;
    TYChemin chemin1;
 
    etape1._pt = rayon._ptA;
    etape1._type = TYSOURCE;
    etape1._Absorption =
        source.directivity->lwAdjustment(OVector3D(rayon._ptA, rayon._ptB), rayon.longueur());
 
    chemin1.setType(TYTypeChemin::CHEMIN_DIRECT);
 
    tabEtapes.push_back(etape1);   // Ajout de l'etape directe
    chemin1.setLongueur(distance); // Dans ce cas, la longueur = la distance source/recepteur
    chemin1.setDistance(distance);
    chemin1.calcAttenuation(tabEtapes, *pSolverAtmos);
 
    TabChemins.push_back(chemin1);
    tabEtapes.clear(); // Vide le tableau des etapes
 
    // Liaison avec reflexion
    //              1. Calcul du point de reflexion
 
    // Ajout du chemin reflechi
    TYEtape etape2;
    TYChemin chemin2;
    chemin2.setType(TYTypeChemin::CHEMIN_SOL);
 
    etape2.setPoint(rayon._ptA);
    etape2._type = TYSOURCE;
 
    OPoint3D ptSym;
    int symOK = 0;
    if (penteMoyenne.longueur() > 0) // Si la pente moyenne est definie, on prend le point symetrique
    {
        symOK = penteMoyenne.symetrieOf(rayon._ptA, ptSym);
    }
 
    if (symOK == 0) // Sinon on prend une simple symetrie par rapport a z
    {
        ptSym = rayon._ptA;
        ptSym._z = 2 * penteMoyenne._ptA._z - ptSym._z;
    }
 
    // Calcul du point de reflexion
    OPoint3D ptReflex;
    penteMoyenne.intersects(OSegment3D(ptSym, rayon._ptB), ptReflex, TYSEUILCONFONDUS);
 
    //              2. Etape avant la reflexion
    OSegment3D seg1(rayon._ptA, ptReflex); // On part de la source vers le point de reflexion
    double rr = seg1.longueur();
 
    // Directivite de la source
    etape2._Absorption = source.directivity->lwAdjustment(OVector3D(seg1._ptA, seg1._ptB), seg1.longueur());
 
    tabEtapes.push_back(etape2); // Ajout de l'etape avant reflexion
 
    //              3. Etape apres la reflexion
 
    TYEtape etape3;
    etape3._pt = ptReflex;
    etape3._type = TYREFLEXIONSOL;
 
    OSegment3D seg2 = OSegment3D(ptReflex, rayon._ptB); // Segment Point de reflexion->Point de reception
    rr += seg2.longueur();                              // Longueur parcourue sur le trajet reflechi
 
    if (_useSol)
    {
        etape3._Absorption = getReflexionSpectrumAt(seg1, rr, penteMoyenne, source);
    }
    else // Sol totalement reflechissant
    {
        etape3._Absorption = _absoNulle;
    }
 
    tabEtapes.push_back(etape3); // Ajout de l'etape apres reflexion
 
    chemin2.setLongueur(rr);
    chemin2.setDistance(distance);
    chemin2.calcAttenuation(tabEtapes, *pSolverAtmos);
    TabChemins.push_back(chemin2);
 
    tabEtapes.clear(); // Vide le tableau des etapes
}
 
bool TYAcousticModel::computeCheminsAvecEcran(const OSegment3D& ray, const tympan::AcousticSource& source,
                                              const TabPoint3D& pts, const bool vertical,
                                              TYTabChemin& tabPaths, double distance, const bool left) const
{
    /* ============================================================================================================
        In 9613, no reflexion on the ground is computed
     ==============================================================================================================*/
    if (pts.size() <= 1)
    {
        tabPaths[0].setAttenuationBarWhenNoPath(vertical, left);
        return false;
    }
 
    double dss{0.0}; // Length between source and first edge of diffraction
    double dsr{0.0}; // Length between last edge of diffraction and receptor
 
    OPoint3D firstPt(pts[1]);
    OPoint3D lastPt(pts[pts.size() - 1]);
 
    TYTabEtape tabSteps;
    double pathLength = 0.0;
 
    /*--- BEFORE OBSTACLE ---*/
 
    TYTabEtape steps;
    OSegment3D curSegment(ray._ptA, firstPt);
    double tempLength = curSegment.longueur();
 
    bool bPathOk = addStep(ray._ptA, firstPt, source, true, steps); // Add step before obstacle
 
    // If a problem has occurred, stop path creation
    if (!bPathOk)
    {
        return true;
    }
 
    tabSteps.push_back(steps[0]); // Add source step to table of steps
    pathLength += tempLength;
    dss = tempLength;
 
    steps.clear(); // Clear steps content
 
    /*--- BYPASS OF THE OBSTACLE ---*/
 
    double width = 0.0;
    TYEtape step;
 
    for (unsigned int i = 1; i < pts.size() - 1; i++)
    {
        width += (OSegment3D(pts[i], pts[i + 1])).longueur();
 
        step._pt = pts[i];
        step._type = TYDIFFRACTION;
        step._Absorption = _absoNulle;
 
        tabSteps.push_back(step);
    }
 
    pathLength += width;
 
    /*--- AFTER OBSTACLE ---*/
    curSegment = OSegment3D(lastPt, ray._ptB);
    tempLength = curSegment.longueur();
 
    addStep(lastPt, ray._ptB, source, false, steps);
 
    tabSteps.push_back(steps[0]);
    pathLength += tempLength;
    dsr = tempLength;
 
    steps.clear();
 
    /*--- COMPUTE SCREEN EFFECT ATTENUATION ON THE CURRENT PATH ---*/
 
    OSpectreOctave Dz;
 
    step._pt = ray._ptB;
    step._Absorption = _absoNulle;
 
    Dz = calculAttDiffraction(ray, pathLength, dss, dsr, width, vertical);
    step._Attenuation = Dz;
    tabSteps.push_back(step);
 
    /*--- ADD PATH TO the table of paths ---*/
 
    TYChemin path;
    path.setType(TYTypeChemin::CHEMIN_ECRAN);
    path.setDistance(distance);
    path.setLongueur(pathLength);
 
    tabPaths.push_back(path);
 
    // Compute barrier attenuation
    tabPaths[0].computeBarAttenuation(Dz, vertical, left);
 
    // Build equivalent path for rays
    path.build_eq_path(tabSteps);
 
    tabSteps.clear();
    steps.clear();
 
    return true;
}
 
bool TYAcousticModel::addStep(const OPoint3D& ptStart, const OPoint3D& ptEnd,
                              const tympan::AcousticSource& source, const bool& fromSource,
                              TYTabEtape& steps) const
{
    bool res = true;
 
    TYEtape curStep;
 
    // === BUILD DIRECT TRIP ptStart-ptEnd
    curStep._pt = ptStart;
 
    if (fromSource) // If we start from source, its directivity is considered
    {
        curStep._type = TYSOURCE;
        curStep._Absorption =
            source.directivity->lwAdjustment(OVector3D(ptStart, ptEnd), ptStart.distFrom(ptEnd));
    }
    else
    {
        curStep._type = TYDIFFRACTION;
        curStep._Absorption = _absoNulle;
    }
 
    steps.push_back(curStep);
 
    return res;
}
 
void TYAcousticModel::computeCheminSansEcran(const std::deque<TYSIntersection>& tabIntersect,
                                             const OSegment3D& ray, const tympan::AcousticSource& source,
                                             TYTabChemin& TabChemin, double distance, bool conditionFav) const
{
    /*
        COMPUTATION FOR A ROUTE WITHOUT OBSTACLE CONSISTS OF ONE DIRECT PATH
    */
    TYTabEtape tabEtapes;
 
    // Compute mean slope on source-receptor route
    OSegment3D segMeanSlope;
    meanSlope(ray, segMeanSlope);
 
    OPoint3D S2D{ray._ptA._x, ray._ptA._y, 0.0};
    OPoint3D R2D{ray._ptB._x, ray._ptB._y, 0.0};
    OSegment3D SR2D{S2D, R2D};
    double dp = SR2D.longueur();
 
    TYTabEtape Etapes;
    double Gs{0.0}, Gm{0.0}, Gr{0.0};
    double hs{0.0}, hr{0.0};
 
    computeSegmentEdgesHeights(hs, hr, segMeanSlope, ray);
 
    addStep(ray._ptA, ray._ptB, source, true, Etapes);
    getGroundfactors(tabIntersect, SR2D, hs, hr, Gs, Gm, Gr);
 
    // Add direct path
    TYChemin chemin;
    chemin.setType(TYTypeChemin::CHEMIN_DIRECT);
    tabEtapes.push_back(Etapes[0]); // Add direct step
    chemin.setLongueur(distance);   // In this case, lenght = source/receptor distance
    chemin.setDistance(distance);
    chemin.calcAttenuation(tabEtapes, *pSolverAtmos, dp, hs, hr, Gs, Gm, Gr);
    TabChemin.push_back(chemin); // Add path in array of paths
 
    tabEtapes.clear(); // Empty array of steps
}
 
bool TYAcousticModel::getGroundfactors(const std::deque<TYSIntersection>& tabIntersect,
                                       const OSegment3D& ray2D, double hs, double hr, double& Gs, double& Gm,
                                       double& Gr) const
{
    double heightRatio = 30.0;
 
    bool res = true;
 
    // === CONSTRUCTION OF DIRECT ROUTE ABOVE HORIZONTAL PLANE
    OPoint3D ptStartProj = ray2D._ptA; // PtDebut projected on horizontal plane
    OPoint3D ptEndProj = ray2D._ptB;   // PtFin projected on horizontal plane
 
    // === COMPUTE SOURCE, MIDDLE AND RECEPTOR ZONES
    OVector3D SR{ray2D._ptA, ray2D._ptB};
    double dp = SR.norme(); // Distance in meter, between source and receptor, projected on ground plan
    SR.normalize();
    OPoint3D ptGSrcZone;
    OPoint3D ptGRcpZone;
    OPoint3D ptGMidZone;
 
    // Source zone must not exceed receptor
    if (heightRatio * hs < dp)
    {
        ptGSrcZone = ptStartProj + SR * heightRatio * hs;
    }
    else
    {
        ptGSrcZone = ptStartProj + SR * dp;
    }
 
    // Receptor zone must not exceed source
    if (heightRatio * hr < dp)
    {
        ptGRcpZone = ptEndProj + (-1.0) * SR * heightRatio * hr;
    }
    else
    {
        ptGRcpZone = ptEndProj + (-1.0) * SR * dp;
    }
 
    // === COMPUTE GROUND FACTOR FOR EACH ZONE
    double GZone{0.0}, dpZone{0.0};
    computeGZone(ptStartProj, ptGSrcZone, GZone, dpZone, tabIntersect);
    if (dpZone != 0.0)
    {
        Gs = GZone / dpZone;
    }
    else
    {
        Gs = 0.5;
    }
    computeGZone(ptGRcpZone, ptEndProj, GZone, dpZone, tabIntersect);
    if (dpZone != 0.0)
    {
        Gr = GZone / dpZone;
    }
    else
    {
        Gr = 0.5;
    }
    computeGZone(ptGSrcZone, ptGRcpZone, GZone, dpZone, tabIntersect);
    if (dpZone != 0.0)
    {
        Gm = GZone / dpZone;
    }
    else
    {
        Gm = 0.5;
    }
 
    return res;
}
 
bool TYAcousticModel::getGroundfactors(const std::deque<TYSIntersection>& tabIntersectUpSegment,
                                       const std::deque<TYSIntersection>& tabIntersectDownSegment,
                                       const OSegment3D& SO2D, const OSegment3D& OR2D, double hs, double hr,
                                       double& Gs, double& Gm, double& Gr) const
{
    double heightRatio = 30.0;
    bool res = true;
 
    // Construction of points from segments parameters
    OPoint3D ptSrc2D = SO2D._ptA;
    OPoint3D ptO2D = SO2D._ptB;
    OPoint3D ptRcp2D = OR2D._ptB;
 
    // Computation of 2D distance and construction of unit vectors
    OVector3D SO{SO2D._ptA, SO2D._ptB};
    double dpSO =
        SO.norme(); // Distance in meter, between source and reflexion point, projected on ground plan
    SO.normalize();
 
    OVector3D OR{OR2D._ptA, OR2D._ptB};
    double dpOR =
        OR.norme(); // Distance in meter, between reflexion point and receptor, projected on ground plan
    OR.normalize();
 
    double dp = dpSO + dpOR; // Distance in meter, between source and receptor, projected on ground plan
 
    // === COMPUTE GROUND FACTOR FOR EACH ZONE
    OPoint3D ptGSrcZone;
    OPoint3D ptGRcpZone;
    OPoint3D ptGMidZone;
 
    bool bPtGSrcZoneInSO{false}; // True if ptGSrcZone in [SO] segment
    bool bPtGRcpZoneInOR{false}; // True if ptGRcpZone in [OR] segment
 
    // == COMPUTE GROUND FACTOR FOR SOURCE ZONE
    if (heightRatio * hs < dpSO)
    {
        // ptGSrcZone belongs to [SO]
        bPtGSrcZoneInSO = true;
        ptGSrcZone = ptSrc2D + (SO)*heightRatio * hs;
    }
    else if (heightRatio * hs < dp)
    {
        // ptGSrcZone belongs to [OR]
        ptGSrcZone = ptO2D + (OR) * (heightRatio * hs - dpSO);
    }
    else
    { // Source Zone must not exceed receptor
        ptGSrcZone = ptO2D + (OR)*dpOR;
    }
 
    // Compute Gs
    double GZone{0.0}, dpZone{0.0};
    if (bPtGSrcZoneInSO)
    {
        computeGZone(ptSrc2D, ptGSrcZone, GZone, dpZone, tabIntersectUpSegment);
    }
    else
    {
        double Gs1{0.0}, dp1{0.0}, Gs2{0.0}, dpOGs{0.0};
        computeGZone(ptSrc2D, ptO2D, Gs1, dp1, tabIntersectUpSegment);        // Gs1 = G(Src2D->ptO2D) / dpSO
        computeGZone(ptO2D, ptGSrcZone, Gs2, dpOGs, tabIntersectDownSegment); // Gs2 = G(ptO2D->Gs2D) / dpOGs
        GZone = Gs1 + Gs2;
        dpZone = dp1 + dpOGs;
    }
 
    if (dpZone != 0)
    {
        Gs = GZone / dpZone;
    }
    else
    {
        Gs = 0.5;
    }
 
    // == COMPUTE GROUND FACTOR FOR RECEPTOR ZONE
    if (heightRatio * hr < dpOR)
    {
        // ptGRcpZone belongs to [OR]
        bPtGRcpZoneInOR = true;
        ptGRcpZone = ptRcp2D + (-1.0) * OR * heightRatio * hr;
    }
    else if (heightRatio * hr < dp)
    {
        // ptGSrcZone belongs to [SO]
        ptGRcpZone = ptO2D + (-1.0) * SO * (heightRatio * hr - dpOR);
    }
    else
    { // Source Zone must not exceed receptor
        ptGRcpZone = ptO2D + (-1.0) * SO * dpSO;
    }
 
    // Compute Gr
    dpZone = 0.0;
    if (bPtGRcpZoneInOR)
    {
        computeGZone(ptGRcpZone, ptRcp2D, GZone, dpZone, tabIntersectDownSegment);
    }
    else
    {
        double Gr1{0.0}, dp2{0.0}, Gr2{0.0}, dpGrO{0.0};
        computeGZone(ptGRcpZone, ptO2D, Gr1, dpGrO, tabIntersectUpSegment); // Gr1 = G(Gr2D->ptO2D) / dpGrO
        computeGZone(ptO2D, ptRcp2D, Gr2, dp2, tabIntersectDownSegment);    // Gr2 = G(ptO2D->ptRcp2D) / dpOR
 
        GZone = Gr1 + Gr2;
        dpZone = dpGrO + dp2;
    }
 
    if (dpZone != 0)
    {
        Gr = GZone / dpZone;
    }
    else
    {
        Gr = 0.5;
    }
 
    // == COMPUTE GROUND FACTOR FOR MIDDLE ZONE
    // If GSrc and GRcp belong to [SO]
    if (bPtGSrcZoneInSO && !bPtGRcpZoneInOR)
    {
        computeGZone(ptGSrcZone, ptGRcpZone, GZone, dpZone, tabIntersectUpSegment);
    }
    // Else, if GSrc and GRcp belong to [OR]
    else if (!bPtGSrcZoneInSO && bPtGRcpZoneInOR)
    {
        computeGZone(ptGSrcZone, ptGRcpZone, GZone, dpZone, tabIntersectDownSegment);
    }
    // Else, if GSrc belongs to [SO], therefore GRcp belongs to [OR]
    else if (bPtGSrcZoneInSO)
    {
        double Gm1{0.0}, dpm1{0.0}, Gm2{0.0}, dpm2{0.0};
        computeGZone(ptGSrcZone, ptO2D, Gm1, dpm1, tabIntersectUpSegment);
        computeGZone(ptO2D, ptGRcpZone, Gm2, dpm2, tabIntersectDownSegment);
        GZone = Gm1 + Gm2;
        dpZone = dpm1 + dpm2;
    }
    // Else if GSrc belongs to [OR], therefore GRcp belongs to [SO]
    else
    {
        double Gm1{0.0}, dpm1{0.0}, Gm2{0.0}, dpm2{0.0};
        computeGZone(ptGRcpZone, ptO2D, Gm1, dpm1, tabIntersectUpSegment);
        computeGZone(ptO2D, ptGSrcZone, Gm2, dpm2, tabIntersectDownSegment);
        GZone = Gm1 + Gm2;
        dpZone = dpm1 + dpm2;
    }
 
    if (dpZone != 0)
    {
        Gm = GZone / dpZone;
    }
    else
    {
        Gm = 0.5;
    }
 
    return res;
}
 
OSpectreOctave TYAcousticModel::computeEffectiveBarAttenuation(const OSpectreOctave& Abar_top,
                                                               const OSpectreOctave& Abar_left,
                                                               const OSpectreOctave& Abar_right)
{
    OSpectreOctave Abar;
    for (unsigned int i = 0; i < TY_SPECTRE_OCT_NB_ELMT; i++)
    {
        Abar.getTabValReel()[i] = -10 * log10(pow(10, -0.1 * Abar_top.getTabValReel()[i]) +
                                              pow(10, -0.1 * Abar_left.getTabValReel()[i]) +
                                              pow(10, -0.1 * Abar_right.getTabValReel()[i]));
        if (Abar.getTabValReel()[i] < 0)
        {
            Abar.getTabValReel()[i] = 0;
        }
    }
    return Abar;
}
 
OPlan TYAcousticModel::buildMeanSlopePlan(const OSegment3D& penteMoyenne) const
{
    OPoint3D pt1{penteMoyenne._ptA};
    OPoint3D pt2{penteMoyenne._ptB};
    OPoint3D pt3{};
    pt3._z = pt1._z;
 
    if (pt2._x != pt1._x)
    {
        pt3._y = pt1._y + 1;
        pt3._x = (pt1._y - pt2._y) * (pt3._y - pt1._y) / (pt2._x - pt1._x) + (pt1._x);
    }
    else
    {
        if (pt1._y != pt2._y)
        {
            pt3._x = pt1._x + 1;
            pt3._y = (pt2._x - pt1._x) * (pt3._x - pt1._x) / (pt1._y - pt2._y) + (pt1._y);
        }
        else // pt1 and pt2 coincide, we shift pt2, building an horizontal plane
        {
            pt2._x = pt1._x + 1;
            pt2._y = pt1._y - 1;
            pt3._y = pt1._y + 1;
            pt3._x = (pt1._y - pt2._y) * (pt3._y - pt1._y) / (pt2._x - pt1._x) + (pt1._x);
        }
    }
 
    return OPlan{pt1, pt2, pt3};
}
 
bool TYAcousticModel::computeGZone(const OPoint3D& ptDebut, const OPoint3D& ptFin, double& GZone,
                                   double& dpZone, const std::deque<TYSIntersection>& tabIntersect) const
{
    bool ret = true;
    OSegment3D segZone(ptDebut, ptFin);
    OVector3D DF(ptDebut, ptFin);
 
    // Loop on intersections of Topography triangles with EV plane.
    // For each triangle, we search the intersection between the intersecting segment and the zone segment.
    // The intersecting segment has an homogeneous ground factor G.
    // It is used to compute a balanced ground factor over the zone segment.
    size_t nbTriangles = tabIntersect.size();
    double currentG = 0.0;
    GZone = 0.0;
 
    // Edges of current intersecting segment
    OPoint3D ptDebutResult;
    OPoint3D ptFinResult;
 
    for (unsigned int i = 0; i < nbTriangles; i++)
    {
        TYSIntersection inter = tabIntersect[i];
 
        // If triangle is not topography or does not intersect EV plane, then continue with following triangle
        if ((inter.isInfra) || !(inter.bIntersect[0]))
        {
            continue;
        }
 
        OSegment3D currentSeg = tabIntersect[i].segInter[0];
        currentG = tabIntersect[i].material->get_ISO9613_G();
 
        // We build AB segment the projection of the topography segment on the horizontal plane
        OPoint3D ptDebutCurrentProj{currentSeg._ptA._x, currentSeg._ptA._y, 0.0};
        OPoint3D ptFinCurrentProj{currentSeg._ptB._x, currentSeg._ptB._y, 0.0};
 
        OSegment3D segAB{ptDebutCurrentProj, ptFinCurrentProj};
        OVector3D AB{ptDebutCurrentProj, ptFinCurrentProj};
        // Orientate AB segment as zone segment DF
        if (AB.scalar(DF) <= 0)
        {
            segAB = segAB.swap();
        }
        OVector3D AD(segAB._ptA, segZone._ptA);
        OVector3D AF(segAB._ptA, segZone._ptB);
        OVector3D BD(segAB._ptB, segZone._ptA);
        OVector3D BF(segAB._ptB, segZone._ptB);
 
        bool intersect = true;
 
        // If A belongs to zone segment DF
        if (AD.scalar(AF) <= 0)
        {
            // then A is the starting edge of the intersecting segment
            ptDebutResult = segAB._ptA;
 
            // If B belongs to zone segment DF
            if (BD.scalar(BF) <= 0)
            {
                // then B is the ending edge of the intersecting segment
                ptFinResult = segAB._ptB;
            }
            else
            // else F is the ending edge of the intersecting segment
            {
                ptFinResult = segZone._ptB;
            }
        }
        else
        {
            // Else A does not belong to zone segment DF
            // If B does not belong to zone segment either
            if (BD.scalar(BF) >= 0)
            {
                // and if A and B are from each side of zone segment
                if (AD.scalar(BF) <= 0)
                {
                    // then intersecting segment is the zone segment DF itself
                    ptDebutResult = segZone._ptA;
                    ptFinResult = segZone._ptB;
                }
                else
                {
                    // Else A and B are on the same side of segment zone, so no intersection
                    intersect = false;
                }
            }
            else
            // Else B belong to segment zone DF but not A
            {
                ptDebutResult = segZone._ptA;
                ptFinResult = segAB._ptB;
            }
        }
 
        if (intersect)
        {
            OVector3D result{ptDebutResult, ptFinResult};
            GZone = GZone + result.norme() * currentG;
        }
    }
    OVector3D zone{ptDebut, ptFin};
    dpZone = zone.norme();
    return ret;
}
 
void TYAcousticModel::computeCheminReflexion(const std::deque<TYSIntersection>& tabIntersect,
                                             const OSegment3D& ray, const tympan::AcousticSource& source,
                                             TYTabChemin& TabChemins, double distance) const
{
    if (!_useReflex)
    {
        return;
    }
 
    OSegment3D segInter;
    OSegment3D rayonTmp;
    OPoint3D ptSym;
    OSpectreOctave SpectreAbso;
 
    OSegment3D seg;         // Image source -> receptor segment
    OSegment3D upwardSeg;   // Source -> reflexion point segment
    OSegment3D downwardSeg; // Reflexion point -> receptor segment
 
    OPoint3D pt; // Intersection (reflexion) point
 
    size_t nbFaces = tabIntersect.size();
 
    // For each face test reflexion
    for (unsigned int i = 0; i < nbFaces; i++)
    {
        TYSIntersection inter = tabIntersect[i];
 
        // If face cannot interact skip it
        if ((!inter.isInfra) || !(inter.bIntersect[1]))
        {
            continue;
        }
 
        segInter = inter.segInter[1];
 
        // Compute symmetric of A with respect to the segment
        segInter.symetrieOf(ray._ptA, ptSym); // We don't deal with this function return value
        seg._ptA = ptSym;
        seg._ptB = ray._ptB; // Image source -> receptor segment
 
        if (segInter.intersects(seg, pt, TYSEUILCONFONDUS))
        {
            // Construction of reflexion point -> source segment
            upwardSeg._ptA = ray._ptA;
            upwardSeg._ptB = pt;
            // Construction of reflexion point -> receptor segment
            downwardSeg._ptA = upwardSeg._ptB;
            downwardSeg._ptB = ray._ptB;
 
            bool intersect = false;
            size_t j = 0;
 
            // If we cross another face, which can be topography, the reflexion path is not taken into account
            while ((j < nbFaces) && (!intersect))
            {
                if (j == i)
                {
                    j++;
                    continue; // If face cannot interact skip it
                }
 
                segInter = tabIntersect[j].segInter[1];
 
                // We test whether segInter intersects upward segment or
                // downward segment in global plane.
                // Point pt is not use, we only care about testing intersection
                if ((segInter.intersects(upwardSeg, pt, TYSEUILCONFONDUS)) ||
                    (segInter.intersects(downwardSeg, pt, TYSEUILCONFONDUS)))
                {
                    // Intersection found, exit from the loop
                    intersect = true;
                    break;
                }
 
                j++;
            }
 
            // If reflected path is not intersected, reflexion can be computed
            if (!intersect)
            {
                SpectreAbso = dynamic_cast<tympan::AcousticBuildingMaterial*>(inter.material)->spectrum;
 
                TYTabEtape tabEtapes;
 
                double pathLength = upwardSeg.longueur() + downwardSeg.longueur();
 
                TYEtape Etape;
                // First step : from source to reflexion point
                Etape._pt = ray._ptA;
                Etape._type = TYSOURCE;
                Etape._Absorption = source.directivity->lwAdjustment(
                    OVector3D(upwardSeg._ptA, upwardSeg._ptB),
                    upwardSeg.longueur()); // Directivity factor toward receptor image
 
                tabEtapes.push_back(Etape);
 
                // Second step : from reflexion point to end of ray
                Etape._pt = downwardSeg._ptA;
                Etape._type = TYREFLEXION;
                Etape._Absorption = SpectreAbso;
 
                tabEtapes.push_back(Etape);
 
                // Compute mean slope on source-receptor route
                OSegment3D segMeanSlope;
                meanSlope(ray, segMeanSlope);
 
                OPoint3D S2D{upwardSeg._ptA._x, upwardSeg._ptA._y, 0.0};
                OPoint3D O2D{upwardSeg._ptB._x, upwardSeg._ptB._y, 0.0};
                OPoint3D R2D{downwardSeg._ptB._x, downwardSeg._ptB._y, 0.0};
                OSegment3D raySO{S2D, O2D};
                OSegment3D rayOR{O2D, R2D};
                double dp = raySO.longueur() + rayOR.longueur();
 
                TYTabEtape Etapes;
                double Gs{0.0}, Gm{0.0}, Gr{0.0};
                double hs{0.0}, hr{0.0};
 
                computeSegmentEdgesHeights(hs, hr, segMeanSlope, ray);
 
                // Compute intersecting segments for reflected ray
                std::deque<TYSIntersection> tabIntersectUpSegment, tabIntersectDownSegment;
                _solver.getFaceSelector()->selectFaces(tabIntersectUpSegment, upwardSeg, source.volume_id);
                _solver.getFaceSelector()->selectFaces(tabIntersectDownSegment, downwardSeg,
                                                       source.volume_id);
 
                // Compute ground factors for reflected ray
                getGroundfactors(tabIntersectUpSegment, tabIntersectDownSegment, raySO, rayOR, hs, hr, Gs, Gm,
                                 Gr);
 
                TYChemin Chemin;
                Chemin.setType(TYTypeChemin::CHEMIN_REFLEX);
                Chemin.setLongueur(pathLength);
                Chemin.setDistance(distance);
                Chemin.calcAttenuation(tabEtapes, *pSolverAtmos, dp, hs, hr, Gs, Gm, Gr);
 
                TabChemins.push_back(Chemin); // Put the reflected path in paths table
                tabEtapes.clear();
            }
        }
    }
}
 
OSpectreOctave TYAcousticModel::calculC3(const double& width) const
{
    // C3 = (1 + (5 * lambda / width)^2) / (1 / 3 + (5 * lambda / width)^2)
 
    OSpectreOctave C3 = OSpectreOctave::getEmptyLinSpectre();
    OSpectreOctave opLambda;
 
    if (width < 0.5)
    {
        C3.setDefaultValue(1.0);
    }
    else
    {
        const double oneThird = 1.0 / 3.0;
 
        opLambda = _lambda * (5.0 / width); // (5*lambda/e)
        opLambda = opLambda * opLambda;     // (5*lambda/e)^2
 
        C3 = opLambda + 1.0;              // 1 + (5*lambda/e)^2
        C3 = C3.div(opLambda + oneThird); // (1 + (5*lambda/e)^2) / (1/3 + (5*lambda/e)^2)
    }
 
    C3.setType(SPECTRE_TYPE_AUTRE); // Neither Attenuation, nor Absorption
 
    return C3;
}
 
OSpectreOctave TYAcousticModel::calculAttDiffraction(const OSegment3D& ray, const double& re,
                                                     const double& dss, const double& dsr,
                                                     const double& width, const bool& vertical) const
{
    double rd;
 
    OSpectreOctave s, Dz;
 
    OSpectreOctave C3 = calculC3(width); // Corrective factor linked with screen width
 
    double C2{20.0};
 
    rd = ray.longueur();
 
    double z = re - rd; // Path-length difference
    z = z <= 0 ? 0.0 : z;
 
    double Kmeteo = 1.0;
    if (z > 0.0 && vertical)
    {
        Kmeteo = exp(-(1.0 / 2000.0) * sqrt(dss * dsr * rd / (2 * z)));
    }
 
    // Attenuation brought by diffraction Dz = 10 * log [3 + (C2 / lambda) * C3 * delta * Kmeteo] dB
 
    s = _lambda.invMult(C2 * z); // =C2*z/lambda
    s = s * C3 * Kmeteo;         // C2*z*C3*Kmeteo/lambda
    s = s + 3.0;
    Dz = s.log() * 10.0;
    // If diffraction occurs in vertical plane (horizontal edge)
    // minimal and amaximal attenuations are limited respectively
    // to 0 dB and 20 or 25 dB (whether screen is thin or wide).
    if (vertical)
    {
        Dz = limAttDiffraction(Dz, C3);
    }
 
    return Dz;
}
 
OSpectreOctave TYAcousticModel::limAttDiffraction(const OSpectreOctave& sNC, const OSpectreOctave& C) const
{
    OSpectreOctave s;
 
    double lim20dB = 20.0;
    double lim25dB = 25.0;
    double lim0dB = 0.0;
 
    double valeur;
 
    for (unsigned int i = 0; i < sNC.getNbValues(); i++)
    {
        valeur = sNC.getTabValReel()[i];
 
        valeur = valeur < lim0dB ? lim0dB : valeur; // L'attenuation ne peut etre inferieure a 0 dB
 
        if ((C.getTabValReel()[i] - 1) <= 1e-2) // Comportement ecran mince
        {
            valeur = valeur > lim20dB ? lim20dB : valeur;
        }
        else // Comportement ecran epais ou multiple
        {
            valeur = valeur > lim25dB ? lim25dB : valeur;
        }
 
        s.getTabValReel()[i] = valeur;
    }
 
    return s;
}
 
bool TYAcousticModel::solve(TYTrajet& trajet)
{
    // Get results for each path
#ifdef _DEBUG
    std::vector<PathResults> pathsResults;
#endif
    PathResults currentPathResults;
    // Global sound level pressure
    OSpectreOctave& SLp = trajet.getSpectreOct();
    SLp.setType(SPECTRE_TYPE_LP); // Receptor spectrum is a pressure spectrum
 
    for (unsigned int i = 0; i < trajet.getNbChemins(); i++)
    {
        currentPathResults.path_id = i;
        currentPathResults.pathType = trajet.getChemin(i).getType();
        double longueur = trajet.getChemin(i).getLongueur();
 
        // Screen and ground paths results are held by direct path
        if (currentPathResults.pathType == TYTypeChemin::CHEMIN_ECRAN ||
            currentPathResults.pathType == TYTypeChemin::CHEMIN_SOL)
        {
            continue;
        }
        // Direct Ray computations
        OSpectreOctave L = OSpectreOctave(0);
 
        // Compute attenuations
        currentPathResults.Adiv = OSpectreOctave(20.0 * log10(longueur) + 11.0);
        currentPathResults.Aatm = trajet.getChemin(i).getAttenuation(TYTypeAttenuation::ATTENUATION_ATM);
        currentPathResults.Agr_s = trajet.getChemin(i).getAttenuation(TYTypeAttenuation::ATTENUATION_GND_S);
        currentPathResults.Agr_r = trajet.getChemin(i).getAttenuation(TYTypeAttenuation::ATTENUATION_GND_R);
        currentPathResults.Agr_m = trajet.getChemin(i).getAttenuation(TYTypeAttenuation::ATTENUATION_GND_M);
        currentPathResults.Dz_top = trajet.getChemin(i).getAttenuation(TYTypeAttenuation::DZ_TOP);
        currentPathResults.Dz_left = trajet.getChemin(i).getAttenuation(TYTypeAttenuation::DZ_LEFT);
        currentPathResults.Dz_right = trajet.getChemin(i).getAttenuation(TYTypeAttenuation::DZ_RIGHT);
        currentPathResults.Abar_top =
            trajet.getChemin(i).getAttenuation(TYTypeAttenuation::ATTENUATION_BAR_TOP);
        currentPathResults.Abar_left =
            trajet.getChemin(i).getAttenuation(TYTypeAttenuation::ATTENUATION_BAR_LEFT);
        currentPathResults.Abar_right =
            trajet.getChemin(i).getAttenuation(TYTypeAttenuation::ATTENUATION_BAR_RIGHT);
        currentPathResults.Abar = computeEffectiveBarAttenuation(
            currentPathResults.Abar_top, currentPathResults.Abar_left, currentPathResults.Abar_right);
 
        currentPathResults.A = currentPathResults.Adiv + currentPathResults.Aatm + currentPathResults.Agr_s +
                               currentPathResults.Agr_r + currentPathResults.Agr_m + currentPathResults.Abar;
 
        // Get source power level and source directivity LW + DC
        currentPathResults.LW = OSpectreOctave(trajet.asrc.spectrum).round(); // Round spectrum to 2 digits
                                                                              // for compliance with ISO
        // TR 17534-3 values
 
        // Get source directivity correction
        currentPathResults.Dc =
            trajet.getChemin(i).getAttenuation(TYTypeAttenuation::DIRECTIVITY_INDEX).log() * 10.0;
 
        // If path is reflected one, then compute LW_image
        if (currentPathResults.pathType == TYTypeChemin::CHEMIN_REFLEX)
        {
            currentPathResults.LW =
                currentPathResults.LW +
                trajet.getChemin(i).getAttenuation(TYTypeAttenuation::ATTENUATION_REFLEX).log() * 10.0;
        }
        currentPathResults.L = currentPathResults.LW + currentPathResults.Dc - currentPathResults.A;
        currentPathResults.L.setType(SPECTRE_TYPE_LP); // L is a pressure spectrum
        SLp = SLp.sumdB(currentPathResults.L);
 
#ifdef _DEBUG
        pathsResults.push_back(currentPathResults);
#endif
    }
 
    // Trace results
    // TODO Remove trace or keep it only in Debug
#ifdef _DEBUG
    TYSolverHelper::exportResults17534(
        pathsResults, SLp,
        *pSolverAtmos); // Export results in a csv file for comparison with 17534-3 standard
#endif
 
    trajet.build_tab_rays();
    trajet.reset(); // Erase paths array to (try to) spare memory
    return true;
}
 
OSpectreOctave TYAcousticModel::getReflexionSpectrumAt(const OSegment3D& incident, double length,
                                                       const OSegment3D& segPente,
                                                       const tympan::AcousticSource& source) const
{
    OSpectreOctave spectre;
 
    // Search for material at reflexion point
    // Set position of ray at the same high as the source
    vec3 start = OPoint3Dtovec3(incident._ptB);
    start.z = (decimal)incident._ptA._z;
    Ray ray1(start, vec3(0., 0., -1.));
    ray1.setMaxt(20000);
 
    std::list<Intersection> LI;
 
    static_cast<double>(_solver.getScene()->getAccelerator()->traverse(&ray1, LI));
 
    if (LI.empty())
    {
        start.z = (decimal)(incident._ptB._z + 1000);
        Ray ray1(start, vec3(0., 0., -1.));
        ray1.setMaxt(20000);
        static_cast<double>(_solver.getScene()->getAccelerator()->traverse(&ray1, LI));
    }
 
    assert(!LI.empty());
    unsigned int indexFace = LI.begin()->p->getPrimitiveId();
    tympan::AcousticMaterialBase* mat = _solver.getTabPolygon()[indexFace].material;
 
    // Avoid cases where the reflexion point is below a "floating" volumic source
    while (_solver.getTabPolygon()[indexFace].is_infra() &&
           source.volume_id == _solver.getTabPolygon()[indexFace].volume_id)
    {
        start.z = (decimal)min(min(_solver.getTabPolygon()[indexFace].tabPoint[0]._z,
                                   _solver.getTabPolygon()[indexFace].tabPoint[1]._z),
                               _solver.getTabPolygon()[indexFace].tabPoint[2]._z);
        Ray ray(start, vec3(0, 0, -1));
        ray.setMaxt(20000);
        std::list<Intersection> LI2;
        static_cast<double>(_solver.getScene()->getAccelerator()->traverse(&ray, LI2));
        // assert( !LI2.empty() );
        if (LI2.empty())
            break;
        indexFace = LI2.begin()->p->getPrimitiveId();
        mat = _solver.getTabPolygon()[indexFace].material;
    }
 
    // Angle estimation
    OVector3D direction(incident._ptA, incident._ptB);
    direction.normalize();
 
    // This is kept commented for the time being, even though it's the best solution
    // double angle = ( direction * -1 ).angle( _solver.getTabPolygon()[indexFace].normal );
    // angle = ABS(M_PI/2. - angle);
    double angle = (direction * -1).angle(OVector3D(incident._ptB, segPente._ptA));
    // Compute reflexion spectrum
    spectre = mat->get_absorption(angle, length);
 
    return spectre;
}
 
void TYAcousticModel::meanSlope(const OSegment3D& director, OSegment3D& slope) const
{
    // Search for primitives under the two segment extremities
 
    // To begin : initialize slope
    slope = director;
 
    // first one
    OPoint3D pt = director._ptA;
    pt._z += 1000.;
    vec3 start = OPoint3Dtovec3(pt);
    Ray ray1(start, vec3(0., 0., -1.));
 
    std::list<Intersection> LI;
 
    double distance1 = static_cast<double>(_solver.getScene()->getAccelerator()->traverse(&ray1, LI));
    assert(distance1 > 0.);
    assert(!LI.empty());
 
    unsigned int indexFace = LI.begin()->p->getPrimitiveId();
 
    // Avoid cases where the extremities are above infrastructure elements
    while (_solver.getTabPolygon()[indexFace].is_infra())
    {
        start.z = (decimal)min(min(_solver.getTabPolygon()[indexFace].tabPoint[0]._z,
                                   _solver.getTabPolygon()[indexFace].tabPoint[1]._z),
                               _solver.getTabPolygon()[indexFace].tabPoint[2]._z);
        Ray ray(start, vec3(0, 0, -1));
        ray.setMaxt(20000);
        std::list<Intersection> LI2;
        double distance = static_cast<double>(_solver.getScene()->getAccelerator()->traverse(&ray, LI2));
        // assert(distance > 0.);
        if (LI2.empty() || distance < 0.)
            break;
        distance1 += distance;
        indexFace = LI2.begin()->p->getPrimitiveId();
    }
 
    // Second one
    LI.clear();
    pt = director._ptB;
    pt._z += 1000.;
    start = OPoint3Dtovec3(pt);
    Ray ray2(start, vec3(0., 0., -1.));
 
    double distance2 = static_cast<double>(_solver.getScene()->getAccelerator()->traverse(&ray2, LI));
    // An error can occur if some elements are outside of the grip (emprise)
    assert(distance2 > 0.);
    assert(!LI.empty());
 
    indexFace = LI.begin()->p->getPrimitiveId();
 
    // Avoid cases where the extremities are above infrastructure elements
    while (_solver.getTabPolygon()[indexFace].is_infra())
    {
        start.z = (decimal)min(min(_solver.getTabPolygon()[indexFace].tabPoint[0]._z,
                                   _solver.getTabPolygon()[indexFace].tabPoint[1]._z),
                               _solver.getTabPolygon()[indexFace].tabPoint[2]._z);
        Ray ray(start, vec3(0, 0, -1));
        ray.setMaxt(20000);
        std::list<Intersection> LI2;
        double distance = static_cast<double>(_solver.getScene()->getAccelerator()->traverse(&ray, LI2));
        // assert(distance > 0.);
        if (LI2.empty() || distance < 0.)
            break;
        distance2 += distance;
        indexFace = LI2.begin()->p->getPrimitiveId();
    }
 
    // Compute projection on the ground of segment points suppose sol is under the points ...
    slope._ptA._z = director._ptA._z - (distance1 - 1000.);
    slope._ptB._z = director._ptB._z - (distance2 - 1000.);
}
 
bool TYAcousticModel::computeSegmentEdgesHeights(double& hauteurA, double& hauteurB,
                                                 const OSegment3D& meanSlope, const OSegment3D& ray) const
{
    bool res = true;
    hauteurA = ray._ptA._z - meanSlope._ptA._z;
    hauteurB = ray._ptB._z - meanSlope._ptB._z;
    return res;
}