98 #include "eos_bifluid.h" 101 #include "utilitaires.h" 107 double puis(
double,
double) ;
108 double enthal1(
const double x,
const Param& parent) ;
109 double enthal23(
const double x,
const Param& parent) ;
110 double enthal(
const double x,
const Param& parent) ;
123 name =
"bi-fluid polytropic non-relativistic EOS" ;
129 double gamma4,
double gamma5,
double gamma6,
130 double kappa1,
double kappa2,
double kappa3,
131 double bet,
double mass1,
double mass2,
132 double l_relax,
double l_precis,
double l_ecart) :
133 Eos_bf_poly(gamma1, gamma2, gamma3, gamma4, gamma5, gamma6,
134 kappa1, kappa2, kappa3, bet, mass1, mass2, l_relax, l_precis,
136 name =
"bi-fluid polytropic non-relativistic EOS" ;
184 cout <<
"The second EOS is not of type Eos_bf_poly_newt !" << endl ;
195 <<
"The two Eos_bf_poly_newt have different gammas : " <<
gam1 <<
" <-> " 196 << eos.
gam1 <<
", " <<
gam2 <<
" <-> " 197 << eos.
gam2 <<
", " <<
gam3 <<
" <-> " 198 << eos.
gam3 <<
", " <<
gam4 <<
" <-> " 199 << eos.
gam4 <<
", " <<
gam5 <<
" <-> " 200 << eos.
gam5 <<
", " <<
gam6 <<
" <-> " 201 << eos.
gam6 << endl ;
207 <<
"The two Eos_bf_poly_newt have different kappas : " <<
kap1 <<
" <-> " 208 << eos.
kap1 <<
", " <<
kap2 <<
" <-> " 209 << eos.
kap2 <<
", " <<
kap3 <<
" <-> " 210 << eos.
kap3 << endl ;
216 <<
"The two Eos_bf_poly_newt have different betas : " <<
beta <<
" <-> " 217 << eos.
beta << endl ;
223 <<
"The two Eos_bf_poly_newt have different masses : " <<
m_1 <<
" <-> " 224 << eos.
m_1 <<
", " <<
m_2 <<
" <-> " 253 ost <<
"EOS of class Eos_bf_poly_newt (non-relativistic polytrope) : " << endl ;
254 ost <<
" Adiabatic index gamma1 : " <<
gam1 << endl ;
255 ost <<
" Adiabatic index gamma2 : " <<
gam2 << endl ;
256 ost <<
" Adiabatic index gamma3 : " <<
gam3 << endl ;
257 ost <<
" Adiabatic index gamma4 : " <<
gam4 << endl ;
258 ost <<
" Adiabatic index gamma5 : " <<
gam5 << endl ;
259 ost <<
" Adiabatic index gamma6 : " <<
gam6 << endl ;
260 ost <<
" Pressure coefficient kappa1 : " <<
kap1 <<
261 " rho_nuc c^2 / n_nuc^gamma" << endl ;
262 ost <<
" Pressure coefficient kappa2 : " <<
kap2 <<
263 " rho_nuc c^2 / n_nuc^gamma" << endl ;
264 ost <<
" Pressure coefficient kappa3 : " <<
kap3 <<
265 " rho_nuc c^2 / n_nuc^gamma" << endl ;
266 ost <<
" Coefficient beta : " <<
beta <<
267 " rho_nuc c^2 / n_nuc^gamma" << endl ;
268 ost <<
" Mean particle 1 mass : " <<
m_1 <<
" m_B" << endl ;
269 ost <<
" Mean particle 2 mass : " <<
m_2 <<
" m_B" << endl ;
270 ost <<
" Parameters for inversion (used if typeos = 4 : " << endl ;
271 ost <<
" relaxation : " <<
relax << endl ;
272 ost <<
" precision for zerosec_b : " <<
precis << endl ;
273 ost <<
" final discrepancy in the result : " <<
ecart << endl ;
289 const double delta2,
double& nbar1,
290 double& nbar2)
const {
295 bool one_fluid =
false;
302 double determ =
kap1*
kap2 - kpd*kpd ;
304 nbar1 = (
kap2*ent1*
m_1 - kpd*ent2*
m_2) / determ ;
305 nbar2 = (
kap1*ent2*
m_2 - kpd*ent1*
m_1) / determ ;
306 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
310 double b1 = ent1*
m_1 ;
311 double b2 = ent2*
m_2 ;
313 if (fabs(denom) < 1.e-15) {
316 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
331 double f0 = enthal1(0.,parent) ;
332 if (fabs(f0)<1.e-15) {
336 assert (fabs(cas) > 1.e-15) ;
338 xmax = cas/fabs(cas) ;
341 if (fabs(xmax) > 1.e10) {
342 cout <<
"Problem in Eos_bf_poly::nbar_ent_p!" << endl ;
343 cout << f0 <<
", " << cas << endl ;
344 cout <<
"typeos = 1" << endl ;
347 }
while (enthal1(xmax,parent)*f0 > 0.) ;
348 double l_precis = 1.e-12 ;
351 nbar1 =
zerosec_b(enthal1, parent, xmin, xmax, l_precis,
354 nbar2 = (b1 - coef1*puis(nbar1,
gam1m1))/denom ;
355 double resu1 = coef1*puis(nbar1,
gam1m1) + denom*nbar2 ;
357 double erreur = fabs((resu1/
m_1-ent1)/ent1) +
358 fabs((resu2/
m_2-ent2)/ent2) ;
359 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)||(erreur > 1.e-4)) ;
364 double b1 = ent1*
m_1 ;
365 double b2 = ent2*
m_2 ;
367 if (fabs(denom) < 1.e-15) {
370 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
373 double coef0 =
beta*delta2 ;
376 double coef3 = 1./
gam1m1 ;
390 double f0 = enthal23(0.,parent) ;
391 if (fabs(f0)<1.e-15) {
394 double pmax = (fabs(
kap3) < 1.e-15 ? 0. :
gam4*(
gam4-1) ) ;
395 double ptemp = (fabs(
kap3*coef0) < 1.e-15 ? 0.
397 pmax = (pmax>ptemp ? pmax : ptemp) ;
398 ptemp = (fabs(
kap3*coef0) < 1.e-15 ? 0. :
gam4*(
gam6-1) ) ;
399 pmax = (pmax>ptemp ? pmax : ptemp) ;
400 ptemp = (fabs(coef0) < 1.e-15 ? 0. :
gam6*(
gam6-1) ) ;
401 pmax = (pmax>ptemp ? pmax : ptemp) ;
404 assert (fabs(cas) > 1.e-15) ;
406 xmax = cas/fabs(cas) ;
409 if (fabs(xmax) > 1.e10) {
410 cout <<
"Problem in Eos_bf_poly::nbar_ent_p!" << endl ;
411 cout <<
"typeos = 2" << endl ;
414 }
while (enthal23(xmax,parent)*f0 > 0.) ;
415 double l_precis = 1.e-12 ;
418 nbar2 =
zerosec_b(enthal23, parent, xmin, xmax, l_precis,
421 nbar1 = (b1 -
kap3*puis(nbar2,
gam4) - coef0*puis(nbar2,
gam6))
423 nbar1 = puis(nbar1,coef3) ;
425 + coef0*puis(nbar2,
gam6) ;
426 double resu2 = coef2*puis(nbar2,
gam2m1)
428 +
gam6*coef0*puis(nbar2,
gam6-1)*nbar1 ;
429 double erreur = fabs((resu1/
m_1-ent1)/ent1) +
430 fabs((resu2/
m_2-ent2)/ent2) ;
432 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)||(erreur > 1.e-4)) ;
437 double b1 = ent1*
m_1 ;
438 double b2 = ent2*
m_2 ;
440 if (fabs(denom) < 1.e-15) {
443 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
446 double coef0 =
beta*delta2 ;
449 double coef3 = 1./
gam2m1 ;
463 double f0 = enthal23(0.,parent) ;
464 if (fabs(f0)<1.e-15) {
467 double pmax = (fabs(
kap3) < 1.e-15 ? 0. :
gam3*(
gam3-1) ) ;
468 double ptemp = (fabs(
kap3*coef0) < 1.e-15 ? 0.
470 pmax = (pmax>ptemp ? pmax : ptemp) ;
471 ptemp = (fabs(
kap3*coef0) < 1.e-15 ? 0. :
gam3*(
gam5-1) ) ;
472 pmax = (pmax>ptemp ? pmax : ptemp) ;
473 ptemp = (fabs(coef0) < 1.e-15 ? 0. :
gam5*(
gam5-1) ) ;
474 pmax = (pmax>ptemp ? pmax : ptemp) ;
477 assert (fabs(cas) > 1.e-15) ;
479 xmax = cas/fabs(cas) ;
482 if (fabs(xmax) > 1.e10) {
483 cout <<
"Problem in Eos_bf_poly::nbar_ent_p!" << endl ;
484 cout <<
"typeos = 3" << endl ;
487 }
while (enthal23(xmax,parent)*f0 > 0.) ;
488 double l_precis = 1.e-12 ;
491 nbar1 =
zerosec_b(enthal23, parent, xmin, xmax, l_precis,
494 nbar2 = (b2 -
kap3*puis(nbar1,
gam3) - coef0*puis(nbar1,
gam5))
496 nbar2 = puis(nbar2,coef3) ;
497 double resu1 = coef1*puis(nbar1,
gam1m1)
499 + coef0*
gam5*puis(nbar1,
gam5-1)*nbar2 ;
500 double resu2 = coef2*puis(nbar2,
gam2m1)
502 double erreur = fabs((resu1/
m_1-ent1)/ent1) +
503 fabs((resu2/
m_2-ent2)/ent2) ;
504 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)||(erreur > 1.e-4)) ;
509 double b1 = ent1*
m_1 ;
510 double b2 = ent2*
m_2 ;
512 if (fabs(denom) < 1.e-15) {
515 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
530 double n1l, n2l, n1s, n2s ;
532 double delta = a11*a22 - (a13+a14)*(a23+a24) ;
533 n1l = (a22*b1 - (a13+a14)*b2)/delta ;
534 n2l = (a11*b2 - (a23+a24)*b1)/delta ;
535 n1s = puis(b1/a11, 1./(
gam1m1)) ;
536 n2s = puis(b2/a22, 1./(
gam2m1)) ;
538 double n1m = (n1l<0. ? n1s :
sqrt(n1l*n1s)) ;
539 double n2m = (n2l<0. ? n2s :
sqrt(n2l*n2s)) ;
543 double c1, c2, c3, c4, c5, c6, c7 ;
548 c5 = a13*puis(n2m,
gam4) ;
549 c6 = a14*puis(n2m,
gam6) ;
559 double d1, d2, d3, d4, d5, d6, d7 ;
564 d5 = a23*puis(n1m,
gam3) ;
565 d6 = a24*puis(n1m,
gam5) ;
575 double xmin = -3*(n1s>n2s ? n1s : n2s) ;
576 double xmax = 3*(n1s>n2s ? n1s : n2s) ;
591 sortie = (fabs((n1m-n1)/(n1m+n1)) + fabs((n2m-n2)/(n2m+n2)) >
ecart) ;
613 }
while (sortie&&(mer<nmermax)) ;
637 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
645 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
668 const double delta2)
const {
670 double n1 = (nbar1>double(0) ? nbar1 : 0) ;
671 double n2 = (nbar2>double(0) ? nbar2 : 0) ;
672 double x2 = ((nbar1>double(0))&&(nbar2>
double(0))) ? delta2 : 0 ;
686 const double delta2)
const {
688 double n1 = (nbar1>double(0) ? nbar1 : 0) ;
689 double n2 = (nbar2>double(0) ? nbar2 : 0) ;
690 double x2 = ((nbar1>double(0))&&(nbar2>
double(0))) ? delta2 : 0 ;
707 if (n1 <= 0.) xx = 0. ;
717 if (n2 <= 0.) xx = 0. ;
727 if ((n1 <= 0.) || (n2 <= 0.)) xx = 0.;
virtual bool operator!=(const Eos_bifluid &) const
Comparison operator (difference)
double m_1
Individual particle mass [unit: ].
double & get_double_mod(int position=0) const
Returns the reference of a stored modifiable double .
double kap2
Pressure coefficient , see Eq.
Cmp sqrt(const Cmp &)
Square root.
double relax
Parameters needed for some inversions of the EOS.
void operator=(const Eos_bf_poly_newt &)
Assignment to another Eos_bf_poly_newt.
virtual bool operator==(const Eos_bifluid &) const
Comparison operator (egality)
Equation of state base class.
virtual int identify() const =0
Returns a number to identify the sub-classe of Eos_bifluid the object belongs to.
void operator=(const Eos_bf_poly &)
Assignment to another Eos_bf_poly.
virtual double get_K12(const double n1, const double n2, const double delta2) const
Computes the derivative of the energy with respect to .
Analytic equation of state for two fluids (relativistic case).
double zerosec_b(double(*f)(double, const Param &), const Param &par, double a, double b, double precis, int nitermax, int &niter)
Finding the zero a function on a bounded domain.
Eos_bf_poly_newt(double kappa1, double kappa2, double kappa3, double beta)
Standard constructor.
virtual double get_K22(const double n1, const double n2, const double delta2) const
Computes the derivative of the energy/(baryonic density 2) .
double ecart
contains the precision required in the relaxation nbar_ent_p
void add_double_mod(double &x, int position=0)
Adds the address of a new modifiable double to the list.
2-fluids equation of state base class.
int typeos
The bi-fluid analytical EOS type:
double beta
Coefficient , see Eq.
double m_2
Individual particle mass [unit: ].
virtual ~Eos_bf_poly_newt()
Destructor.
virtual ostream & operator>>(ostream &) const
Operator >>
virtual double get_K11(const double n1, const double n2, const double delta2) const
Computes the derivative of the energy with respect to (baryonic density 1) .
virtual bool nbar_ent_p(const double ent1, const double ent2, const double delta2, double &nbar1, double &nbar2) const
Computes both baryon densities from the log-enthalpies.
double kap3
Pressure coefficient , see Eq.
Analytic equation of state for two fluids (Newtonian case).
virtual Eos * trans2Eos() const
Makes a translation from Eos_bifluid to Eos .
Cmp pow(const Cmp &, int)
Power .
virtual void sauve(FILE *) const
Save in a file.
Polytropic equation of state (Newtonian case).
virtual double press_nbar_p(const double nbar1, const double nbar2, const double delta2) const
Computes the pressure from the baryonic densities and the relative velocity.
virtual double nbar_ent_p1(const double ent1) const
Computes baryon density out of the log-enthalpy asuming that only fluid 1 is present (virtual functio...
virtual void sauve(FILE *) const
Save in a file.
double precis
contains the precision required in zerosec_b
double gam5
Adiabatic indexes , see Eq.~eeosbfpolye}.
void add_double(const double &x, int position=0)
Adds the the address of a new double to the list.
virtual int identify() const
Returns a number to identify the sub-classe of Eos_bifluid the object belongs to.
double kap1
Pressure coefficient , see Eq.
virtual double ener_nbar_p(const double nbar1, const double nbar2, const double delta2) const
Computes the total energy density from the baryonic densities and the relative velocity.
virtual double nbar_ent_p2(const double ent2) const
Computes baryon density out of the log-enthalpy assuming that only fluid 2 is present.
double gam1
Adiabatic indexes , see Eq.~eeosbfpolye}.
double gam6
Adiabatic indexes , see Eq.~eeosbfpolye}.
double gam3
Adiabatic indexes , see Eq.~eeosbfpolye}.
double gam4
Adiabatic indexes , see Eq.~eeosbfpolye}.
double gam2
Adiabatic indexes , see Eq.~eeosbfpolye}.