30char eos_tabul_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Eos/eos_tabul.C,v 1.16 2015/08/04 14:41:29 j_novak Exp $" ;
117#include "utilitaires.h"
122void interpol_herm(
const Tbl& ,
const Tbl&,
const Tbl&,
double,
int&,
125void interpol_linear(
const Tbl&,
const Tbl&,
double,
int&,
double&) ;
134 const char* path) :
Eos(name_i)
159 char tmp_string[160] ;
160 fread(tmp_string,
sizeof(
char), 160, fich) ;
192 logh(0x0), logp(0x0), dlpsdlh(0x0),
193 lognb(0x0), dlpsdlnb(0x0)
217 char tmp_string[160] ;
219 fwrite(tmp_string,
sizeof(
char), 160, fich) ;
233 cerr <<
"Eos_tabul::read_table(): " << endl ;
234 cerr <<
"Problem in opening the EOS file!" << endl ;
235 cerr <<
"While trying to open " <<
tablename << endl ;
236 cerr <<
"Aborting..." << endl ;
240 fich.ignore(1000,
'\n') ;
243 for (
int i=0; i<3; i++) {
244 fich.ignore(1000,
'\n') ;
248 fich >> nbp ; fich.ignore(1000,
'\n') ;
250 cerr <<
"Eos_tabul::read_table(): " << endl ;
251 cerr <<
"Wrong value for the number of lines!" << endl ;
252 cerr <<
"nbp = " << nbp << endl ;
253 cerr <<
"Aborting..." << endl ;
257 for (
int i=0; i<3; i++) {
258 fich.ignore(1000,
'\n') ;
261 press =
new double[nbp] ;
262 nb =
new double[nbp] ;
263 ro =
new double[nbp] ;
277 double rhonuc_cgs = rhonuc_si * 1e-3 ;
278 double c2_cgs = c_si * c_si * 1e4 ;
281 double nb_fm3, rho_cgs, p_cgs ;
283 for (
int i=0; i<nbp; i++) {
288 fich >> p_cgs ; fich.ignore(1000,
'\n') ;
289 if ( (nb_fm3<0) || (rho_cgs<0) || (p_cgs < 0) ){
290 cout <<
"Eos_tabul::read_table(): " << endl ;
291 cout <<
"Negative value in table!" << endl ;
292 cout <<
"nb = " << nb_fm3 <<
", rho = " << rho_cgs <<
293 ", p = " << p_cgs << endl ;
294 cout <<
"Aborting..." << endl ;
298 press[i] = p_cgs / c2_cgs ;
304 for (
int i=0; i<nbp; i++) {
305 double h =
log( (ro[i] + press[i]) /
306 (10 * nb[i] * rhonuc_cgs) ) ;
308 if (i==0) { ww = h ; }
309 h = h - ww + 1.e-14 ;
313 dlpsdlh->
set(i) = h * (ro[i] + press[i]) / press[i] ;
317 double p0, p1, p2, n0, n1, n2, dpdnb;
329 dpdnb = p0*(2*n0-n1-n2)/(n0-n1)/(n0-n2) +
330 p1*(n0-n2)/(n1-n0)/(n1-n2) +
331 p2*(n0-n1)/(n2-n0)/(n2-n1) ;
335 for(
int i=1;i<nbp-1;i++) {
337 p0 =
log(press[i-1]);
339 p2 =
log(press[i+1]);
345 dpdnb = p0*(n1-n2)/(n0-n1)/(n0-n2) +
346 p1*(2*n1-n0-n2)/(n1-n0)/(n1-n2) +
347 p2*(n1-n0)/(n2-n0)/(n2-n1) ;
355 p0 =
log(press[nbp-3]);
356 p1 =
log(press[nbp-2]);
357 p2 =
log(press[nbp-1]);
363 dpdnb = p0*(n2-n1)/(n0-n1)/(n0-n2) +
364 p1*(n2-n0)/(n1-n0)/(n1-n2) +
365 p2*(2*n2-n0-n1)/(n2-n0)/(n2-n1) ;
387 cout <<
"hmin, hmax : " <<
hmin <<
" " <<
hmax << endl ;
412 cout <<
"Eos_tabul::nbar_ent_p : ent > hmax !" << endl ;
415 double logh0 =
log10( ent ) ;
421 double pp =
pow(
double(10), logp0) ;
423 double resu = pp / ent * dlpsdlh0 *
exp(-ent) ;
440 cout <<
"Eos_tabul::ener_ent_p : ent > hmax !" << endl ;
443 double logh0 =
log10( ent ) ;
449 double pp =
pow(
double(10), logp0) ;
451 double resu = pp / ent * dlpsdlh0 - pp ;
468 cout <<
"Eos_tabul::press_ent_p : ent > hmax !" << endl ;
472 double logh0 =
log10( ent ) ;
477 return pow(
double(10), logp0) ;
491 cout <<
"Eos_tabul::der_nbar_ent_p : ent > hmax !" << endl ;
515 cout <<
"Eos_tabul::der_ener_ent_p : ent > hmax !" << endl ;
536 cout <<
"Eos_tabul::der_press_ent_p : ent > hmax !" << endl ;
540 double logh0 =
log10( ent ) ;
566 cout <<
"Eos_tabul::der_press_nbar_p : ent > hmax !" << endl ;
570 double logh0 =
log10( ent ) ;
573 interpol_linear(*
logh, *
dlpsdlnb, logh0, i_near, dlpsdlnb0) ;
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 der_press_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
void read_table()
Reads the file containing the table and initializes in the arrays logh , logp and dlpsdlh .
virtual ~Eos_tabul()
Destructor.
Eos_tabul(const char *name_i, const char *table, const char *path)
Standard constructor.
virtual double nbar_ent_p(double ent, const Param *par=0x0) const
Computes the baryon density from the log-enthalpy.
virtual double der_nbar_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
virtual double der_ener_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
virtual void sauve(FILE *) const
Save in a file.
string tablename
Name of the file containing the tabulated data.
double hmin
Lower boundary of the enthalpy interval.
virtual double der_press_nbar_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
string authors
Authors - reference for the table.
double hmax
Upper boundary of the enthalpy interval.
Equation of state base class.
virtual void sauve(FILE *) const
Save in a file.
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.
Standard units of space, time and mass.