67 #include "utilitaires.h" 72 void interpol_herm(
const Tbl& ,
const Tbl&,
const Tbl&,
double,
int&,
74 void interpol_linear(
const Tbl&,
const Tbl&,
double,
int&,
double&) ;
112 cout <<
"Computing a thermodynamically - consistent table." << endl ;
116 for (
int i=0; i<5; i++) {
117 tfich.ignore(1000,
'\n') ;
121 tfich >> nbp ; tfich.ignore(1000,
'\n') ;
123 cerr <<
"Eos_consistent::read_table(): " << endl ;
124 cerr <<
"Wrong value for the number of lines!" << endl ;
125 cerr <<
"nbp = " << nbp << endl ;
126 cerr <<
"Aborting..." << endl ;
130 for (
int i=0; i<3; i++) {
131 tfich.ignore(1000,
'\n') ;
134 press =
new double[nbp] ;
135 nb =
new double[nbp] ;
136 ro =
new double[nbp] ;
138 double rhonuc_cgs = rhonuc_si * 1e-3 ;
139 double c2_cgs = c_si * c_si * 1e4 ;
142 double nb_fm3, rho_cgs, p_cgs ;
144 for (
int i=0; i<nbp; i++) {
149 tfich >> p_cgs ; tfich.ignore(1000,
'\n') ;
150 if ( (nb_fm3<0) || (rho_cgs<0) || (p_cgs < 0) ){
151 cout <<
"Eos_consistent(ifstream&): " << endl ;
152 cout <<
"Negative value in table!" << endl ;
153 cout <<
"nb = " << nb_fm3 <<
", rho = " << rho_cgs <<
154 ", p = " << p_cgs << endl ;
155 cout <<
"Aborting..." << endl ;
159 press[i] = p_cgs / c2_cgs ;
166 for (
int i=0; i<nbp; i++) {
167 pp.
set(i) =
log(press[i] / rhonuc_cgs) ;
168 dh.
set(i) = press[i] / (ro[i] + press[i]) ;
173 for (
int i=0; i<nbp; i++) {
196 cout <<
"Computing a thermodynamically - consistent table." << endl ;
201 string file_thermo =
tablename +
".thermo" ;
203 ifstream in_nb(file_nb.data()) ;
208 in_nb >> index1 >> index2 ;
209 int nbp = index2 - index1 + 1 ;
212 press =
new double[nbp] ;
213 nb =
new double[nbp] ;
214 ro =
new double[nbp] ;
218 double nb_fm3, rho_cgs, p_cgs, p_over_nb_comp, eps_comp ;
222 double rhonuc_cgs = rhonuc_si * 1e-3 ;
223 double c2_cgs = c_si * c_si * 1e4 ;
224 double m_neutron_MeV, m_proton_MeV ;
226 ifstream in_p_rho (file_thermo.data()) ;
227 in_p_rho >> m_neutron_MeV >> m_proton_MeV ;
228 in_p_rho.ignore(1000,
'\n') ;
230 double p_convert = mev_si * 1.e45 * 10. ;
231 double eps_convert = mev_si * 1.e42 / (c_si*c_si) ;
235 for (
int i=0; i<nbp; i++) {
237 in_p_rho >> dummy_n >> dummy_n >> dummy_n >> p_over_nb_comp ;
238 in_p_rho >> dummy_x >> dummy_x >> dummy_x >> dummy_x >> dummy_x >> eps_comp ;
239 in_p_rho.ignore(1000,
'\n') ;
240 p_cgs = p_over_nb_comp * nb_fm3 * p_convert ;
241 rho_cgs = ( eps_comp + 1. ) * m_neutron_MeV * nb_fm3 * eps_convert ;
243 press[i] = p_cgs / c2_cgs ;
250 for (
int i=0; i<nbp; i++) {
251 pp.
set(i) =
log(press[i] / rhonuc_cgs) ;
252 dh.
set(i) = press[i] / (ro[i] + press[i]) ;
257 for (
int i=0; i<nbp; i++) {
297 cout <<
"The second EOS is not of type Eos_consistent !" << endl ;
324 cout <<
"Eos_consistent::nbar_ent_p : ent > hmax !" << endl ;
327 double logh0 =
log10( ent ) ;
333 double pp =
pow(
double(10), logp0) ;
335 double resu = pp / ent * dlpsdlh0 *
exp(-ent) ;
338 double pp_near =
pow(
double(10), (*
logp)(i_near)) ;
339 double ent_near =
pow(
double(10), (*
logh)(i_near)) ;
340 resu = pp_near / ent_near * (*dlpsdlh)(i_near) *
exp(-ent_near) ;
358 cout <<
"Eos_consistent::ener_ent_p : ent > hmax !" << endl ;
361 double logh0 =
log10( ent ) ;
367 double pp =
pow(
double(10), logp0) ;
369 double resu = pp / ent * dlpsdlh0 - pp ;
372 double p_near =
pow(
double(10), (*
logp)(i_near)) ;
373 double p_nearp1 =
pow(
double(10), (*
logp)(i_near+1)) ;
374 double ener_near = p_near/
pow(
double(10), (*
logh)(i_near))
375 * (*dlpsdlh)(i_near) - p_near ;
376 double ener_nearp1 = p_nearp1/
pow(
double(10), (*
logh)(i_near+1))
377 * (*dlpsdlh)(i_near+1) - p_nearp1 ;
378 double delta = (*logh)(i_near+1) - (*
logh)(i_near) ;
379 resu = (ener_nearp1*(logh0 - (*logh)(i_near))
380 - ener_near*(logh0 - (*logh)(i_near+1))) / delta ;
398 cout <<
"Eos_consistent::press_ent_p : ent > hmax !" << endl ;
402 double logh0 =
log10( ent ) ;
409 double logp_near = (*logp)(i_near) ;
410 double logp_nearp1 = (*logp)(i_near+1) ;
411 double delta = (*logh)(i_near+1) - (*
logh)(i_near) ;
412 logp0 = (logp_nearp1*(logh0 - (*logh)(i_near))
413 - logp_near*(logh0 - (*logh)(i_near+1))) / delta ;
415 return pow(
double(10), logp0) ;
428 cout <<
"Eos_consistent::csound_square_ent_p : ent>hmax !" << endl ;
431 double log_ent0 =
log10( ent ) ;
434 interpol_linear(*
logh, *
log_cs2, log_ent0, i_near, log_csound0) ;
436 return pow(10., log_csound0) ;
451 ost <<
"EOS of class Eos_consistent." << endl ;
452 ost <<
"Built from file " <<
tablename << endl ;
453 ost <<
"Authors : " <<
authors << endl ;
454 ost <<
"Number of points in file : " <<
logh->
get_taille() << endl ;
455 ost <<
"Table eventually slightly modified to ensure the relation" << endl ;
456 ost <<
"dp = (e+p) dh" << endl ;
virtual ostream & operator>>(ostream &) const
Operator >>
Cmp log(const Cmp &)
Neperian logarithm.
Cmp exp(const Cmp &)
Exponential.
virtual ~Eos_consistent()
Destructor.
Standard units of space, time and mass.
virtual double ener_ent_p(double ent, const Param *par=0x0) const
Computes the total energy density from the log-enthalpy.
Equation of state base class.
double & set(int i)
Read/write of a particular element (index i) (1D case)
virtual int identify() const =0
Returns a number to identify the sub-classe of Eos the object belongs to.
virtual void read_compose_data()
Reads the files containing the table and initializes in the arrays logh , logp and dlpsdlh (CompOSE f...
double hmin
Lower boundary of the enthalpy interval.
virtual void read_table()
Reads the file containing the table and initializes in the arrays logh , logp and dlpsdlh ...
int format
0 for standard (old) LORENE format, 1 for CompOSE format
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
virtual double csound_square_ent_p(double, const Param *) const
Computes the sound speed squared from the enthapy with extra parameters (virtual function implemente...
string tablename
Name of the file containing the tabulated data.
string authors
Authors - reference for the table.
virtual double press_ent_p(double ent, const Param *par=0x0) const
Computes the pressure from the log-enthalpy.
virtual int identify() const
Returns a number to identify the sub-classe of Eos the object belongs to.
double hmax
Upper boundary of the enthalpy interval.
virtual bool operator!=(const Eos &) const
Comparison operator (difference)
Tbl integ1D(const Tbl &xx, const Tbl &ff)
Integrates a function defined on an unequally-spaced grid, approximating it by piecewise parabolae...
Cmp pow(const Cmp &, int)
Power .
virtual double nbar_ent_p(double ent, const Param *par=0x0) const
Computes the baryon density from the log-enthalpy.
Cmp log10(const Cmp &)
Basis 10 logarithm.
int get_taille() const
Gives the total size (ie dim.taille)
virtual bool operator==(const Eos &) const
Comparison operator (egality)
Eos_consistent(const string &files_path)
Constructor from CompOSE data.
Equation of state for the CompOSE database.