TYChemin9613Solver.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/9613Solver/TYChemin9613Solver.h"
#include "Tympan/models/solver/config.h"
 
TYChemin9613Solver::TYChemin9613Solver() : TYChemin()
{
    _attenuations.clear();
}
 
TYChemin9613Solver::TYChemin9613Solver(const TYChemin9613Solver& other)
{
    *this = other;
}
 
TYChemin9613Solver::~TYChemin9613Solver()
{
    _attenuations.clear();
}
 
TYChemin9613Solver& TYChemin9613Solver::operator=(const TYChemin9613Solver& other)
{
    if (this != &other)
    {
        TYChemin::operator=(other);
        _attenuations = other._attenuations;
    }
 
    return *this;
}
 
bool TYChemin9613Solver::operator==(const TYChemin9613Solver& other) const
{
    if (this != &other)
    {
        if (TYChemin::operator!=(other))
        {
            return false;
        }
        if (_attenuations != other._attenuations)
        {
            return false;
        }
    }
 
    return true;
}
 
bool TYChemin9613Solver::operator!=(const TYChemin9613Solver& other) const
{
    return !operator==(other);
}
 
void TYChemin9613Solver::calcAttenuation(const TYTabEtape9613Solver& tabEtapes,
                                         const AtmosphericConditions& atmos, double dp, double hs, double hr,
                                         double Gs, double Gm, double Gr)
{
    unsigned int i = 0;
 
    switch (_typeChemin)
    {
        case TYTypeChemin::CHEMIN_DIRECT:
            _attenuations[TYTypeAttenuation::DIRECTIVITY_INDEX] =
                tabEtapes[0]._Absorption; // S = Source directivity index
            _attenuations[TYTypeAttenuation::ATTENUATION_ATM] =
                atmos.compute_length_absorption_oct(_longueur); // Aatm = alpha * d
            calcGroundAttenuations(dp, hs, hr, Gs, Gm, Gr);
            break;
 
        case TYTypeChemin::CHEMIN_REFLEX:
            _attenuations[TYTypeAttenuation::ATTENUATION_ATM] =
                atmos.compute_length_absorption_oct(_longueur); // Aatm = alpha * d
            calcGroundAttenuations(dp, hs, hr, Gs, Gm, Gr);
            _attenuations[TYTypeAttenuation::DIRECTIVITY_INDEX] =
                tabEtapes[0]._Absorption; // S = Source directivity index
            // We multiply steps absorptions from second
            _attenuations[TYTypeAttenuation::ATTENUATION_REFLEX] = OSpectreOctave(1.0);
            for (i = 1; i < tabEtapes.size(); i++)
            {
                _attenuations[TYTypeAttenuation::ATTENUATION_REFLEX] =
                    _attenuations[TYTypeAttenuation::ATTENUATION_REFLEX] *
                    (OSpectreOctave(1.0) -
                     tabEtapes[i]._Absorption); // Product of directivty factor and reflexion coefficients
            }
            break;
 
        default:
            break;
    }
 
    build_eq_path(tabEtapes);
}
 
void TYChemin9613Solver::calcMinimalExtensionCondition(const OPoint3D& S, const OPoint3D& O,
                                                       const OPoint3D& R, double a, double h,
                                                       const OVector3D& n)
{
    // this code assumes 1st order reflection
    _attenuations[TYTypeAttenuation::REFLEXION_MIN_EXT_COND] =
        calcMinimalExtensionConditionOneReflection(S, O, R, a, h, n);
}
 
void TYChemin9613Solver::computeBarAttenuation(const OSpectreOctave& Dz, const bool vertical, const bool left)
{
    if (vertical)
    {
        _attenuations[TYTypeAttenuation::DZ_TOP] = Dz;
        _attenuations[TYTypeAttenuation::ATTENUATION_BAR_TOP] = Dz;
        OSpectreOctave Agr = _attenuations[TYTypeAttenuation::ATTENUATION_GND_TOTAL];
        for (unsigned int i = 0; i < TY_SPECTRE_OCT_NB_ELMT; i++)
        {
            if (Agr.getTabValReel()[i] > 0)
            {
                _attenuations[TYTypeAttenuation::ATTENUATION_BAR_TOP].getTabValReel()[i] -=
                    Agr.getTabValReel()[i];
            }
            if (_attenuations[TYTypeAttenuation::ATTENUATION_BAR_TOP].getTabValReel()[i] < 0)
            {
                _attenuations[TYTypeAttenuation::ATTENUATION_BAR_TOP].getTabValReel()[i] = 0;
            }
        }
    }
    else
    {
        if (left)
        {
            _attenuations[TYTypeAttenuation::DZ_LEFT] = Dz;
            _attenuations[TYTypeAttenuation::ATTENUATION_BAR_LEFT] = Dz;
        }
        else
        {
            _attenuations[TYTypeAttenuation::DZ_RIGHT] = Dz;
            _attenuations[TYTypeAttenuation::ATTENUATION_BAR_RIGHT] = Dz;
        }
    }
}
 
OSpectreOctave& TYChemin9613Solver::getAttenuation(TYTypeAttenuation type)
{
    return _attenuations[type];
}
 
void TYChemin9613Solver::setAttenuationBarWhenNoPath(bool vertical, bool left)
{
    if (vertical)
    {
        _attenuations[TYTypeAttenuation::ATTENUATION_BAR_TOP] = OSpectreOctave{100.0};
    }
    else
    {
        if (left)
        {
            _attenuations[TYTypeAttenuation::ATTENUATION_BAR_LEFT] = OSpectreOctave{100.0};
        }
        else
        {
            _attenuations[TYTypeAttenuation::ATTENUATION_BAR_RIGHT] = OSpectreOctave{100.0};
        }
    }
}
 
void TYChemin9613Solver::calcGroundAttenuations(double dp, double hs, double hr, double Gs, double Gm,
                                                double Gr)
{
    // Main method for flat ground or with constant slope
    // Compute zones dimensions
    bool bHasIntermediateZone = dp > 30 * (hs + hr);
 
    double q = 0.0;
    if (bHasIntermediateZone)
    {
        q = 1 - 30 * (hs + hr) / dp;
    }
 
    // Compute ground attenuations for source, receptor and middle zones
    _attenuations[TYTypeAttenuation::ATTENUATION_GND_S] = calcGroundAttenuationSR(dp, hs, Gs);
    _attenuations[TYTypeAttenuation::ATTENUATION_GND_R] = calcGroundAttenuationSR(dp, hr, Gr);
    _attenuations[TYTypeAttenuation::ATTENUATION_GND_M] = calcGroundAttenuationM(q, Gm);
    _attenuations[TYTypeAttenuation::ATTENUATION_GND_TOTAL] =
        calcTotalGroundAttenuation(_attenuations[TYTypeAttenuation::ATTENUATION_GND_S],
                                   _attenuations[TYTypeAttenuation::ATTENUATION_GND_R],
                                   _attenuations[TYTypeAttenuation::ATTENUATION_GND_M], hs, hr, dp);
}
 
OSpectreOctave TYChemin9613Solver::calcGroundAttenuationSR(double dp, double h, double G)
{
    double a = 1.5 + 3.0 * exp(-0.12 * (h - 5) * (h - 5)) * (1 - exp(-dp / 50.0)) +
               5.7 * exp(-0.09 * h * h) * (1 - exp(-2.8 * pow(10, -6) * dp * dp));
    double b = 1.5 + 8.6 * exp(-0.09 * h * h) * (1 - exp(-dp / 50.0));
    double c = 1.5 + 14.0 * exp(-0.46 * h * h) * (1 - exp(-dp / 50.0));
    double d = 1.5 + 5.0 * exp(-0.9 * h * h) * (1 - exp(-dp / 50.0));
 
    const double As[9]{-1.5,         -1.5,           -1.5 + G * a,   -1.5 + G * b,  -1.5 + G * c,
                       -1.5 + G * d, -1.5 * (1 - G), -1.5 * (1 - G), -1.5 * (1 - G)};
    return OSpectreOctave{As, 9, 0};
}
 
OSpectreOctave TYChemin9613Solver::calcGroundAttenuationM(double q, double Gm)
{
    OSpectreOctave Agr_m = OSpectreOctave{0.0};
    if (q != 0.0)
    {
        const double Am[9]{-3.0 * q,          -3.0 * q,          -3 * q * (1 - Gm),
                           -3 * q * (1 - Gm), -3 * q * (1 - Gm), -3 * q * (1 - Gm),
                           -3 * q * (1 - Gm), -3 * q * (1 - Gm), -3 * q * (1 - Gm)};
        Agr_m = OSpectreOctave{Am, 9, 0};
    }
    return Agr_m;
}
 
OSpectreOctave TYChemin9613Solver::calcTotalGroundAttenuation(const OSpectreOctave& Agr_s,
                                                              const OSpectreOctave& Agr_r,
                                                              const OSpectreOctave& Agr_m, double hs,
                                                              double hr, double dp) const
{
    return Agr_s + Agr_r + Agr_m;
}
 
OSpectreOctave TYChemin9613Solver::calcMinimalExtensionConditionOneReflection(
    const OPoint3D& S, const OPoint3D& O, const OPoint3D& R, double a, double h, const OVector3D& n) const
{
    OSpectreOctave allWavelengthsAreReflected{1.};
    return allWavelengthsAreReflected;
}