30char eos_consistent_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Eos/eos_consistent.C,v 1.1 2015/08/04 14:41:29 j_novak Exp $ " ;
49#include "utilitaires.h"
53 void interpol_herm(
const Tbl& ,
const Tbl&,
const Tbl&,
double,
int&,
88 string file_nb = files +
".nb" ;
89 string file_thermo = files +
".thermo" ;
91 ifstream in_nb(file_nb.data()) ;
96 in_nb >> index1 >> index2 ;
97 int nbp = index2 - index1 + 1 ;
100 press =
new double[nbp] ;
101 nb =
new double[nbp] ;
102 ro =
new double[nbp] ;
106 double nb_fm3, rho_cgs, p_cgs, p_over_nb_comp, eps_comp ;
110 double rhonuc_cgs = rhonuc_si * 1e-3 ;
111 double c2_cgs = c_si * c_si * 1e4 ;
112 double m_neutron_MeV, m_proton_MeV ;
114 ifstream in_p_rho (file_thermo.data()) ;
115 in_p_rho >> m_neutron_MeV >> m_proton_MeV ;
116 in_p_rho.ignore(1000,
'\n') ;
118 double p_convert = mev_si * 1.e45 * 10. ;
119 double eps_convert = mev_si * 1.e42 / (c_si*c_si) ;
123 for (
int i=0; i<nbp; i++) {
125 in_p_rho >> dummy_n >> dummy_n >> dummy_n >> p_over_nb_comp ;
126 in_p_rho >> dummy_x >> dummy_x >> dummy_x >> dummy_x >> dummy_x >> eps_comp ;
127 in_p_rho.ignore(1000,
'\n') ;
128 p_cgs = p_over_nb_comp * nb_fm3 * p_convert ;
129 rho_cgs = ( eps_comp + 1. ) * m_neutron_MeV * nb_fm3 * eps_convert ;
131 press[i] = p_cgs / c2_cgs ;
138 for (
int i=0; i<nbp; i++) {
139 pp.
set(i) =
log(press[i] / rhonuc_cgs) ;
140 dh.
set(i) = press[i] / (ro[i] + press[i]) ;
145 for (
int i=0; i<nbp; i++) {
187 cout <<
"The second EOS is not of type Eos_consistent !" << endl ;
214 cout <<
"Eos_consistent::nbar_ent_p : ent > hmax !" << endl ;
217 double logh0 =
log10( ent ) ;
223 double pp =
pow(
double(10), logp0) ;
225 double resu = pp / ent * dlpsdlh0 *
exp(-ent) ;
228 double pp_near =
pow(
double(10), (*
logp)(i_near)) ;
229 double ent_near =
pow(
double(10), (*
logh)(i_near)) ;
230 resu = pp_near / ent_near * (*dlpsdlh)(i_near) *
exp(-ent_near) ;
248 cout <<
"Eos_consistent::ener_ent_p : ent > hmax !" << endl ;
251 double logh0 =
log10( ent ) ;
257 double pp =
pow(
double(10), logp0) ;
259 double resu = pp / ent * dlpsdlh0 - pp ;
262 double pp_near =
pow(
double(10), (*
logp)(i_near)) ;
263 double ent_near =
pow(
double(10), (*
logh)(i_near)) ;
264 resu = pp_near / ent_near * (*dlpsdlh)(i_near) - pp_near ;
282 cout <<
"Eos_consistent::press_ent_p : ent > hmax !" << endl ;
286 double logh0 =
log10( ent ) ;
293 double logp_near = (*logp)(i_near) ;
294 double logp_nearp1 = (*logp)(i_near+1) ;
295 double delta = (*logh)(i_near+1) - (*
logh)(i_near) ;
296 logp0 = (logp_nearp1*(logh0 - (*logh)(i_near))
297 - logp_near*(logh0 - (*logh)(i_near+1))) / delta ;
299 return pow(
double(10), logp0) ;
313 ost <<
"EOS of class Eos_consistent." << endl ;
314 ost <<
"Built from file " <<
tablename << endl ;
315 ost <<
"Authors : " <<
authors << endl ;
316 ost <<
"Number of points in file : " <<
logh->
get_taille() << endl ;
317 ost <<
"Table eventually slightly modified to ensure the relation" << endl ;
318 ost <<
"dp = (e+p) dh" << endl ;
Equation of state for the CompOSE database.
virtual bool operator==(const Eos &) const
Comparison operator (egality)
virtual double press_ent_p(double ent, const Param *par=0x0) const
Computes the pressure from the log-enthalpy.
virtual bool operator!=(const Eos &) const
Comparison operator (difference)
virtual ostream & operator>>(ostream &) const
Operator >>
Eos_consistent(const string &files_path)
Constructor from CompOSE data.
virtual double nbar_ent_p(double ent, const Param *par=0x0) const
Computes the baryon density from the log-enthalpy.
virtual int identify() const
Returns a number to identify the sub-classe of Eos the object belongs to.
virtual ~Eos_consistent()
Destructor.
virtual double ener_ent_p(double ent, const Param *par=0x0) const
Computes the total energy density from the log-enthalpy.
string tablename
Name of the file containing the tabulated data.
double hmin
Lower boundary of the enthalpy interval.
string authors
Authors - reference for the table.
double hmax
Upper boundary of the enthalpy interval.
Equation of state base class.
virtual int identify() const =0
Returns a number to identify the sub-classe of Eos the object belongs to.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
double & set(int i)
Read/write of a particular element (index i) (1D case)
int get_taille() const
Gives the total size (ie dim.taille)
Cmp log10(const Cmp &)
Basis 10 logarithm.
Cmp exp(const Cmp &)
Exponential.
Cmp pow(const Cmp &, int)
Power .
Cmp log(const Cmp &)
Neperian logarithm.
Tbl integ1D(const Tbl &xx, const Tbl &ff)
Integrates a function defined on an unequally-spaced grid, approximating it by piece parabolae.
Standard units of space, time and mass.