56 #include "eos_compose_fit.h" 58 #include "utilitaires.h" 72 Eos(
"EoS fitted from CompOSE data"), p_eos_low(nullptr), p_eos_high(nullptr),
73 mg(nullptr), mp(nullptr), log_p(nullptr), log_e(nullptr), log_nb(nullptr),
76 ifstream para_file(par_file_name) ;
77 para_file.ignore(1000,
'\n') ;
78 para_file.ignore(1000,
'\n') ;
87 Eos(fich), p_eos_low(nullptr), p_eos_high(nullptr), mg(nullptr),
88 mp(nullptr), log_p(nullptr), log_e(nullptr), log_nb(nullptr), log_cs2(nullptr)
91 char tmp_string[160] ;
92 fread(tmp_string,
sizeof(
char), 160, fich) ;
100 fread_be(&hmax,
sizeof(
double), 1, fich) ;
101 double gamma_low, gamma_high, kappa_low, kappa_high;
102 double m0_low, m0_high, mu0_low, mu0_high ;
103 fread_be(&gamma_low,
sizeof(
double), 1, fich) ;
104 fread_be(&gamma_high,
sizeof(
double), 1, fich) ;
105 fread_be(&kappa_low,
sizeof(
double), 1, fich) ;
106 fread_be(&kappa_high,
sizeof(
double), 1, fich) ;
107 fread_be(&m0_low,
sizeof(
double), 1, fich) ;
108 fread_be(&m0_high,
sizeof(
double), 1, fich) ;
109 fread_be(&mu0_low,
sizeof(
double), 1, fich) ;
110 fread_be(&mu0_high,
sizeof(
double), 1, fich) ;
130 Eos(para_file), p_eos_low(nullptr), p_eos_high(nullptr),
131 mg(nullptr), mp(nullptr), log_p(nullptr), log_e(nullptr), log_nb(nullptr),
145 if (
mg !=
nullptr)
delete mg ;
146 if (
mp !=
nullptr)
delete mp ;
163 cout <<
"The second EOS is not of type Eos_compose_fit !" << endl ;
182 char tmp_string[160] ;
184 fwrite(tmp_string,
sizeof(
char), 160, fich) ;
190 fwrite_be(&hmax,
sizeof(
double), 1, fich) ;
199 fwrite_be(&gamma_low,
sizeof(
double), 1, fich) ;
200 fwrite_be(&gamma_high,
sizeof(
double), 1, fich) ;
201 fwrite_be(&kappa_low,
sizeof(
double), 1, fich) ;
202 fwrite_be(&kappa_high,
sizeof(
double), 1, fich) ;
203 fwrite_be(&m0_low,
sizeof(
double), 1, fich) ;
204 fwrite_be(&m0_high,
sizeof(
double), 1, fich) ;
205 fwrite_be(&mu0_low,
sizeof(
double), 1, fich) ;
206 fwrite_be(&mu0_high,
sizeof(
double), 1, fich) ;
220 ost <<
"EOS of class Eos_compose_fit." << endl ;
221 ost <<
"Built from files in directory " <<
tablename << endl ;
222 ost <<
"Number of coefficients used for the polynomial fit : " 225 <<
", nb_max = " <<
nb_max <<
" [fm^-3]" << endl ;
227 <<
", hmax = " << hmax << endl ;
228 ost <<
"EoS for low density part: " << *
p_eos_low ;
229 ost <<
"EoS for high density part: " << *
p_eos_high ;
241 cout << para_file.good() << endl ;
244 para_file.ignore(1000,
'\n') ;
247 para_file.ignore(1000,
'\n') ;
248 int nr ; para_file >> nr ;
249 para_file.ignore(1000,
'\n') ;
251 Tbl* logh_read = nullptr ;
252 Tbl* logp_read = nullptr ;
253 Tbl* loge_read = nullptr ;
254 Tbl* lognb_read = nullptr ;
255 Tbl* dlpsdlnb = nullptr ;
260 int i_min = 0 ;
int i_mid = 0 ;
262 Tbl nb_read =
exp((*lognb_read)) ;
264 for (
int i=0; i<nbp; i++) {
265 if (nb_read(i) < 10.*
nb_min) i_min = i ;
266 if (nb_read(i) < 10.*
nb_mid) i_mid = i ;
269 int i_max = nbp - 1 ;
272 double x_min = (*logh_read)(i_min) ;
273 double x_mid = (*logh_read)(i_mid) ;
274 double x_max = (*logh_read)(i_max) ;
282 xtab.
set(0) = x_min ;
283 xtab.
set(1) = x_mid ;
284 xtab.
set(nz) = x_max ;
287 mg =
new Mg3d(nz, nr, 1, 1, SYM, SYM) ;
294 if (logh_read !=
nullptr)
delete logh_read ;
295 if (logp_read !=
nullptr)
delete logp_read ;
296 if (loge_read !=
nullptr)
delete loge_read ;
297 if (lognb_read !=
nullptr)
delete lognb_read ;
298 if (dlpsdlnb !=
nullptr)
delete dlpsdlnb ;
304 cout <<
"Writing the Eos_compose_fit object into a Lorene-format file (" 305 << name_file <<
") ..." << flush ;
307 cerr <<
"Eos_compose_fit::write_lorene_table" << endl ;
308 cerr <<
"The number of lines to be outputted is too small!" << endl ;
309 cerr <<
" nlines = " << nlines << endl ;
310 cerr <<
"Aborting..." << endl ;
313 double rhonuc_cgs = Unites::rhonuc_si * 1.e-3 ;
314 double c2_cgs = Unites::c_si * Unites::c_si * 1.e4 ;
316 ofstream lorene_file(name_file) ;
321 lorene_file <<
"#" << endl ;
322 lorene_file <<
"# Built from Eos_compose_fit object" << endl ;
323 lorene_file <<
"# From the table: " <<
tablename << endl ;
324 lorene_file <<
"#\n#" << endl ;
325 lorene_file << nlines <<
"\t Number of lines" << endl ;
326 lorene_file <<
"#\n#\t n_B [fm^{-3}] rho [g/cm^3] p [dyn/cm^2]" << endl ;
327 lorene_file <<
"#" << endl ;
334 double xmin0 =
log(1.e-14) ;
336 int nlines0 = nlines / 10 ;
337 double dx = (xmax0 - xmin0)/
double(nlines0-1) ;
339 double logh = xmin0 ;
340 for (
int i=0; i<nlines0; i++) {
341 double ent =
exp(logh) ;
342 lorene_file << setprecision(16) ;
343 lorene_file << i <<
'\t' <<
nbar_ent_p(ent)/10. <<
'\t' 344 <<
ener_ent_p(ent)*rhonuc_cgs / c2_cgs <<
'\t' 350 xmax0 =
log(10.*hmax) ;
351 dx = (xmax0 - xmin0) /
double(nlines - nlines0-1) ;
352 for (
int i=nlines0; i<nlines; i++) {
353 double ent =
exp(logh) ;
354 lorene_file << setprecision(16) ;
355 lorene_file << i <<
'\t' <<
nbar_ent_p(ent)/10. <<
'\t' 356 <<
ener_ent_p(ent)*rhonuc_cgs / c2_cgs <<
'\t' 360 lorene_file.close() ;
361 cout <<
" done!" << endl ;
378 double logh0 =
log( ent ) ;
399 double logh0 =
log( ent ) ;
420 double logh0 =
log( ent ) ;
452 return (ener + press)*ent / (ener * cs2) ;
464 return ent*(ener + press)/press ;
477 return dlpsdlh0 / dlnsdlh0 ;
485 double logh0 =
log( ent ) ;
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
Higher bound in density.
void read_compose_data(int &nbp, Tbl *&logh, Tbl *&logp, Tbl *&loge, Tbl *&lognb, Tbl *&gam1)
Reads Compose data and stores the values of thermodynamic quantities.
virtual void sauve(FILE *) const
Save in a file.
Cmp log(const Cmp &)
Neperian logarithm.
Cmp exp(const Cmp &)
Exponential.
Eos_compose_fit(const string ¶m_file)
Constructor from a parameter file.
virtual double der_press_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.
const Mg3d * mg
Multi-grid defining the number of domains and spectral coefficients.
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.
Tensor field of valence 0 (or component of a tensorial field).
virtual double der_nbar_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
virtual bool operator!=(const Eos &) const
Comparison operator (difference)
string tablename
Name of the file containing the tabulated data.
virtual double press_ent_p(double ent, const Param *par=0x0) const
Computes the pressure from the log-enthalpy.
virtual double der_press_nbar_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
virtual double csound_square_ent_p(double, const Param *par=0x0) const
Computes the sound speed squared from the enthapy with extra parameters (virtual function implemente...
virtual double val_r_jk(int l, double xi, int j, int k) const
Returns the value of the radial coordinate r for a given and a given collocation point in in a give...
double get_gam() const
Returns the adiabatic index (cf. Eq. (3))
void read_and_compute(ifstream &)
Reads the Compose data and makes the fit.
virtual void sauve(FILE *) const
Save in a file.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
virtual double ener_ent_p(double ent, const Param *par=0x0) const
Computes the total energy density from the log-enthalpy.
void adiabatic_index_fit(int i_min, int i_max, const Tbl &logh_read, const Tbl &logp_read, const Tbl &loge_read, const Tbl &lognb_read, const Tbl &gam1_read)
From the read values, makes the fit on the adiabatic index and deduces the other quantities from ther...
virtual void sauve(FILE *) const
Save in a file.
Scalar * log_cs2
Table of .
void sauve(FILE *fd, bool save_base=false) const
Saves into a file.
virtual ostream & operator>>(ostream &) const
Operator >>
double get_kap() const
Returns the pressure coefficient (cf.
double get_m_0() const
Return the individual particule mass (cf.
Scalar * log_nb
Table of .
int n_coefs
Number of coeficients for polynomial regression.
const Map_af * mp
Mapping in .
double nb_max
Higher bound on the density, above which the EoS is assumed to be a polytrope.
virtual double nbar_ent_p(double ent, const Param *par=0x0) const
Computes the baryon density from the log-enthalpy.
Polytropic equation of state (relativistic case).
const Eos_poly * p_eos_low
Pointer on a polytropic EoS for the description of low densities (nb<nb_min)
int get_nzone() const
Returns the number of domains.
virtual double csound_square_ent_p(double ent, const Param *par=0x0) const
Computes the sound speed squared from the enthapy with extra parameters (virtual function implemente...
int fwrite_be(const int *aa, int size, int nb, FILE *fich)
Writes integer(s) into a binary file according to the big endian convention.
virtual double press_ent_p(double ent, const Param *par=0x0) const
Computes the pressure from the log-enthalpy.
double get_mu_0() const
Return the relativistic chemical potential at zero pressure [unit: , with ].
const Eos_poly * p_eos_high
Pointer on a polytropic EoS for the description of high densities (nb>nb_max)
virtual double der_ener_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
double val_point(double r, double theta, double phi) const
Computes the value of the field at an arbitrary point , by means of the spectral expansion.
double nb_min
Lower bound in baryon density, below which the EoS is assumed to be a polytrope.
int fread_be(int *aa, int size, int nb, FILE *fich)
Reads integer(s) from a binary file according to the big endian convention.
double hmin
Values of enthalpy corresponding to nb_min and nb_max.
double nb_mid
Middle bound in baryon density, above which which the EoS is determined from the polynomial fit to th...
void write_lorene_table(const string &, int nlines=200) const
Save into a table in Lorene format.
virtual bool operator==(const Eos &) const
Comparison operator (egality)
virtual int identify() const
Returns a number to identify the sub-classe of Eos the object belongs to.
virtual ~Eos_compose_fit()
Destructor.