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_HADGAS_EOS_H_
#define SRC_INCLUDE_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, size_t n_e, size_t n_b);
  struct table_element {
    double p;
    double T;
    double mub;
    double mus;
  };
  void compile_table(HadronGasEos& eos,
                     const std::string& eos_savefile_name = "hadgas_eos.dat");
  void get(table_element& res, double e, double nb) const;

 private:
  size_t index(size_t ie, size_t inb) const { return ie * n_nb_ + inb; }
  std::vector<table_element> table_;
  double de_;
  double dnb_;
  size_t n_e_;
  size_t n_nb_;
};

class HadronGasEos {
 public:
  HadronGasEos(bool tabulate, bool account_for_widths);
  ~HadronGasEos();

  static double energy_density(double T, double mub, double mus);

  static double density(double T, double mub, double mus,
                        bool account_for_resonance_widths = false);

  static double pressure(double T, double mub, double mus,
                         bool account_for_resonance_widths = false) {
    return T * density(T, mub, mus, account_for_resonance_widths);
  }

  static double net_baryon_density(double T, double mub, double mus,
                                   bool account_for_resonance_widths = false);

  static double net_strange_density(double T, double mub, double mus,
                                    bool account_for_resonance_widths = false);

  static double partial_density(const ParticleType& ptype, double T, double mub,
                                double mus,
                                bool account_for_resonance_widths = false);
  static double sample_mass_thermal(const ParticleType& ptype, double beta);
  std::array<double, 3> solve_eos(double e, double nb, double ns,
                                  std::array<double, 3> initial_approximation);

  std::array<double, 3> solve_eos(double e, double nb, double ns) {
    return solve_eos(e, nb, ns, solve_eos_initial_approximation(e, nb));
  }

  std::array<double, 3> solve_eos_initial_approximation(double e, double nb);

  static double mus_net_strangeness0(double T, double mub);

  void from_table(EosTable::table_element& res, double e, double nb) const {
    eos_table_.get(res, e, nb);
  }

  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;
    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,
                                       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_ = 3;

  EosTable eos_table_ = EosTable(1.e-2, 1.e-2, 900, 900);

  gsl_vector* x_;

  gsl_multiroot_fsolver* solver_;

  const bool tabulate_;

  const bool account_for_resonance_widths_;
};

}  // namespace smash

#endif  // SRC_INCLUDE_HADGAS_EOS_H_