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& S, const OPoint3D& O, const OPoint3D& R, 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, S};
    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);
 
    // distances
    const double d_SO = OSvec.norme();
    const OVector3D ORvec{O, R};
    const double d_OR = ORvec.norme();
 
    // 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;
}