ThermalModelEVCrossterms.h

Go to the documentation of this file.

/*
 * Thermal-FIST package
 * 
 * Copyright (c) 2015-2018 Volodymyr Vovchenko
 *
 * GNU General Public License (GPLv3 or later)
 */
#ifndef THERMALMODELEVCROSSTERMS_H
#define THERMALMODELEVCROSSTERMS_H

#include "HRGBase/ThermalModelBase.h"

namespace thermalfist {

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

    virtual ~ThermalModelEVCrossterms(void);

    // Override functions begin

    virtual 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);

    void SetRadius(double rad);

    double VirialCoefficient(int i, int j) const;

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

    virtual void ChangeTPS(ThermalParticleSystem *TPS);
    
    virtual void CalculatePrimordialDensities();

    void CalculateTwoParticleCorrelations();

    void CalculateFluctuations();

    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 m_skewprim[part]; }

    // TODO properly with excluded volume
    virtual double ParticleKurtosis(int part) { return m_kurtprim[part]; }

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

    // Override functions end

  protected:
    virtual void SolvePressure(bool resetPartials = true);  // Using Broyden's method

    void SolvePressureIter();  // Using iteration method
    
    virtual void CalculatePrimordialDensitiesNoReset();

    virtual void CalculatePrimordialDensitiesIter();

    virtual void SolveDiagonal();

    double PartialPressureDiagonal(int i, double P);

    double PressureDiagonalTotal(double P);

    virtual double DensityId(int i, const std::vector<double>& pstars = std::vector<double>());

    virtual double Pressure(int i, const std::vector<double>& pstars = std::vector<double>());

    double ScaledVarianceId(int ind);
    
    virtual double MuShift(int i) const;

    std::vector<double> m_densitiesid;            
    std::vector<double> m_Ps;                     
    std::vector< std::vector<double> > m_Virial;  
    double m_Pressure;                            
    double m_TotalEntropyDensity;                 
  private:
    class BroydenEquationsCRS : public BroydenEquations
    {
    public:
      BroydenEquationsCRS(ThermalModelEVCrossterms *model) : BroydenEquations(), m_THM(model) { m_N = model->Densities().size(); }
      std::vector<double> Equations(const std::vector<double> &x);
    private:
      ThermalModelEVCrossterms *m_THM;
    };

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

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

    class BroydenSolutionCriteriumCRS : public Broyden::BroydenSolutionCriterium
    {
    public:
      BroydenSolutionCriteriumCRS(ThermalModelEVCrossterms *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:
      ThermalModelEVCrossterms *m_THM;
    };
  };

} // namespace thermalfist

#endif