TYChemin9613Solver2024.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 "Tympan/solvers/9613/9613Solver2024/TYChemin9613Solver2024.h"
 
OSpectreOctave TYChemin9613Solver2024::calcTotalGroundAttenuation(const OSpectreOctave& Agr_s,
                                                                  const OSpectreOctave& Agr_r,
                                                                  const OSpectreOctave& Agr_m, double hs,
                                                                  double hr, double dp) const
{
    OSpectreOctave Agr9613 = TYChemin9613Solver::calcTotalGroundAttenuation(Agr_s, Agr_r, Agr_m, hs, hr, dp);
    double Kgeo = (dp * dp + (hs - hr) * (hs - hr)) / (dp * dp + (hs + hr) * (hs + hr));
    OSpectreOctave s = Agr9613 * (-1.0);
    s = s.div(10.0);
    s = s.exp(std::log(10));
    s = (s + (-1.0)) * Kgeo + 1.0;
    s = s.log() * (-10.0);
 
    return s;
}
 
OSpectreOctave TYChemin9613Solver2024::calcMinimalExtensionConditionOneReflection(const OPoint3D& Oprev,
                                                                                  const OPoint3D& O,
                                                                                  double d_SO, double d_OR,
                                                                                  double a, double h,
                                                                                  const OVector3D& n) const
{
    // first we need to set a convention to the normal of the surface:
    // (1) it is a normal of a vertical plane (vertical component = 0)
    // (2) it is normalized
    // (3) from O, it points towards S and R
    OVector3D nCopy{n}; // a copy is necessary because coordinates cannot be retrieved for a const OVector3D
    const double nx = nCopy.getCoords()[0];
    const double ny = nCopy.getCoords()[1];
    OVector3D nConventional{nx, ny, 0.}; // condition (1)
    nConventional.normalize();           // condition (2)
    const OVector3D OSvec{O, Oprev};
    if (OSvec.scalar(nConventional) < 0.) // condition (3)
    {
        nConventional.invert();
    }
 
    // to compute cos(beta_a) and cos(beta_h), we need to perform 2D geometry in two planes:
    // - a horizontal plane passing by O --> cos(beta_a)
    // - a vertical plane passing by O --> cos(beta_h)
    // for that, we define normalized tangential vectors of the reflecting surface
    // - tHorizontal
    // - tVertical
    const OVector3D tVertical{0., 0., 1.};
    OVector3D tHorizontal = tVertical.cross(nConventional);
    tHorizontal.normalize();
 
    // cos(beta_a)
    const double OSvecHorizontalT = OSvec.scalar(tHorizontal);
    const double OSvecN = OSvec.scalar(nConventional);
    const OVector3D OSvecHorizontal = tHorizontal * OSvecHorizontalT + nConventional * OSvecN;
    const double dSOHorizontal = OSvecHorizontal.norme();
    const double cos_beta_a = OSvecN / dSOHorizontal;
 
    // cos(beta_h)
    const double OSvecVerticalT = OSvec.scalar(tVertical);
    const OVector3D OSvecVertical = tVertical * OSvecVerticalT + nConventional * OSvecN;
    const double dSOVertical = OSvecVertical.norme();
    const double cos_beta_h = OSvecN / dSOVertical;
 
    // l_eff (Eq. 27)
    const double l_eff = std::min(a * cos_beta_a, h * cos_beta_h);
 
    // f_threshold (Eq. 26 multiplied by c)
    const double c_9613_2 = 340.0; // TODO: to factorize by using somehow the constant defined in
                                   // TYAcousticModel::computeWaveLength
    const double fThreshold = c_9613_2 * (2. / (l_eff * l_eff)) * (d_SO * d_OR) / (d_SO + d_OR);
 
    // filter spectrum
    OSpectreOctave filter{0.};
    OTabFreq frequencies = OSpectreOctave::getTabFreqExact();
    for (unsigned int i = 0; i < frequencies.size(); i++)
    {
        if (frequencies[i] > fThreshold)
        {
            filter.getTabValReel()[i] = 1.;
        }
    }
 
    return filter;
}
 
void TYChemin9613Solver2024::calcCylinderReflectionAttenuation(const OPoint3D& S, const OPoint3D& R,
                                                               const OPoint3D& P, const OPoint3D& M, double r,
                                                               const OVector3D& axis)
{
    // in 9613-2:2024 standard, geometric computation have to be performed "on a plane perpendicular to the
    // cylinder-axis"
 
    assert(r > 0.);
 
    // build an orthonormal basis of a plane orthogonal to axis
    OVector3D u1, u2;
    if (abs(axis._x) < EPSILON_13 && abs(axis._y) < EPSILON_13)
    {
        // axis is colinear to (0, 0, 1)
        // in this case, we get directly an orthonormal basis using e1 and e2
        u1 = OVector3D{1., 0., 0.};
        u2 = OVector3D{0., 1., 0.};
    }
    else
    {
        // axis is not colinear to (0, 0, 1)
        // in this case, we have to build an orthonormal basis:
        // - define a candidate u1_prime, non-colinear with axis
        // - build u1, orthogonal to axis using the Gram-Schmidt process
        // - build u2 using cross product
        // - normalize u1 and u2
        const OVector3D u1_prime{0., 0., 1.};
        const OVector3D u1_non_normalized{u1_prime - axis * (u1_prime.scalar(axis) / axis.scalar(axis))};
        const OVector3D u2_non_normalized{u1_non_normalized.cross(axis)};
        u1 = u1_non_normalized * (1. / u1_non_normalized.norme());
        u2 = u2_non_normalized * (1. / u2_non_normalized.norme());
    }
 
    // project S, R, P and M in this basis with origin (0, 0, 0)
    const OPoint3D S2D = projectOnPlane(S, u1, u2);
    const OPoint3D R2D = projectOnPlane(R, u1, u2);
    const OPoint3D P2D = projectOnPlane(P, u1, u2);
    const OPoint3D M2D = projectOnPlane(M, u1, u2);
 
    // compute intermediate variables dS, dR, d, and k
    const double dS{S2D.dist2DFrom(P2D)};
    const double dR{R2D.dist2DFrom(P2D)};
    // d is the norm of a vector v fulfilling these properties:
    //    v . SP = 0
    //    MP + v = mu SP
    // thus:
    //    v = mu SP - MP, with mu = MP . SP / SP^2
    const OVector3D MP{M2D, P2D};
    const OVector3D SP{S2D, P2D};
    const OVector3D v{SP * (MP.scalar(SP) / (dS * dS)) - MP};
    const double d{v.norme()};
    const double k{d / r};
 
    assert(k <= 1.);
 
    // compute Acurv using Eq. 30 of 9613-2:2024 standard
    const double Acurv = 10. * log10(1. + 2 * dS * dR / r / (dS + dR) / sqrt(1. - k * k));
 
    OSpectreOctave AcurvSpectrum{Acurv};
    setAttenuation(TYTypeAttenuation::ATTENUATION_CYLINDER, AcurvSpectrum);
}
 
OPoint3D TYChemin9613Solver2024::projectOnPlane(const OPoint3D& P, const OVector3D& u1, const OVector3D& u2)
{
    return OPoint3D{u1.scalar(P), u2.scalar(P), 0.};
}