hadgas_eos.h

Go to the documentation of this file.

/*
 *
 *    Copyright (c) 2016-
 *      SMASH Team
 *
 *    GNU General Public License (GPLv3 or later)
 *
 */
#ifndef SRC_INCLUDE_SMASH_HADGAS_EOS_H_
#define SRC_INCLUDE_SMASH_HADGAS_EOS_H_
 
#include <gsl/gsl_multiroots.h>
#include <gsl/gsl_roots.h>
#include <gsl/gsl_vector.h>
 
#include <array>
#include <string>
#include <vector>
 
#include "constants.h"
#include "particletype.h"
 
namespace smash {
 
// Forward declaration of HadronGasEos - it is used in EosTable
class HadronGasEos;
 
class EosTable {
 public:
  EosTable(double de, double dnb, double dq, size_t n_e, size_t n_b,
           size_t n_q);
  struct table_element {
    double p;
    double T;
    double mub;
    double mus;
    double muq;
  };
  void compile_table(HadronGasEos& eos,
                     const std::string& eos_savefile_name = "hadgas_eos.dat");
  void get(table_element& res, double e, double nb, double nq) const;
 
 private:
  size_t index(size_t ie, size_t inb, size_t inq) const {
    return n_q_ * (ie * n_nb_ + inb) + inq;
  }
  std::vector<table_element> table_;
  double de_;
  double dnb_;
  double dq_;
  size_t n_e_;
  size_t n_nb_;
  size_t n_q_;
};
 
class HadronGasEos {
 public:
  HadronGasEos(bool tabulate, bool account_for_widths);
  ~HadronGasEos();
 
  static double energy_density(double T, double mub, double mus, double muq);
 
  static double density(double T, double mub, double mus, double muq,
                        bool account_for_resonance_widths = false);
 
  static double pressure(double T, double mub, double mus, double muq,
                         bool account_for_resonance_widths = false) {
    return T * density(T, mub, mus, muq, account_for_resonance_widths);
  }
 
  static double net_baryon_density(double T, double mub, double mus, double muq,
                                   bool account_for_resonance_widths = false);
 
  static double net_strange_density(double T, double mub, double mus,
                                    double muq,
                                    bool account_for_resonance_widths = false);
 
  static double net_charge_density(double T, double mub, double mus, double muq,
                                   bool account_for_resonance_widths = false);
 
  static double partial_density(const ParticleType& ptype, double T, double mub,
                                double mus, double muq,
                                bool account_for_resonance_widths = false);
  static double sample_mass_thermal(const ParticleType& ptype, double beta);
  std::array<double, 4> solve_eos(double e, double nb, double ns, double nq,
                                  std::array<double, 4> initial_approximation);
 
  std::array<double, 4> solve_eos(double e, double nb, double ns, double nq) {
    return solve_eos(e, nb, ns, nq, solve_eos_initial_approximation(e, nb, nq));
  }
 
  std::array<double, 4> solve_eos_initial_approximation(double e, double nb,
                                                        double nq);
 
  static double mus_net_strangeness0(double T, double mub, double muq);
 
  void from_table(EosTable::table_element& res, double e, double nb,
                  double nq) const {
    eos_table_.get(res, e, nb, nq);
  }
 
  static bool is_eos_particle(const ParticleType& ptype) {
    return ptype.is_hadron() && ptype.pdgcode().charmness() == 0;
  }
 
  bool is_tabulated() const { return tabulate_; }
 
  bool account_for_resonance_widths() const {
    return account_for_resonance_widths_;
  }
 
 private:
  struct rparams {
    double e;
    double nb;
    double ns;
    double nq;
    bool account_for_width;
  };
 
  struct eparams {
    double edens;
  };
 
  static double scaled_partial_density_auxiliary(double m_over_T,
                                                 double mu_over_T);
  static double scaled_partial_density(const ParticleType& ptype, double beta,
                                       double mub, double mus, double muq,
                                       bool account_for_width = false);
 
  static int set_eos_solver_equations(const gsl_vector* x, void* params,
                                      gsl_vector* f);
 
  static double e_equation(double T, void* params);
 
  std::string print_solver_state(size_t iter) const;
 
  static constexpr double prefactor_ =
      0.5 * M_1_PI * M_1_PI / (hbarc * hbarc * hbarc);
 
  static constexpr double tolerance_ = 1.e-8;
 
  static constexpr size_t n_equations_ = 4;
 
  EosTable eos_table_ = EosTable(1.e-1, 1.e-1, 1.e-1, 90, 90, 90);
 
  gsl_vector* x_;
 
  gsl_multiroot_fsolver* solver_;
 
  const bool tabulate_;
 
  const bool account_for_resonance_widths_;
};
 
}  // namespace smash
 
#endif  // SRC_INCLUDE_SMASH_HADGAS_EOS_H_