52 void huntm(
const Tbl& xx,
double& x,
int& i_low) ;
55 const Tbl& AA,
const Tbl& N_lim,
const Tbl& Ent_lim,
double gamma,
double kappa,
double n_lim0):
56 gamma_high(0x0), kappa_high(0x0), Lambda_high(0x0), a_high(0x0), n_lim_high(0x0), ent_lim_high(0x0), gamma_low(gamma), kappa_low(kappa), n_lim(n_lim0) {
58 set_name(
"Pseudo-polytropic fit of cold, beta-equilibrated EoS") ;
66 gamma_high =
new Tbl(Gamma) ;
67 n_param_high = gamma_high->get_taille() ;
69 kappa_high =
new Tbl(Kappa) ;
71 Lambda_high =
new Tbl(Lambda) ;
73 a_high =
new Tbl(AA) ;
75 n_lim_high =
new Tbl(N_lim) ;
77 ent_lim_high =
new Tbl(Ent_lim) ;
79 eos_low =
new Eos_poly(gamma_low, kappa_low) ;
86 Lambda_high(eosi.Lambda_high), a_high(eosi.a_high), n_lim_high(eosi.n_lim_high), ent_lim_high(eosi.ent_lim_high), gamma_low(eosi.gamma_low),
87 kappa_low(eosi.kappa_low), n_lim(eosi.n_lim), ent_lim(eosi.ent_lim), n_param_high(eosi.n_param_high) {}
94 fread_be(&n_param_high,
sizeof(
int), 1, fich) ;
95 gamma_high =
new Tbl(fich) ;
96 kappa_high =
new Tbl(fich) ;
97 Lambda_high =
new Tbl(fich) ;
98 a_high =
new Tbl(fich) ;
99 n_lim_high =
new Tbl(fich) ;
100 ent_lim_high =
new Tbl(fich) ;
101 fread_be(&kappa_low,
sizeof(
double), 1, fich) ;
102 fread_be(&gamma_low,
sizeof(
double), 1, fich) ;
103 fread_be(&n_lim,
sizeof(
double), 1, fich) ;
104 fread_be(&ent_lim,
sizeof(
double), 1, fich) ;
106 eos_low =
new Eos_poly(gamma_low, kappa_low) ;
113 fich >> n_param_high ; fich.ignore(1000,
'\n') ;
116 for (
int i=0; i < n_param_high; i++)
117 fich >> gamma_high->
set(i) ;
118 fich.ignore(1000,
'\n') ;
121 for (
int i=0; i < n_param_high; i++)
122 fich >> kappa_high->
set(i) ;
123 fich.ignore(1000,
'\n') ;
126 for (
int i=0; i < n_param_high; i++)
127 fich >> Lambda_high->
set(i) ;
128 fich.ignore(1000,
'\n') ;
131 for (
int i=0; i < n_param_high; i++)
132 fich >> a_high->
set(i) ;
133 fich.ignore(1000,
'\n') ;
136 for (
int i=0; i < n_param_high; i++)
137 fich >> n_lim_high->
set(i) ;
138 fich.ignore(1000,
'\n') ;
141 for (
int i=0; i < n_param_high; i++)
142 fich >> ent_lim_high->
set(i) ;
143 fich.ignore(1000,
'\n') ;
145 n_lim = (*n_lim_high)(0) ;
146 ent_lim = (*ent_lim_high)(0) ;
148 fich >> kappa_low ; fich.ignore(1000,
'\n') ;
150 fich >> gamma_low ; fich.ignore(1000,
'\n') ;
157 eos_low =
new Eos_poly(gamma_low, kappa_low) ;
168 if (gamma_high != 0x0)
delete gamma_high ;
169 if (kappa_high != 0x0)
delete kappa_high ;
170 if (Lambda_high != 0x0)
delete Lambda_high ;
171 if (a_high != 0x0)
delete a_high ;
172 if (n_lim_high != 0x0)
delete n_lim_high ;
173 if (ent_lim_high != 0x0)
delete ent_lim_high ;
174 if (eos_low != 0x0)
delete eos_low ;
185 gamma_high = eosi.gamma_high ;
186 kappa_high = eosi.kappa_high ;
187 Lambda_high = eosi.Lambda_high ;
188 a_high = eosi.a_high ;
189 n_lim_high = eosi.n_lim_high ;
190 ent_lim_high = eosi.ent_lim_high ;
192 n_param_high = eosi.n_param_high ;
193 gamma_low = eosi.gamma_low ;
194 kappa_low = eosi.kappa_low ;
195 ent_lim = eosi.ent_lim ;
196 eos_low = eosi.eos_low ;
208 cout <<
"The second EOS is not of type Piecewise_polytrope_1D !" << endl ;
215 for (
int i=0; i < n_param_high; i++){
216 if ((*eos.gamma_high)(i) != (*gamma_high)(i)) {
218 <<
"The two Piecewise_polytrope_1D have different coefficients: " << (*gamma_high)(i) <<
" <-> " 219 << (*eos.gamma_high)(i) << endl ;
223 if ((*eos.kappa_high)(i) != (*kappa_high)(i)) {
225 <<
"The two Piecewise_polytrope_1D have different coefficients: " << (*kappa_high)(i) <<
" <-> " 226 << (*eos.kappa_high)(i) << endl ;
230 if ((*eos.Lambda_high)(i) != (*Lambda_high)(i)) {
232 <<
"The two Piecewise_polytrope_1D have different coefficients: " << (*Lambda_high)(i) <<
" <-> " 233 << (*eos.Lambda_high)(i) << endl ;
237 if ((*eos.a_high)(i) != (*a_high)(i)) {
239 <<
"The two Piecewise_polytrope_1D have different coefficients: " << (*a_high)(i) <<
" <-> " 240 << (*eos.a_high)(i) << endl ;
244 if ((*eos.n_lim_high)(i) != (*n_lim_high)(i)) {
246 <<
"The two Piecewise_polytrope_1D have different coefficients: " << (*n_lim_high)(i) <<
" <-> " 247 << (*eos.n_lim_high)(i) << endl ;
251 if ((*eos.ent_lim_high)(i) != (*ent_lim_high)(i)) {
253 <<
"The two Piecewise_polytrope_1D have different coefficients: " << (*ent_lim_high)(i) <<
" <-> " 254 << (*eos.ent_lim_high)(i) << endl ;
259 if (eos.n_lim != n_lim) {
261 <<
"The two Piecewise_polytrope_1D have different limiting densities n_lim: " << n_lim <<
" <-> " 262 << eos.n_lim << endl ;
266 if (eos.ent_lim != ent_lim) {
268 <<
"The two Piecewise_polytrope_1D have different limiting enthalpies ent_lim: " << ent_lim <<
" <-> " 269 << eos.ent_lim << endl ;
273 if (eos.n_param_high != n_param_high) {
275 <<
"The two Piecewise_polytrope_1D have different number of coefficients: " << n_param_high <<
" <-> " 276 << eos.n_param_high << endl ;
280 if (eos.kappa_low != kappa_low) {
282 <<
"The two Piecewise_polytrope_1D have different kappa in low polytrope: " << kappa_low <<
" <-> " 283 << eos.kappa_low << endl ;
287 if (eos.gamma_low != gamma_low) {
289 <<
"The two Piecewise_polytrope_1D have different kappa in low polytrope: " << gamma_low <<
" <-> " 290 << eos.gamma_low << endl ;
315 fwrite_be(&n_param_high,
sizeof(
int), 1, fich) ;
316 gamma_high->
sauve(fich) ;
317 kappa_high->
sauve(fich) ;
318 Lambda_high->
sauve(fich) ;
319 a_high->
sauve(fich) ;
320 n_lim_high->
sauve(fich) ;
321 ent_lim_high->
sauve(fich) ;
322 fwrite_be(&kappa_low,
sizeof(
double), 1, fich) ;
323 fwrite_be(&gamma_low,
sizeof(
double), 1, fich) ;
324 fwrite_be(&n_lim,
sizeof(
double), 1, fich) ;
325 fwrite_be(&ent_lim,
sizeof(
double), 1, fich) ;
331 ost << setprecision(16) <<
"EOS of class Piecewise_polytrope_1D (analytical fit of cold EoS): " <<
'\n' ;
332 ost <<
" Gamma_high coefficients: " << *gamma_high <<
'\n' ;
333 ost <<
" Kappa_high coefficients: " << *kappa_high <<
'\n' ;
334 ost <<
" Lambda_high coefficients: " << *Lambda_high <<
'\n' ;
335 ost <<
" a_high coefficients: " << *a_high <<
'\n' ;
336 ost <<
" n_lim_high coefficients: " << *n_lim_high <<
'\n' ;
337 ost <<
" ent_lim_high coefficients: " << *ent_lim_high <<
'\n' ;
338 ost <<
"Low densities polytrope :" <<
'\n' ;
339 cout << *eos_low <<
'\n' ;
340 ost << setprecision(16) <<
" Limiting density n_lim: " << 0.1*n_lim <<
" [fm^-3]" <<
'\n' ;
341 ost << setprecision(16) <<
" Limiting enthalpy h_lim: " << ent_lim <<
" [c^2]" << endl;
357 if (ent > ent_lim ) {
360 huntm(*ent_lim_high, ent, i_low) ;
361 Tbl& Gamma = *gamma_high ;
362 Tbl& Kappa = *kappa_high ;
364 double nn = 1./(Gamma(i_low)-1.) ;
366 return pow((
exp(ent)-1.-AA(i_low))/Kappa(i_low)/(nn+1.), nn) ;
369 else if ( (0 < ent) && (ent < ent_lim)){
382 if (ent > ent_lim ) {
385 huntm(*ent_lim_high, ent, i_low) ;
386 Tbl& Gamma = *gamma_high ;
387 Tbl& Kappa = *kappa_high ;
389 Tbl& Lambda = *Lambda_high ;
390 double nn = 1./(Gamma(i_low)-1.) ;
391 double nbar =
pow((
exp(ent)-1.-AA(i_low))/Kappa(i_low)/(nn+1.), nn) ;
394 return Kappa(i_low)*nn *
pow(nbar, Gamma(i_low)) + (1.+AA(i_low))*nbar - Lambda(i_low) ;
396 else if ( (0 < ent) && (ent < ent_lim)){
409 if (ent > ent_lim ) {
412 huntm(*ent_lim_high, ent, i_low) ;
413 Tbl& Gamma = *gamma_high ;
414 Tbl& Kappa = *kappa_high ;
416 Tbl& Lambda = *Lambda_high ;
417 double nn = 1./(Gamma(i_low)-1.) ;
418 double nbar =
pow((
exp(ent)-1.-AA(i_low))/Kappa(i_low)/(nn+1.), nn) ;
420 return Kappa(i_low) *
pow(nbar, Gamma(i_low)) + Lambda(i_low) ;
422 else if ( (0 < ent) && (ent < ent_lim)){
456 if (ent > ent_lim ) {
459 huntm(*ent_lim_high, ent, i_low) ;
461 Tbl& Gamma = *gamma_high ;
463 return (Gamma(i_low)-1.)*(1.-
exp(-ent)*(1.+AA(i_low))) ;
465 else if ( (0 < ent) && (ent < ent_lim)){
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 specific enthalpy.
virtual double ener_ent_p(double ent, const Param *par=0x0) const
Computes the total energy density from the specific enthalpy.
virtual bool operator!=(const Eos &) const
Comparison operator (difference)
Cmp exp(const Cmp &)
Exponential.
Equation of state base class.
void operator=(const Piecewise_polytrope_1D &)
Assignment to another Piecewise_polytrope_1D.
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 sauve(FILE *) const
Save in a file.
virtual double der_ener_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the specific enthalpy.
Piecewise_polytrope_1D(const Tbl &, const Tbl &, const Tbl &, const Tbl &, const Tbl &, const Tbl &, double, double, double)
Standard constructor.
virtual void sauve(FILE *) const
Save in a file.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
int get_ndim() const
Gives the number of dimensions (ie dim.ndim)
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).
virtual double press_ent_p(double ent, const Param *par=0x0) const
Computes the pressure from the specific enthalpy.
virtual double csound_square_ent_p(double, const Param *) const
Computes the sound speed squared from the enthapy with extra parameters (virtual function implemente...
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...
virtual double der_nbar_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the specific enthalpy.
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 bool operator==(const Eos &) const
Comparison operator (egality)
virtual double press_ent_p(double ent, const Param *par=0x0) const
Computes the pressure from the log-enthalpy.
Cmp pow(const Cmp &, int)
Power .
int fread_be(int *aa, int size, int nb, FILE *fich)
Reads integer(s) from a binary file according to the big endian convention.
virtual double der_press_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the specific enthalpy.
void c_est_pas_fait(const char *)
Helpful function to say something is not implemented yet.
void set_name(const char *name_i)
Sets the EOS name.
int get_taille() const
Gives the total size (ie dim.taille)
void sauve(FILE *) const
Save in a file.
virtual ~Piecewise_polytrope_1D()
Destructor.
virtual ostream & operator>>(ostream &) const
Operator >>
virtual int identify() const
Returns a number to identify the sub-classe of Eos the object belongs to.