137 #include "eos_bifluid.h" 141 #include "utilitaires.h" 147 double puis(
double,
double) ;
148 double enthal1(
const double x,
const Param& parent) ;
149 double enthal23(
const double x,
const Param& parent) ;
150 double enthal(
const double x,
const Param& parent) ;
162 gam1(2), gam2(2),gam3(1),gam4(1),gam5(1),
163 gam6(1),kap1(kappa1), kap2(kappa2), kap3(kappa3),beta(bet),
164 relax(0.5), precis(1.e-9), ecart(1.e-8)
175 double gamma4,
double gamma5,
double gamma6,
176 double kappa1,
double kappa2,
double kappa3,
177 double bet,
double mass1,
double mass2,
178 double rel,
double prec,
double ec) :
179 Eos_bifluid(
"bi-fluid polytropic EOS", mass1, mass2),
180 gam1(gamma1),gam2(gamma2),gam3(gamma3),gam4(gamma4),gam5(gamma5),
181 gam6(gamma6),kap1(kappa1),kap2(kappa2),kap3(kappa3),beta(bet),
182 relax(rel), precis(prec), ecart(ec)
194 gam1(eosi.gam1), gam2(eosi.gam2), gam3(eosi.gam3),
195 gam4(eosi.gam4), gam5(eosi.gam5), gam6(eosi.gam6),
196 kap1(eosi.kap1), kap2(eosi.kap2), kap3(eosi.kap3),
198 typeos(eosi.typeos), relax(eosi.relax), precis(eosi.precis),
235 Eos_bifluid(fname), relax(0.5), precis(1.e-8), ecart(1.e-7)
253 cerr <<
"ERROR: could not read all required eos-paramters for Eos_bf_poly()" << endl;
259 if (get_typeos() == 4)
265 else if (get_typeos() == 0)
267 bool slowrot =
false;
268 read_variable (fname, const_cast<char*>(
"slow_rot_style"), slowrot);
276 cerr <<
"ERROR: could not read all required eos-paramters for Eos_bf_poly()" << endl;
331 cout <<
"WARNING!: Eos_bf_poly: the parameters are degenerate!" << endl ;
339 if ((
gam1 ==
double(2)) && (
gam2 ==
double(2)) && (
gam3 ==
double(1))
340 && (
gam4 ==
double(1)) && (
gam5 ==
double(1))
341 && (
gam6 ==
double(0))) {
343 cout <<
"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" <<endl ;
344 cout <<
"WARNING!! The entrainment factor does not depend on" <<endl ;
345 cout <<
" density 2! This may be incorrect and should only be used"<<endl ;
346 cout <<
" for testing purposes..." << endl ;
347 cout <<
"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" <<endl ;
349 else if ((
gam1 ==
double(2)) && (
gam2 ==
double(2)) && (
gam3 ==
double(1))
350 && (
gam4 ==
double(1)) && (
gam5 ==
double(0))
351 && (
gam6 ==
double(1))) {
353 cout <<
"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" <<endl ;
354 cout <<
"WARNING!! The entrainment factor does not depend on" << endl ;
355 cout <<
" density 1! This may be incorrect and should only be used"<<endl ;
356 cout <<
" for testing purposes..." << endl ;
357 cout <<
"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" <<endl ;
359 else if ((
gam1 ==
double(2)) && (
gam2 ==
double(2)) && (
gam3 ==
double(1))
360 && (
gam4 ==
double(1)) && (
gam5 ==
double(1))
361 && (
gam6 ==
double(1))) {
364 else if ((
gam3 ==
double(1)) && (
gam4 ==
double(1)) && (
gam5 ==
double(1))
365 && (
gam6 ==
double(1))) {
368 else if ((
gam3 ==
double(1)) && (
gam5 ==
double(1))) {
371 else if ((
gam4 ==
double(1)) && (
gam6 ==
double(1))) {
379 cout <<
"DEBUG: EOS-type was determined as typeos = " <<
typeos << endl;
394 cout <<
"The second EOS is not of type Eos_bf_poly !" << endl ;
404 <<
"The two Eos_bf_poly have different gammas : " <<
gam1 <<
" <-> " 405 << eos.
gam1 <<
", " <<
gam2 <<
" <-> " 406 << eos.
gam2 <<
", " <<
gam3 <<
" <-> " 407 << eos.
gam3 <<
", " <<
gam4 <<
" <-> " 408 << eos.
gam4 <<
", " <<
gam5 <<
" <-> " 409 << eos.
gam5 <<
", " <<
gam6 <<
" <-> " 410 << eos.
gam6 << endl ;
416 <<
"The two Eos_bf_poly have different kappas : " <<
kap1 <<
" <-> " 417 << eos.
kap1 <<
", " <<
kap2 <<
" <-> " 418 << eos.
kap2 <<
", " <<
kap3 <<
" <-> " 419 << eos.
kap3 << endl ;
425 <<
"The two Eos_bf_poly have different betas : " <<
beta <<
" <-> " 426 << eos.
beta << endl ;
432 <<
"The two Eos_bf_poly have different masses : " <<
m_1 <<
" <-> " 433 << eos.
m_1 <<
", " <<
m_2 <<
" <-> " 477 ost <<
"EOS of class Eos_bf_poly (relativistic polytrope) : " << endl ;
478 ost <<
" Adiabatic index gamma1 : " <<
gam1 << endl ;
479 ost <<
" Adiabatic index gamma2 : " <<
gam2 << endl ;
480 ost <<
" Adiabatic index gamma3 : " <<
gam3 << endl ;
481 ost <<
" Adiabatic index gamma4 : " <<
gam4 << endl ;
482 ost <<
" Adiabatic index gamma5 : " <<
gam5 << endl ;
483 ost <<
" Adiabatic index gamma6 : " <<
gam6 << endl ;
484 ost <<
" Pressure coefficient kappa1 : " <<
kap1 <<
485 " rho_nuc c^2 / n_nuc^gamma" << endl ;
486 ost <<
" Pressure coefficient kappa2 : " <<
kap2 <<
487 " rho_nuc c^2 / n_nuc^gamma" << endl ;
488 ost <<
" Pressure coefficient kappa3 : " <<
kap3 <<
489 " rho_nuc c^2 / n_nuc^gamma" << endl ;
490 ost <<
" Coefficient beta : " <<
beta <<
491 " rho_nuc c^2 / n_nuc^gamma" << endl ;
492 ost <<
" Type for inversion : " <<
typeos << endl ;
493 ost <<
" Parameters for inversion (used if typeos = 4) : " << endl ;
494 ost <<
" relaxation : " <<
relax << endl ;
495 ost <<
" precision for zerosec_b : " <<
precis << endl ;
496 ost <<
" final discrepancy in the result : " <<
ecart << endl ;
511 const double delta2,
double& nbar1,
512 double& nbar2)
const {
517 bool one_fluid =
false;
524 double determ =
kap1*
kap2 - kpd*kpd ;
528 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
532 double b1 =
exp(ent1) -
m_1 ;
533 double b2 =
exp(ent2) -
m_2 ;
535 if (fabs(denom) < 1.e-15) {
538 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
553 double f0 = enthal1(0.,parent) ;
554 if (fabs(f0)<1.e-15) {
558 assert (fabs(cas) > 1.e-15) ;
560 xmax = cas/fabs(cas) ;
563 if (fabs(xmax) > 1.e10) {
564 cout <<
"Problem in Eos_bf_poly::nbar_ent_p!" << endl ;
565 cout << f0 <<
", " << cas << endl ;
566 cout <<
"typeos = 1" << endl ;
569 }
while (enthal1(xmax,parent)*f0 > 0.) ;
570 double l_precis = 1.e-12 ;
573 nbar1 =
zerosec_b(enthal1, parent, xmin, xmax, l_precis,
576 nbar2 = (b1 - coef1*puis(nbar1,
gam1m1))/denom ;
577 double resu1 = coef1*puis(nbar1,
gam1m1) + denom*nbar2 ;
579 double erreur = fabs((
log(fabs(1+resu1))-ent1)/ent1) +
580 fabs((
log(fabs(1+resu2))-ent2)/ent2) ;
581 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)||(erreur > 1.e-4)) ;
586 double b1 =
exp(ent1) -
m_1 ;
587 double b2 =
exp(ent2) -
m_2 ;
589 if (fabs(denom) < 1.e-15) {
592 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
595 double coef0 =
beta*delta2 ;
598 double coef3 = 1./
gam1m1 ;
612 double f0 = enthal23(0.,parent) ;
613 if (fabs(f0)<1.e-15) {
616 double pmax = (fabs(
kap3) < 1.e-15 ? 0. :
gam4*(
gam4-1) ) ;
617 double ptemp = (fabs(
kap3*coef0) < 1.e-15 ? 0.
619 pmax = (pmax>ptemp ? pmax : ptemp) ;
620 ptemp = (fabs(
kap3*coef0) < 1.e-15 ? 0. :
gam4*(
gam6-1) ) ;
621 pmax = (pmax>ptemp ? pmax : ptemp) ;
622 ptemp = (fabs(coef0) < 1.e-15 ? 0. :
gam6*(
gam6-1) ) ;
623 pmax = (pmax>ptemp ? pmax : ptemp) ;
626 assert (fabs(cas) > 1.e-15) ;
628 xmax = cas/fabs(cas) ;
631 if (fabs(xmax) > 1.e10) {
632 cout <<
"Problem in Eos_bf_poly::nbar_ent_p!" << endl ;
633 cout <<
"typeos = 2" << endl ;
636 }
while (enthal23(xmax,parent)*f0 > 0.) ;
637 double l_precis = 1.e-12 ;
640 nbar2 =
zerosec_b(enthal23, parent, xmin, xmax, l_precis,
643 nbar1 = (b1 -
kap3*puis(nbar2,
gam4) - coef0*puis(nbar2,
gam6))
645 nbar1 = puis(nbar1,coef3) ;
647 + coef0*puis(nbar2,
gam6) ;
648 double resu2 = coef2*puis(nbar2,
gam2m1)
650 +
gam6*coef0*puis(nbar2,
gam6-1)*nbar1 ;
651 double erreur = fabs((
log(fabs(1+resu1))-ent1)/ent1) +
652 fabs((
log(fabs(1+resu2))-ent2)/ent2) ;
654 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)||(erreur > 1.e-4)) ;
659 double b1 =
exp(ent1) -
m_1 ;
660 double b2 =
exp(ent2) -
m_2 ;
662 if (fabs(denom) < 1.e-15) {
665 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
668 double coef0 =
beta*delta2 ;
671 double coef3 = 1./
gam2m1 ;
685 double f0 = enthal23(0.,parent) ;
686 if (fabs(f0)<1.e-15) {
689 double pmax = (fabs(
kap3) < 1.e-15 ? 0. :
gam3*(
gam3-1) ) ;
690 double ptemp = (fabs(
kap3*coef0) < 1.e-15 ? 0.
692 pmax = (pmax>ptemp ? pmax : ptemp) ;
693 ptemp = (fabs(
kap3*coef0) < 1.e-15 ? 0. :
gam3*(
gam5-1) ) ;
694 pmax = (pmax>ptemp ? pmax : ptemp) ;
695 ptemp = (fabs(coef0) < 1.e-15 ? 0. :
gam5*(
gam5-1) ) ;
696 pmax = (pmax>ptemp ? pmax : ptemp) ;
699 assert (fabs(cas) > 1.e-15) ;
701 xmax = cas/fabs(cas) ;
704 if (fabs(xmax) > 1.e10) {
705 cout <<
"Problem in Eos_bf_poly::nbar_ent_p!" << endl ;
706 cout <<
"typeos = 3" << endl ;
709 }
while (enthal23(xmax,parent)*f0 > 0.) ;
710 double l_precis = 1.e-12 ;
713 nbar1 =
zerosec_b(enthal23, parent, xmin, xmax, l_precis,
716 nbar2 = (b2 -
kap3*puis(nbar1,
gam3) - coef0*puis(nbar1,
gam5))
718 nbar2 = puis(nbar2,coef3) ;
719 double resu1 = coef1*puis(nbar1,
gam1m1)
721 + coef0*
gam5*puis(nbar1,
gam5-1)*nbar2 ;
722 double resu2 = coef2*puis(nbar2,
gam2m1)
724 double erreur = fabs((
log(fabs(1+resu1))-ent1)/ent1) +
725 fabs((
log(fabs(1+resu2))-ent2)/ent2) ;
726 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)||(erreur > 1.e-4)) ;
731 double b1 =
exp(ent1) -
m_1 ;
732 double b2 =
exp(ent2) -
m_2 ;
734 if (fabs(denom) < 1.e-15) {
737 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
752 double n1l, n2l, n1s, n2s ;
754 double delta = a11*a22 - (a13+a14)*(a23+a24) ;
755 n1l = (a22*b1 - (a13+a14)*b2)/delta ;
756 n2l = (a11*b2 - (a23+a24)*b1)/delta ;
757 n1s = puis(b1/a11, 1./(
gam1m1)) ;
758 n2s = puis(b2/a22, 1./(
gam2m1)) ;
760 double n1m = (n1l<0. ? n1s :
sqrt(n1l*n1s)) ;
761 double n2m = (n2l<0. ? n2s :
sqrt(n2l*n2s)) ;
765 double c1, c2, c3, c4, c5, c6, c7 ;
770 c5 = a13*puis(n2m,
gam4) ;
771 c6 = a14*puis(n2m,
gam6) ;
781 double d1, d2, d3, d4, d5, d6, d7 ;
786 d5 = a23*puis(n1m,
gam3) ;
787 d6 = a24*puis(n1m,
gam5) ;
797 double xmin = -3*(n1s>n2s ? n1s : n2s) ;
798 double xmax = 3*(n1s>n2s ? n1s : n2s) ;
813 sortie = (fabs((n1m-n1)/(n1m+n1)) + fabs((n2m-n2)/(n2m+n2)) >
ecart) ;
835 }
while (sortie&&(mer<nmermax)) ;
861 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
871 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
893 const double delta2)
const {
895 if (( nbar1 >
double(0) ) || ( nbar2 >
double(0))) {
897 double n1 = (nbar1>double(0) ? nbar1 : double(0)) ;
898 double n2 = (nbar2>double(0) ? nbar2 : double(0)) ;
899 double x2 = ((nbar1>double(0))&&(nbar2>
double(0))) ? delta2 : 0 ;
913 const double delta2)
const {
915 if (( nbar1 >
double(0) ) || ( nbar2 >
double(0))) {
917 double n1 = (nbar1>double(0) ? nbar1 : double(0)) ;
918 double n2 = (nbar2>double(0) ? nbar2 : double(0)) ;
919 double x2 = ((nbar1>double(0))&&(nbar2>
double(0))) ? delta2 : 0 ;
967 if ((n1 <= 0.) || (n2 <= 0.)) { xx = 0.; }
969 double gamma_delta3 =
pow(1-delta2,-1.5) ;
992 double puis(
double x,
double p) {
994 if (p==0.)
return (x>=0 ? 1 : -1) ;
996 if (x<0.)
return (-
pow(-x,p)) ;
997 else return pow(x,p) ;
1003 double enthal1(
const double x,
const Param& parent) {
1004 assert(parent.get_n_double() == 7) ;
1006 return parent.get_double(0)*puis(parent.get_double(1) - parent.get_double(2)
1007 *puis(x,parent.get_double(3)), parent.get_double(4))
1008 + parent.get_double(5)*x - parent.get_double(6) ;
1012 double enthal23(
const double x,
const Param& parent) {
1013 assert(parent.get_n_double() == 10) ;
1015 double nx = (parent.get_double(0) - parent.get_double(1)*
1016 puis(x,parent.get_double(2)) - parent.get_double(3)*
1017 puis(x,parent.get_double(4)) )/parent.get_double(5) ;
1018 nx = puis(nx,parent.get_double(6)) ;
1019 return parent.get_double(7)*puis(x,parent.get_double(8))
1020 + parent.get_double(1)*parent.get_double(2)*nx*
1021 puis(x,parent.get_double(2) - 1)
1022 + parent.get_double(3)*parent.get_double(4)*nx*
1023 puis(x,parent.get_double(4) - 1)
1024 - parent.get_double(9) ;
1028 double enthal(
const double x,
const Param& parent) {
1029 assert(parent.get_n_double_mod() == 7) ;
1031 double alp1 = parent.get_double_mod(0) ;
1032 double alp2 = parent.get_double_mod(1) ;
1033 double alp3 = parent.get_double_mod(2) ;
1034 double cc1 = parent.get_double_mod(3) ;
1035 double cc2 = parent.get_double_mod(4) ;
1036 double cc3 = parent.get_double_mod(5) ;
1037 double cc4 = parent.get_double_mod(6) ;
1039 return (cc1*puis(x,alp1) + cc2*puis(x,alp2) + cc3*puis(x,alp3) - cc4) ;
virtual double get_K12(const double n1, const double n2, const double delta2) const
Computes the derivative of the energy with respect to .
virtual void sauve(FILE *) const
Save in a file.
double m_1
Individual particle mass [unit: ].
Cmp log(const Cmp &)
Neperian logarithm.
double & get_double_mod(int position=0) const
Returns the reference of a stored modifiable double .
Cmp exp(const Cmp &)
Exponential.
double kap2
Pressure coefficient , see Eq.
Cmp sqrt(const Cmp &)
Square root.
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.
double relax
Parameters needed for some inversions of the EOS.
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.
Analytic equation of state for two fluids (relativistic case).
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 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.
void determine_type()
Determines the type of the analytical EOS (see typeos )
virtual Eos * trans2Eos() const
Makes a translation from Eos_bifluid to Eos .
virtual ~Eos_bf_poly()
Destructor.
void operator=(const Eos_bifluid &)
Assignment to another Eos_bifluid.
virtual ostream & operator>>(ostream &) const
Operator >>
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: ].
void set_auxiliary()
Computes the auxiliary quantities gam1m1 , gam2m1 and gam3m1.
virtual double get_K22(const double n1, const double n2, const double delta2) const
Computes the derivative of the energy/(baryonic density 2) .
virtual int identify() const
Returns a number to identify the sub-classe of Eos_bifluid the object belongs to.
double kap3
Pressure coefficient , see Eq.
Polytropic equation of state (relativistic case).
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.
Cmp pow(const Cmp &, int)
Power .
virtual double nbar_ent_p2(const double ent2) const
Computes baryon density out of the log-enthalpy assuming that only fluid 2 is present.
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 void sauve(FILE *) const
Save in a file.
double precis
contains the precision required in zerosec_b
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 bool operator!=(const Eos_bifluid &) const
Comparison operator (difference)
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) .
int read_variable(const char *fname, const char *var_name, char *fmt, void *varp)
Reads a variable from file.
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.
double kap1
Pressure coefficient , see Eq.
virtual double nbar_ent_p1(const double ent1) const
Computes baryon density out of the log-enthalpy asuming that only fluid 1 is present.
double gam1
Adiabatic indexes , see Eq.~eeosbfpolye}.
double gam6
Adiabatic indexes , see Eq.~eeosbfpolye}.
virtual bool operator==(const Eos_bifluid &) const
Comparison operator (egality)
double gam3
Adiabatic indexes , see Eq.~eeosbfpolye}.
Eos_bf_poly(double kappa1, double kappa2, double kappa3, double beta)
Standard constructor.
double gam4
Adiabatic indexes , see Eq.~eeosbfpolye}.
double gam2
Adiabatic indexes , see Eq.~eeosbfpolye}.