56 #include "eos_compose_fit.h" 58 #include "utilitaires.h" 60 #include "graphique.h" 64 Tbl*& logp_read,
Tbl*& loge_read,
65 Tbl*& lognb_read,
Tbl*& gam1_read) {
70 string file_thermo =
tablename +
"eos.thermo" ;
72 cout <<
"opening " << file_nb <<
" and " << file_thermo <<
" ... " << flush ;
73 ifstream in_nb(file_nb.data()) ;
75 cerr <<
"Eos_compose_fit::read_compose_data:" << endl ;
76 cerr <<
"Problem in opening the EOS file!" << endl ;
77 cerr <<
"While trying to open " << file_nb << endl ;
78 cerr <<
"Aborting..." << endl ;
86 in_nb >> index1 >> index2 ;
87 nbp = index2 - index1 + 1 ;
94 if (logh_read !=
nullptr)
delete logh_read ;
95 logh_read =
new Tbl(nbp) ;
97 if (logp_read !=
nullptr)
delete logp_read ;
98 logp_read =
new Tbl(nbp) ;
100 if (loge_read !=
nullptr)
delete loge_read ;
101 loge_read =
new Tbl(nbp) ;
103 if (lognb_read !=
nullptr)
delete lognb_read ;
104 lognb_read =
new Tbl(nbp) ;
106 if (gam1_read !=
nullptr)
delete gam1_read ;
107 gam1_read =
new Tbl(nbp) ;
112 double nb_fm3, rho_cgs, p_cgs, p_over_nb_comp, eps_comp ;
117 double rhonuc_cgs = Unites::rhonuc_si * 1e-3 ;
118 double c2_cgs = Unites::c_si * Unites::c_si * 1e4 ;
119 double m_neutron_MeV, m_proton_MeV ;
121 ifstream in_p_rho (file_thermo.data()) ;
123 cerr <<
"Reading CompOSE data: " << endl ;
124 cerr <<
"Problem in opening the EOS file!" << endl ;
125 cerr <<
"While trying to open " << file_thermo << endl ;
126 cerr <<
"Aborting..." << endl ;
129 in_p_rho >> m_neutron_MeV >> m_proton_MeV ;
130 in_p_rho.ignore(1000,
'\n') ;
133 double p_convert = Unites::mev_si * 1.e45 * 10. ;
135 double eps_convert = Unites::mev_si * 1.e42 / (Unites::c_si*Unites::c_si) ;
137 cout <<
" done. " << endl ;
138 cout <<
"Reading the table ... " << flush ;
143 for (
int i=0; i<nbp; i++) {
145 in_p_rho >> dummy_n >> dummy_n >> dummy_n >> p_over_nb_comp ;
146 in_p_rho >> dummy_x >> dummy_x >> dummy_x >> dummy_x >> dummy_x
148 in_p_rho.ignore(1000,
'\n') ;
149 p_cgs = p_over_nb_comp * nb_fm3 * p_convert ;
150 rho_cgs = ( eps_comp + 1. ) * m_neutron_MeV * nb_fm3 * eps_convert ;
152 if ( (nb_fm3 < 0) || (rho_cgs < 0) || (p_cgs < 0) ){
153 cerr <<
"Eos_compose_fit::read_compose_data: " << endl ;
154 cerr <<
"Negative value in table!" << endl ;
155 cerr <<
"nb = " << nb_fm3 <<
", rho = " << rho_cgs <<
156 ", p = " << p_cgs << endl ;
157 cerr <<
"Aborting..." << endl ;
161 press.
set(i) = p_cgs / c2_cgs ;
163 ener.
set(i) = rho_cgs ;
166 logp_read->
set(i) =
log( press(i) / rhonuc_cgs ) ;
167 loge_read->
set(i) =
log( ener(i) / rhonuc_cgs ) ;
168 lognb_read->
set(i) =
log( 10.* nb(i) ) ;
171 cout <<
" done." << endl ;
176 cout <<
"Computation of derived quantities and file outputs.. " << flush ;
183 for (
int i=0; i<nbp; i++) {
184 double h =
log( (ener(i) + press(i)) / (10 * nb(i) * rhonuc_cgs) ) ;
186 if (i==0) { ww = h ; }
187 h = h - ww + 1.e-14 ;
189 logh_read->
set(i) =
log( h ) ;
195 cout <<
"done" << endl ;
201 logh_read,
const Tbl& logp_read,
203 lognb_read,
const Tbl& gam1_read) {
205 int nbp = logh_read.
get_dim(0) ;
206 double nb0_max =
exp(lognb_read(nbp-1)) ;
208 assert((
mp !=
nullptr) && (
mg !=
nullptr)) ;
214 for (
int l=1; l<nz; l++)
218 double x_max = (+xx)(nz-1, 0, 0, nr-1) ;
219 double x_mid = (+xx)(1, 0, 0, 0) ;
220 double x_min = (+xx)(0, 0, 0, 0) ;
228 double mean_gam = 0. ;
229 for (
int i=0; i<i_min; i++) {
230 mean_gam += gam1_read(i) ;
241 mean_gam /= double(np_gam) ;
242 cout <<
"Mean Gamma_1 = " << mean_gam << endl ;
246 gam1_spec = mean_gam ;
251 int ndata = nbp - i_mid ;
254 for (
int i=i_mid; i<nbp; i++) {
255 xdata.
set(i-i_mid) = 2.*(logh_read(i) - x_mid)
256 / (x_max - x_mid) - 1. ;
257 ydata.
set(i-i_mid) = gam1_read(i) ;
264 for (
int i=0; i<=
n_coefs; i++) {
281 * (xx - x_min) / (x_mid - x_min) ;
290 double p_min =
exp(logp_read(i_min)) ;
291 double ccc = (mean_gam - 1.) / mean_gam ;
292 double kappa_low0 =
pow(ccc*expm1(
exp(x_min)), mean_gam)
293 /
pow(p_min, mean_gam*ccc) ;
294 double kappa_low1 = 0.1*kappa_low0 ;
297 double logp_max = logp_read(nbp-1) ;
301 while(fabs(logp_max0 - logp_max)/logp_max > 1.e-3) {
305 = (kappa_low0*(logp_max1 -logp_max)- kappa_low1*(logp_max0 - logp_max))
306 / (logp_max1 - logp_max0) ;
307 cout <<
"Kappa(new) = " << kappa_low_new <<
", logp = " << logp_max1 << endl ;
308 if (kappa_low_new < 0.) kappa_low_new = 0.01*kappa_low0 ;
309 kappa_low0 = kappa_low1 ;
310 logp_max0 = logp_max1 ;
311 kappa_low1 = kappa_low_new ;
330 Scalar spec_nbar = spec_press / (YYY * expexpx) ;
331 cout <<
"Relative difference in baryon density at the last point" << endl ;
332 cout <<
"(fitted computed / original tabulated) : " 333 << fabs(1. - spec_nbar.
val_grid_point(nz-1, 0, 0, nr-1) / nb0_max)
347 double gam_high = gam1_spec.val_grid_point(1, 0, 0, nr-1) ;
352 - kappa_high/(gam_high - 1.)
359 cout <<
"... done (analytic representation)" << endl ;
364 (
double kappa_low,
double mean_gam,
const Scalar& gam1_spec,
371 if (p_eos_low !=
nullptr)
delete p_eos_low ;
372 p_eos_low =
new Eos_poly(mean_gam, kappa_low) ;
375 double p_min = p_eos_low->press_ent_p(ent_min) ;
376 double e_min = p_eos_low->ener_ent_p(ent_min) ;
381 Scalar fff = ((gam1_spec-1.)/gam1_spec)*expx*expexpx ;
390 double Y_min = p_min / (e_min + p_min) ;
395 YYY = (A_min + FFF)/expexpx ;
400 Scalar Pidot = expx / YYY ;
401 if (log_p !=
nullptr)
delete log_p ;
404 double lnp0 =
log(p_min) - log_p->val_grid_point(0, 0, 0, 0) ;
405 (*log_p) = (*log_p) + lnp0 ;
408 return log_p->val_grid_point(1, 0, 0, mg->get_nr(1)-1) ;
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.
Cmp log(const Cmp &)
Neperian logarithm.
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Tbl poly_regression(const Tbl &, const Tbl &, int)
Polynomial regression, giving Chebyshev coefficients.
Cmp exp(const Cmp &)
Exponential.
void coef() const
Computes the coeffcients of *this.
const Mg3d * mg
Multi-grid defining the number of domains and spectral coefficients.
const Tbl & domain(int l) const
Read-only of the value in a given domain.
double & set(int i)
Read/write of a particular element (index i) (1D case)
Tensor field of valence 0 (or component of a tensorial field).
void coef_i() const
Computes the physical value of *this.
string tablename
Name of the file containing the tabulated data.
Scalar primr(bool null_infty=true) const
Computes the radial primitive which vanishes for .
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field. ...
Tbl & set_domain(int l)
Read/write of the value in a given domain.
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
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...
Scalar * log_cs2
Table of .
double integrate_equations(double kappa, double mean_gam, const Scalar &gamma1, const Scalar &log_ent, const Scalar &ent, Scalar &YYY)
Integrates the differential system giving all thermo quantities from the adiabatic index...
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.
Mtbl * c
Values of the function at the points of the multi-grid.
Polytropic equation of state (relativistic case).
int get_nzone() const
Returns the number of domains.
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...
int get_dim(int i) const
Gives the i-th dimension (ie dim.dim[i])
Cmp pow(const Cmp &, int)
Power .
const Eos_poly * p_eos_high
Pointer on a polytropic EoS for the description of high densities (nb>nb_max)
Active physical coordinates and mapping derivatives.
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Valeur & set_spectral_va()
Returns va (read/write version)
Coord r
r coordinate centered on the grid