TYTrajetDefaultSolver.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/DefaultSolver/TYTrajetDefaultSolver.h"
#include "Tympan/models/solver/config.h"
 
TYTrajetDefaultSolver::TYTrajetDefaultSolver(tympan::AcousticSource& asrc_, tympan::AcousticReceptor& arcpt_)
    : TYTrajet(asrc_, arcpt_)
{
}
 
TYTrajetDefaultSolver::~TYTrajetDefaultSolver()
{
    reset();
}
 
void TYTrajetDefaultSolver::reset()
{
    _chemins.clear();
    _cheminsDirect.clear();
}
 
TYTrajetDefaultSolver& TYTrajetDefaultSolver::operator=(const TYTrajetDefaultSolver& other)
{
    if (this != &other)
    {
        TYTrajet::operator=(other);
        _chemins = other._chemins;
        _sLP = other._sLP;
    }
    return *this;
}
 
bool TYTrajetDefaultSolver::operator==(const TYTrajetDefaultSolver& other) const
{
    if (this != &other)
    {
        if (TYTrajet::operator!=(other))
        {
            return false;
        }
        if (_chemins != other._chemins)
        {
            return false;
        }
        if (_sLP != other._sLP)
        {
            return false;
        };
        // if (asrc != other.asrc) { return false; };
        // if (arcpt != other.arcpt) ;
    }
 
    return true;
}
 
bool TYTrajetDefaultSolver::operator!=(const TYTrajetDefaultSolver& other) const
{
    return !operator==(other);
}
 
void TYTrajetDefaultSolver::addChemin(const TYCheminDefaultSolver& chemin)
{
    _chemins.push_back(chemin);
}
 
void TYTrajetDefaultSolver::addCheminDirect(const TYCheminDefaultSolver& chemin)
{
    _cheminsDirect.push_back(chemin);
}
 
OSpectre TYTrajetDefaultSolver::getPEnergetique(const AtmosphericConditions& atmos)
{
    OSpectre s = OSpectre::getEmptyLinSpectre();
    OSpectreComplex sTemp;
    int firstReflex = -1;
    unsigned int indiceDebutEffetEcran = 0;
    unsigned int i = 0;
 
    // On calcule l'attenuation sur le trajet direct (sauf chemins reflechis).
    for (i = 0; i < this->_chemins.size(); i++)
    {
        // Si un ecran est present, on ne traite pas les reflexions (dans un premier temp ...)
        if ((_chemins[0].getType() == TYTypeChemin::CHEMIN_ECRAN) &&
            (_chemins[i].getType() == TYTypeChemin::CHEMIN_REFLEX))
        {
            firstReflex = i;
            break;
        }
        sTemp = _chemins[i].getAttenuation();
        s = s.sum(sTemp.mult(sTemp)); // somme des carres des modules
    }
 
    // Dans le cas d'un ecran, on compare l'attenuation obtenue a celle du trajet direct
    // pour eviter les effets d'amplification (plus de bruit avec l'ecran que sans ecran ...)
    if (_chemins[0].getType() == TYTypeChemin::CHEMIN_ECRAN)
    {
        OSpectre attDirect = OSpectre::getEmptyLinSpectre();
 
        for (i = 0; i < _cheminsDirect.size(); i++)
        {
            sTemp = _cheminsDirect[i].getAttenuation();
            attDirect = attDirect.sum(sTemp.mult(sTemp));
        }
 
        // On regarde l'attenuation globale obtenue pour chaque frequence,
        // on la compare a celle obtenue sur le trajet sans ecran,
        // si elle est superieure a 1 alors on prend la valeur obtenue pour le trajet sans ecran
        for (i = 0; i < s.getNbValues(); i++)
        {
            if (s.getTabValReel()[i] < attDirect.getTabValReel()[i])
            {
                indiceDebutEffetEcran = i; // On prend note de l'indice
                break; // Si l'ecran commence a attenuer plus que le trajet direct, il faut sortir de la
                       // boucle
            }
        } //*/
 
        if (firstReflex != -1) // S'il y a une reflexion sur un ecran
        {
            // On rajoute la contribution des chemins reflechis
            // 1. Aux chemins normaux et aux chemins directs
            for (i = firstReflex; i < _chemins.size(); i++)
            {
                sTemp = _chemins[i].getAttenuation().mult(_chemins[i].getAttenuation());
                s = s.sum(sTemp);
                attDirect = attDirect.sum(sTemp);
            }
 
            // On remplace la contribution du trajet direct pour toutes les frequences ou cela est necessaire
            for (i = 0; i < indiceDebutEffetEcran; i++)
            {
                s.getTabValReel()[i] = attDirect.getTabValReel()[i];
            }
        }
    }
    build_tab_rays();
    _chemins.clear(); // On efface le tableau des chemins pour (essayer de) gagner de la place en memoire
    _cheminsDirect.clear();
    return s;
}
 
OSpectre TYTrajetDefaultSolver::getPInterference(const AtmosphericConditions& atmos)
{
    unsigned int i = 0, j = 0;
    int firstReflex = -1;
    unsigned int indiceDebutEffetEcran = 0;
    bool ecranFound = false;
 
    if (_chemins.size())
    {
        ecranFound = (_chemins[0].getType() == TYTypeChemin::CHEMIN_ECRAN);
    }
 
    // On recupere prealablement l'ensemble des attenuations et les longueurs des chemins
    std::vector<OSpectreComplex> tabSpectreAtt;
    std::vector<OSpectreComplex> tabSpectreAttDirect;
    std::vector<double> tabLongueur;
    std::vector<double> tabLongueurDirect;
 
    for (i = 0; i < _chemins.size(); i++)
    {
        tabSpectreAtt.push_back(_chemins[i].getAttenuation());
        tabLongueur.push_back(_chemins[i].getLongueur());
    }
 
    // On calcule la somme des carres des modules et la somme des produits croises
 
    OSpectre sCarreModule = OSpectre::getEmptyLinSpectre();
    OSpectre sProduitCroise = OSpectre::getEmptyLinSpectre();
    OSpectre sTemp;
 
    for (i = 0; i < _chemins.size(); i++)
    {
 
        // Si un ecran est present, on ne traite pas les reflexions (dans un premier temp ...)
        if (ecranFound && (_chemins[i].getType() == TYTypeChemin::CHEMIN_REFLEX))
        {
            firstReflex = i;
            break;
        }
        // on fait la somme du carre des modules
        sCarreModule = sCarreModule.sum(tabSpectreAtt[i].mult(tabSpectreAtt[i]));
 
        // on calcule les produits croises avec les autres chemins
        for (j = i + 1; j < _chemins.size(); j++)
        {
            // On procedera aux produits croise avec les chemins reflechis plus loin ...
            if (ecranFound && (_chemins[j].getType() == TYTypeChemin::CHEMIN_REFLEX))
            {
                continue;
            }
 
            sTemp = tabSpectreAtt[i].mult(tabSpectreAtt[j].mult(2.0));
            sTemp = sTemp.mult(
                correctTiers(tabSpectreAtt[i], tabSpectreAtt[j], atmos, tabLongueur[i], tabLongueur[j]));
            sProduitCroise = sProduitCroise.sum(sTemp);
        }
    }
 
    OSpectre s = sCarreModule.sum(sProduitCroise); //.abs() ;
 
    // Dans le cas d'un ecran, on compare l'attenuation obtenue a celle du trajet direct
    // pour eviter les effets d'amplification (plus de bruit avec l'ecran que sans ecran ...)
 
    if (ecranFound) // On comparera au carre des modules des trajets directs
    {
        // On calcule l'attenuation sur le trajet direct
        OSpectre sCarreModuleDirect = OSpectre::getEmptyLinSpectre();
        OSpectre sProduitCroiseDirect = OSpectre::getEmptyLinSpectre();
 
        for (i = 0; i < _cheminsDirect.size(); i++)
        {
            tabSpectreAttDirect.push_back(_cheminsDirect[i].getAttenuation());
            tabLongueurDirect.push_back(_cheminsDirect[i].getLongueur());
        }
 
        for (i = 0; i < _cheminsDirect.size(); i++)
        {
            // on fait la somme du carre des modules
            sCarreModuleDirect = sCarreModuleDirect.sum(tabSpectreAttDirect[i].mult(tabSpectreAttDirect[i]));
 
            // on calcule les produits croises avec les autres chemins
            for (j = i + 1; j < _cheminsDirect.size(); j++)
            {
                sTemp = tabSpectreAttDirect[i].mult(tabSpectreAttDirect[j].mult(2.0));
                sTemp = sTemp.mult(correctTiers(tabSpectreAttDirect[i], tabSpectreAttDirect[j], atmos,
                                                tabLongueurDirect[i], tabLongueurDirect[j]));
                sProduitCroiseDirect = sProduitCroiseDirect.sum(sTemp);
            }
        }
 
        OSpectre attDirect = sCarreModuleDirect.sum(sProduitCroiseDirect); //.abs() ;
 
        // On compare l'attenuation sur le trajet "normal" en energetique a
        // l'attenuation sur le trajet direct en energetique.
        // Si elle est superieure, alors on prend l'attenuation sur le trajet direct
 
        for (j = 0; j < s.getNbValues(); j++)
        {
            if (s.getTabValReel()[j] < attDirect.getTabValReel()[j])
            {
                indiceDebutEffetEcran = j; // On prend note de l'indice
                break; // Si l'ecran commence a attenuer plus que le trajet direct, il faut sortir de la
                       // boucle
            }
        }
 
        if (firstReflex != -1) // S'il y a une reflexion sur un ecran
        {
            // On rajoute la contribution des chemins reflechis:
            // 1. aux chemins "normaux"
            for (i = firstReflex; i < _chemins.size(); i++)
            {
                sCarreModule = sCarreModule.sum(tabSpectreAtt[i].mult(tabSpectreAtt[i]));
 
                // on calcule les produits croises avec les autres chemins
                for (j = 0; j < _chemins.size(); j++)
                {
                    if (j == i)
                    {
                        continue;
                    } // pas avec lui meme
 
                    sTemp = tabSpectreAtt[i].mult(tabSpectreAtt[j].mult(2.0));
                    sTemp = sTemp.mult(correctTiers(tabSpectreAtt[i], tabSpectreAtt[j], atmos, tabLongueur[i],
                                                    tabLongueur[j]));
                    sProduitCroise = sProduitCroise.sum(sTemp);
                }
            }
 
            s = sCarreModule.sum(sProduitCroise);
 
            //      // 2. au chemins "direct"
            for (i = firstReflex; i < _chemins.size(); i++)
            {
                sCarreModuleDirect = sCarreModuleDirect.sum(tabSpectreAtt[i].mult(tabSpectreAtt[i]));
 
                // Produit croise avec les chemins directs
                for (j = 0; j < _cheminsDirect.size(); j++)
                {
                    sTemp = tabSpectreAtt[i].mult(tabSpectreAttDirect[j].mult(2.0));
                    sTemp = sTemp.mult(correctTiers(tabSpectreAtt[i], tabSpectreAttDirect[j], atmos,
                                                    tabLongueur[i], tabLongueurDirect[j]));
                    sProduitCroiseDirect = sProduitCroiseDirect.sum(sTemp);
                }
                // Produit croise avec les autres chemins reflechis
                for (j = i + 1; j < _chemins.size(); j++)
                {
                    sTemp = tabSpectreAtt[i].mult(tabSpectreAtt[j].mult(2.0));
                    sTemp = sTemp.mult(correctTiers(tabSpectreAtt[i], tabSpectreAtt[j], atmos, tabLongueur[i],
                                                    tabLongueur[j]));
                    sProduitCroiseDirect = sProduitCroiseDirect.sum(sTemp);
                }
            }
 
            attDirect = sCarreModuleDirect.sum(sProduitCroiseDirect); //.abs();
        }
 
        // On remplace la contribution du trajet direct pour toutes les frequences ou cela est necessaire
        for (i = 0; i < indiceDebutEffetEcran; i++)
        {
            s.getTabValReel()[i] = attDirect.getTabValReel()[i];
        }
    }
    build_tab_rays();
    _chemins.clear(); // On efface le tableau des chemins pour (essayer de) gagner de la place en memoire
    _cheminsDirect.clear();
 
    return s;
}
 
OSpectre TYTrajetDefaultSolver::correctTiers(const OSpectreComplex& si, const OSpectreComplex& sj,
                                             const AtmosphericConditions& atmos, const double& ri,
                                             const double& rj) const
{
    const double dp6 = pow(2, (1.0 / 6.0));
    const double invdp6 = 1.0 / dp6;
    const double dfSur2f = (dp6 - invdp6) / 2.0; // df/2f
    OSpectre cosTemp;
    OSpectre s;
 
    OSpectre sTemp = atmos.get_k().mult(ri - rj); // k(ri-rj)
 
    if (ri == rj)
    {
        s = (si.getPhase().subst(sj.getPhase()).subst(sTemp)).cos(); // cos(EPS_i - EPS_j)
    }
    else
    {
 
        s = si.getPhase().subst(sj.getPhase()); // thetaI - thetaJ
 
        double df = sqrt(1 + dfSur2f * dfSur2f);
        cosTemp = sTemp.mult(df);       // k(ri-rj) * sqrt(1 + (df/2f)²)
        cosTemp = cosTemp.sum(s);       // k(ri-rj) * sqrt(1 + (df/2f)²) + (thetaI - thetaJ)
        cosTemp = cosTemp.subst(sTemp); // k(ri-rj)*sqrt(1 + (df/2f)²) + EPS_i - EPS_j
        cosTemp = cosTemp.cos();        // cos(k(ri-rj)*sqrt(1 + (df/2f)²) + EPS_i - EPS_j)
 
        sTemp = sTemp.mult(dfSur2f); // k(ri-rj)*df/2f
 
        s = sTemp.sin();     // sin(k(ri-rj)*df/2f)
        s = s.div(sTemp);    // sin(k(ri-rj)*df/2f) / (k(ri-rj)*df/2f)
        s = s.mult(cosTemp); // sin(k(ri-rj)*df/2f) / (k(ri-rj)*df/2f) * cos(k(ri-rj)*sqrt(1 + (df/2f)²) +
                             // EPS_i - EPS_j)
    }
 
    return s;
}
 
void TYTrajetDefaultSolver::build_tab_rays()
{
    _tabRays.clear();
    for (size_t i = 0; i < _chemins.size(); i++)
    {
        _tabRays.push_back(_chemins[i].get_ray(_ptR));
    }
}