147 #include "utilitaires.h" 152 void interpol_herm(
const Tbl& ,
const Tbl&,
const Tbl&,
double,
int&,
155 void interpol_linear(
const Tbl&,
const Tbl&,
double,
int&,
double&) ;
157 void interpol_quad(
const Tbl&,
const Tbl&,
double,
int&,
double&) ;
168 const char* path) :
Eos(name_i)
193 char tmp_string[160] ;
194 fread(tmp_string,
sizeof(
char), 160, fich) ;
226 logh(0x0), logp(0x0), dlpsdlh(0x0),
227 lognb(0x0), dlpsdlnb(0x0), log_cs2(0x0)
252 char tmp_string[160] ;
254 fwrite(tmp_string,
sizeof(
char), 160, fich) ;
265 ifstream is_meos(
"meos_is_being_built.d") ;
270 cerr <<
"Eos_tabul::read_table(): " << endl ;
271 cerr <<
"Problem in opening the EOS file!" << endl ;
272 cerr <<
"While trying to open " <<
tablename << endl ;
273 cerr <<
"Aborting..." << endl ;
277 fich.ignore(1000,
'\n') ;
280 for (
int i=0; i<3; i++) {
281 fich.ignore(1000,
'\n') ;
285 fich >> nbp ; fich.ignore(1000,
'\n') ;
287 cerr <<
"Eos_tabul::read_table(): " << endl ;
288 cerr <<
"Wrong value for the number of lines!" << endl ;
289 cerr <<
"nbp = " << nbp << endl ;
290 cerr <<
"Aborting..." << endl ;
294 for (
int i=0; i<3; i++) {
295 fich.ignore(1000,
'\n') ;
298 press =
new double[nbp] ;
299 nb =
new double[nbp] ;
300 ro =
new double[nbp] ;
318 double rhonuc_cgs = rhonuc_si * 1e-3 ;
319 double c2_cgs = c_si * c_si * 1e4 ;
322 double nb_fm3, rho_cgs, p_cgs ;
324 for (
int i=0; i<nbp; i++) {
329 fich >> p_cgs ; fich.ignore(1000,
'\n') ;
330 if ( (nb_fm3<0) || (rho_cgs<0) || (p_cgs < 0) ){
331 cout <<
"Eos_tabul::read_table(): " << endl ;
332 cout <<
"Negative value in table!" << endl ;
333 cout <<
"nb = " << nb_fm3 <<
", rho = " << rho_cgs <<
334 ", p = " << p_cgs << endl ;
335 cout <<
"Aborting..." << endl ;
339 press[i] = p_cgs / c2_cgs ;
343 press_tbl.
set(i) = press[i] ;
344 nb_tbl.
set(i) = nb[i] ;
345 ro_tbl.
set(i) = ro[i] ;
349 bool store_offset = false ;
350 bool compute_offset = true ;
355 if (temp.find(ok) != string::npos) {
357 compute_offset = false ;
361 store_offset = true ;
365 for (
int i=0; i<nbp; i++) {
366 double h =
log( (ro[i] + press[i]) /
367 (10 * nb[i] * rhonuc_cgs) ) ;
369 if ((i==0) && compute_offset) { ww = h ; }
370 h = h - ww + 1.e-14 ;
374 dlpsdlh->
set(i) = h * (ro[i] + press[i]) / press[i] ;
378 if (store_offset ==
true) {
379 ofstream write_meos(
"meos_is_being_built.d", ios_base::app) ;
380 write_meos <<
"ok" << endl ;
381 write_meos << setprecision(16) << ww << endl ;
387 Tbl logp0 = ((*logp) +
log10(rhonuc_cgs)) *
log(10.) ;
394 for (
int i=0; i<nbp; i++) {
396 tmp.
set(i) = - tmp(i) ;
403 cout <<
"hmin, hmax : " <<
hmin <<
" " <<
hmax << endl ;
428 cout <<
"Eos_tabul::nbar_ent_p : ent > hmax !" << endl ;
431 double logh0 =
log10( ent ) ;
437 double pp =
pow(
double(10), logp0) ;
439 double resu = pp / ent * dlpsdlh0 *
exp(-ent) ;
456 cout <<
"Eos_tabul::ener_ent_p : ent > hmax !" << endl ;
459 double logh0 =
log10( ent ) ;
465 double pp =
pow(
double(10), logp0) ;
467 double resu = pp / ent * dlpsdlh0 - pp ;
484 cout <<
"Eos_tabul::press_ent_p : ent > hmax !" << endl ;
488 double logh0 =
log10( ent ) ;
493 return pow(
double(10), logp0) ;
507 cout <<
"Eos_tabul::der_nbar_ent_p : ent > hmax !" << endl ;
531 cout <<
"Eos_tabul::der_ener_ent_p : ent > hmax !" << endl ;
552 cout <<
"Eos_tabul::der_press_ent_p : ent > hmax !" << endl ;
556 double logh0 =
log10( ent ) ;
582 cout <<
"Eos_tabul::der_press_nbar_p : ent > hmax !" << endl ;
586 double logh0 =
log10( ent ) ;
589 interpol_linear(*
logh, *
dlpsdlnb, logh0, i_near, dlpsdlnb0) ;
607 cout <<
"Eos_tabul::csound_square_ent_p : ent>hmax !" << endl ;
610 double log_ent0 =
log10( ent ) ;
613 interpol_linear(*
logh, *
log_cs2, log_ent0, i_near, log_csound0) ;
615 return pow(10., log_csound0) ;
Cmp log(const Cmp &)
Neperian logarithm.
Cmp exp(const Cmp &)
Exponential.
Standard units of space, time and mass.
Equation of state base class.
double & set(int i)
Read/write of a particular element (index i) (1D case)
double hmin
Lower boundary of the enthalpy interval.
virtual void sauve(FILE *) const
Save in a file.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
virtual ~Eos_tabul()
Destructor.
virtual double der_ener_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
string tablename
Name of the file containing the tabulated data.
string authors
Authors - reference for the table.
virtual void sauve(FILE *) const
Save in a file.
virtual double der_press_nbar_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
double hmax
Upper boundary of the enthalpy interval.
virtual double der_nbar_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
void compute_derivative(const Tbl &xx, const Tbl &ff, Tbl &dfdx)
Derives a function defined on an unequally-spaced grid, approximating it by piecewise parabolae...
Cmp pow(const Cmp &, int)
Power .
virtual void read_table()
Reads the file containing the table and initializes in the arrays logh , logp and dlpsdlh ...
Eos_tabul(const char *name_i, const char *table, const char *path)
Standard constructor.
virtual double csound_square_ent_p(double, const Param *) const
Computes the sound speed squared from the enthapy with extra parameters (virtual function implemente...
Cmp log10(const Cmp &)
Basis 10 logarithm.
int get_taille() const
Gives the total size (ie dim.taille)
virtual double der_press_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
virtual double press_ent_p(double ent, const Param *par=0x0) const
Computes the pressure from the log-enthalpy.
virtual double ener_ent_p(double ent, const Param *par=0x0) const
Computes the total energy density from the log-enthalpy.
virtual double nbar_ent_p(double ent, const Param *par=0x0) const
Computes the baryon density from the log-enthalpy.