ThermalModelEVDiagonal.h

Go to the documentation of this file.

/*
 * Thermal-FIST package
 * 
 * Copyright (c) 2014-2019 Volodymyr Vovchenko
 *
 * GNU General Public License (GPLv3 or later)
 */
#ifndef THERMALMODELEVDIAGONAL_H
#define THERMALMODELEVDIAGONAL_H

#include "HRGBase/ThermalModelBase.h"

namespace thermalfist {

  struct EVSolution {
    double T, mu;
    double n, P, en, s;
    double mutil;
    double nid, eid, sid, pid;
    double fav, f2av, f3av, wav, wfav, w2fav, w2av;
    double scv, scvid;
    double dnmuid, dndTid, demuid, deTid, dsmuid, dsTid, dmunT, dPnT, dmuTn, dPTn, denT, deTn, dsignT, dsigTn, dTnsig;
    double dwiddmu, dfdmu, df2dmu, d2widdmu2;
    double dndmu, dndT, dendmu, dendT;
    double skew, kurtosis;
    double Cs;
    double PT4, eT4, sT3;
    double deltaEN, sigmaEN;
  };

  class ThermalModelEVDiagonal : public ThermalModelBase
  {
  public:
    ThermalModelEVDiagonal(ThermalParticleSystem *TPS, const ThermalModelParameters& params = ThermalModelParameters());

    virtual ~ThermalModelEVDiagonal(void);

    void FillVirialEV(const std::vector<double> & vi = std::vector<double>(0));
    
    double ExcludedVolume(int i) const;
    
    double CommonSuppressionFactor();

    
    virtual void DisableMesonMesonVirial() {}
    virtual void DisableMesonMesonAttraction() {}
    virtual void DisableMesonBaryonVirial() {}
    virtual void DisableMesonBaryonAttraction() {}
    virtual void DisableBaryonBaryonVirial() {}
    virtual void DisableBaryonBaryonAttraction() {}
    virtual void DisableBaryonAntiBaryonVirial() {}
    virtual void DisableBaryonAntiBaryonAttraction() {}

    // Override functions begin

    void SetRadius(double rad);
    
    void SetRadius(int i, double rad);

    void FillVirial(const std::vector<double> & ri = std::vector<double>(0));

    virtual void ReadInteractionParameters(const std::string &filename);

    virtual void WriteInteractionParameters(const std::string &filename);

    virtual double CalculateEigenvolumeFraction();

    double VirialCoefficient(int i, int j) const;

    void SetVirial(int i, int j, double b);

    virtual void ChangeTPS(ThermalParticleSystem *TPS);

    virtual void CalculatePrimordialDensities();
    
    virtual void CalculateFluctuations();

    virtual void CalculateTwoParticleCorrelations();

    virtual std::vector<double> CalculateChargeFluctuations(const std::vector<double> &chgs, int order = 4);

    virtual double CalculatePressure();

    virtual double CalculateEnergyDensity();

    virtual double CalculateEntropyDensity();

    // Dummy
    virtual double CalculateBaryonMatterEntropyDensity() { return 0.; }

    virtual double CalculateMesonMatterEntropyDensity() { return 0.; }

    // TODO properly with excluded volume
    virtual double ParticleScaledVariance(int /*part*/) { return 1.; }

    // TODO properly with excluded volume
    virtual double ParticleSkewness(int /*part*/) { return 1.; }

    // TODO properly with excluded volume
    virtual double ParticleKurtosis(int /*part*/) { return 1.; }

    virtual double ParticleScalarDensity(int part);

    // Override functions end

  protected:
    virtual void SolvePressure();

    double Pressure(double P);

    double DensityId(int i, double Pressure);

    double PressureId(int i, double Pressure);

    double ScaledVarianceId(int i, double Pressure);

    virtual double MuShift(int i) const;

    std::vector<double> m_densitiesid;             
    std::vector<double> m_densitiesidnoshift;      
    std::vector<double> m_v;                       
    double m_Suppression;                          
    double m_Pressure;                             
    double m_Densityid;
    double m_TotalDensity;
    EVSolution m_sol;                              
  private:
    class BroydenEquationsDEV : public BroydenEquations
    {
    public:
      BroydenEquationsDEV(ThermalModelEVDiagonal *model) : BroydenEquations(), m_THM(model) { m_N = 1; m_mnc = 1.; }
      std::vector<double> Equations(const std::vector<double> &x);
      void SetMnc(double mnc) { m_mnc = mnc; }
    private:
      ThermalModelEVDiagonal *m_THM;
      double m_mnc;
    };

    class BroydenJacobianDEV : public BroydenJacobian
    {
    public:
      BroydenJacobianDEV(ThermalModelEVDiagonal *model) : BroydenJacobian(), m_THM(model) { m_mnc = 1.; }
      std::vector<double> Jacobian(const std::vector<double> &x);
      void SetMnc(double mnc) { m_mnc = mnc; }
    private:
      ThermalModelEVDiagonal *m_THM;
      double m_mnc;
    };

    class BroydenSolutionCriteriumDEV : public Broyden::BroydenSolutionCriterium
    {
    public:
      BroydenSolutionCriteriumDEV(ThermalModelEVDiagonal *model, double relative_error = Broyden::TOL) : Broyden::BroydenSolutionCriterium(relative_error), m_THM(model) { }
      virtual bool IsSolved(const std::vector<double>& x, const std::vector<double>& f, const std::vector<double>& xdelta = std::vector<double>()) const;
    protected:
      ThermalModelEVDiagonal *m_THM;
    };

    class BroydenEquationsDEVOrig : public BroydenEquations
    {
    public:
      BroydenEquationsDEVOrig(ThermalModelEVDiagonal *model) : BroydenEquations(), m_THM(model) { m_N = 1; }
      std::vector<double> Equations(const std::vector<double> &x);
    private:
      ThermalModelEVDiagonal *m_THM;
    };

    class BroydenJacobianDEVOrig : public BroydenJacobian
    {
    public:
      BroydenJacobianDEVOrig(ThermalModelEVDiagonal *model) : BroydenJacobian(), m_THM(model) {}
      std::vector<double> Jacobian(const std::vector<double> &x);
    private:
      ThermalModelEVDiagonal *m_THM;
    };
  };

} // namespace thermalfist

#endif