TYAcousticModel9613Solver.cpp

Go to the documentation of this file.

/*
 * Copyright (C) <2012-2024> <EDF-DTG> <FRANCE>
 * This file is part of Code_TYMPAN (R).
 * Code_TYMPAN (R) 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 3 of the License, or
 * (at your option) any later version.
 * Code_TYMPAN (R) 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 Code_TYMPAN (R). If not, see <https://www.gnu.org/licenses/>.
 */
 
#include <deque>
#include <list>
#include <cmath>
#include <algorithm>
#include "Tympan/core/defines.h"
#include "Tympan/models/solver/config.h"
#include "Tympan/models/common/mathlib.h"
#include "Tympan/geometric_methods/AcousticRaytracer/Geometry/Scene.h"
#include "TYAcousticModel9613Solver.h"
#include "TYSolver9613Solver.h"
#include "TYSolverHelper.h"
 
#include "Tympan/solvers/9613/9613Solver/TYChemin9613Solver.h"
 
TYAcousticModel9613Solver::TYAcousticModel9613Solver(TYSolver9613Solver& solver)
    : TYAcousticModel(), _solver(solver)
{
    _absoNulle = OSpectreOctave(1.0);
    _absoNulle.setType(SPECTRE_TYPE_ABSO); // Spectre d'absorption
}
 
void TYAcousticModel9613Solver::init()
{
    TYAcousticModel::init();
}
 
void TYAcousticModel9613Solver::computeWaveLength()
{
    // Compute wave length
    double c_9613_2 = 340.0;
    _lambda = OSpectreOctave::getLambda(c_9613_2);
}
 
void TYAcousticModel9613Solver::compute(const std::deque<TYSIntersection>& tabIntersect,
                                        TYTrajet9613Solver& 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();
 
    TYTabChemin9613Solver& 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();
}
 
bool TYAcousticModel9613Solver::computeCheminsAvecEcran(const OSegment3D& ray,
                                                        const tympan::AcousticSource& source,
                                                        const TabPoint3D& pts, const bool vertical,
                                                        TYTabChemin9613Solver& 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]);
 
    TYTabEtape9613Solver tabSteps;
    double pathLength = 0.0;
 
    /*--- BEFORE OBSTACLE ---*/
 
    TYTabEtape9613Solver steps;
    OSegment3D curSegment(ray._ptA, firstPt);
    double tempLength = curSegment.longueur();
 
    bool bPathOk = addGroundSteps(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;
    TYEtape9613Solver 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();
 
    addGroundSteps(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._type = TYRECEPTEUR;
    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 ---*/
 
    TYChemin9613Solver 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 TYAcousticModel9613Solver::addGroundSteps(const OPoint3D& ptStart, const OPoint3D& ptEnd,
                                               const tympan::AcousticSource& source, const bool& fromSource,
                                               TYTabEtape9613Solver& steps) const
{
    bool res = true;
 
    TYEtape9613Solver 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 TYAcousticModel9613Solver::computeCheminSansEcran(const std::deque<TYSIntersection>& tabIntersect,
                                                       const OSegment3D& ray,
                                                       const tympan::AcousticSource& source,
                                                       TYTabChemin9613Solver& TabChemin, double distance,
                                                       bool conditionFav) const
{
    /*
        COMPUTATION FOR A ROUTE WITHOUT OBSTACLE CONSISTS OF ONE DIRECT PATH
    */
    // Add direct path
    std::unique_ptr<TYChemin9613Solver> chemin = createChemin();
 
    // 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();
 
    TYTabEtape9613Solver Etapes;
    double Gs{0.0}, Gm{0.0}, Gr{0.0};
    double hs{0.0}, hr{0.0};
 
    computeSegmentEdgesHeights(hs, hr, segMeanSlope, ray);
 
    addGroundSteps(ray._ptA, ray._ptB, source, true, Etapes);
    getGroundfactors(tabIntersect, SR2D, hs, hr, Gs, Gm, Gr);
 
    // Add direct path
    chemin->setType(TYTypeChemin::CHEMIN_DIRECT);
    TYTabEtape9613Solver tabEtapes;
    tabEtapes.push_back(Etapes[0]); // Add direct step
    chemin->setLongueur(distance);  // In this case, length = 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 TYAcousticModel9613Solver::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 TYAcousticModel9613Solver::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 TYAcousticModel9613Solver::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;
}
 
bool TYAcousticModel9613Solver::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);
    std::unordered_map<OVector3D, double, OVector3DHash> mapResultSegFactorG;
 
    // 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;
    dpZone = 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};
            auto it = mapResultSegFactorG.find(result);
            if (it != mapResultSegFactorG.end())
            {
                GZone = GZone + 0.5 * result.norme() * (currentG - it->second);
            }
            else
            {
                GZone = GZone + result.norme() * currentG;
                dpZone += result.norme();
            }
            mapResultSegFactorG[result] = currentG;
        }
    }
 
    return ret;
}
 
void TYAcousticModel9613Solver::computeCheminReflexion(const std::deque<TYSIntersection>& tabIntersect,
                                                       const OSegment3D& ray,
                                                       const tympan::AcousticSource& source,
                                                       TYTabChemin9613Solver& 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();
 
    OSpectreOctave filter{1.};
 
    // 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;
 
                TYTabEtape9613Solver tabEtapes;
 
                double pathLength = upwardSeg.longueur() + downwardSeg.longueur();
 
                TYEtape9613Solver 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.selectFaces(tabIntersectUpSegment, upwardSeg, source.volume_id);
                _solver.selectFaces(tabIntersectDownSegment, downwardSeg, source.volume_id);
 
                // Compute ground factors for reflected ray
                getGroundfactors(tabIntersectUpSegment, tabIntersectDownSegment, raySO, rayOR, hs, hr, Gs, Gm,
                                 Gr);
 
                std::unique_ptr<TYChemin9613Solver> chemin = createChemin();
                chemin->setType(TYTypeChemin::CHEMIN_REFLEX);
                chemin->setLongueur(pathLength);
                chemin->setDistance(distance);
                chemin->calcAttenuation(tabEtapes, *_pSolverAtmos, dp, hs, hr, Gs, Gm, Gr);
                // Set minimal extension condition reflection
                chemin->setMinimalExtensionConditionReflection(filter);
 
                TabChemins.push_back(*chemin); // Put the reflected path in paths table
                tabEtapes.clear();
            }
        }
    }
}
 
OSpectreOctave TYAcousticModel9613Solver::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 TYAcousticModel9613Solver::calculAttDiffraction(const OSegment3D& ray, const double& re,
                                                               const double& dss, const double& dsr,
                                                               const double& width,
                                                               const bool& vertical) const
{
    double rd;
 
    OSpectreOctave C3 = calculC3(width); // Corrective factor linked with screen width
 
    double C2{20.0};
    OSpectreOctave zmin = calculZMin(C2, C3);
 
    rd = ray.longueur();
 
    double z = re - rd; // Path-length difference
    z = z <= 0 ? 0.0 : z;
 
    OSpectreOctave Kmeteo = calculKmeteo(vertical, dss, dsr, rd, z, width, zmin);
 
    OSpectreOctave Dz = calculDz(C2, z, C3, Kmeteo, zmin);
    // 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 TYAcousticModel9613Solver::calculZMin(const double C2, const OSpectreOctave& C3) const
{
    // zmin = 0.0
    OSpectreOctave zmin{0.0};
    return zmin;
}
 
OSpectreOctave TYAcousticModel9613Solver::calculKmeteo(const bool vertical, const double d_SS,
                                                       const double d_SR, const double d, const double z,
                                                       const double e, const OSpectreOctave& z_min) const
{
    // K_meteo = exp[-(1 / 2000) * sqrt(d_SS * d_SR * d/(2z))] for z > 0 and vertical diffraction
    // K_meteo = 1 for z <= 0 or for lateral diffraction
    double Kmeteo = 1.0;
    if (z > 0.0 && vertical)
    {
        Kmeteo = exp(-(1.0 / 2000.0) * sqrt(d_SS * d_SR * d / (2 * z)));
    }
    return OSpectreOctave{Kmeteo};
}
 
OSpectreOctave TYAcousticModel9613Solver::calculDz(const double z, const double C2, const OSpectreOctave& C3,
                                                   const OSpectreOctave& Kmeteo,
                                                   const OSpectreOctave& zmin) const
{
    // Dz = 10 * log10[3 + (C2 / lambda) * C3 * z * Kmeteo] dB for z > zmin
    // Dz = 0 dB for z <= zmin
 
    OSpectreOctave Dz{0.0};
 
    double* dz = Dz.getTabValReel();
    const double* lambda = _lambda.getTabValReel();
    const double* c3 = C3.getTabValReel();
    const double* kmeteo = Kmeteo.getTabValReel();
    const double* zMin = zmin.getTabValReel();
 
    for (std::size_t i = 0; i < TY_SPECTRE_OCT_NB_ELMT; ++i)
    {
        if (z > zMin[i])
        {
            dz[i] = 10.0 * std::log10(3.0 + (C2 / lambda[i]) * c3[i] * z * kmeteo[i]);
        }
    }
 
    return Dz;
}
 
OSpectreOctave TYAcousticModel9613Solver::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 TYAcousticModel9613Solver::solve(TYTrajet9613Solver& 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
    tympan::LPSolverConfiguration config = tympan::SolverConfiguration::get();
 
    for (unsigned int i = 0; i < trajet.getNbChemins(); i++)
    {
        currentPathResults.path_id = i;
        currentPathResults.pathType = trajet.getChemin(i).getType();
        float minSRDistance = config->MinSRDistance;
        double longueur = (minSRDistance > trajet.getChemin(i).getLongueur())
                              ? minSRDistance
                              : trajet.getChemin(i).getLongueur();
 
        // Screen and ground paths results are held by direct path
        if (currentPathResults.pathType == TYTypeChemin::CHEMIN_ECRAN)
        {
            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.Agr = trajet.getChemin(i).getAttenuation(TYTypeAttenuation::ATTENUATION_GND_TOTAL);
        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.Acurv = OSpectreOctave{0.};
        currentPathResults.Amisc = OSpectreOctave{0.};
        if (trajet.getChemin(i).hasAttenuation(TYTypeAttenuation::ATTENUATION_CYLINDER))
        {
            currentPathResults.Acurv =
                trajet.getChemin(i).getAttenuation(TYTypeAttenuation::ATTENUATION_CYLINDER);
            currentPathResults.Amisc = currentPathResults.Amisc + currentPathResults.Acurv;
        }
 
        currentPathResults.A = currentPathResults.Adiv + currentPathResults.Aatm + currentPathResults.Agr +
                               currentPathResults.Abar + currentPathResults.Amisc;
 
        // 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;
        if (currentPathResults.pathType == TYTypeChemin::CHEMIN_REFLEX)
        {
            OSpectreOctave filter =
                trajet.getChemin(i).getAttenuation(TYTypeAttenuation::REFLEXION_MIN_EXT_COND);
            OSpectreOctave filterComplement{filter};
            filterComplement = filterComplement * (-1.0) + 1.0;
            currentPathResults.L = currentPathResults.L * filter + OSpectreOctave{} * filterComplement;
        }
        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;
}
 
bool TYAcousticModel9613Solver::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;
}