29 char eos_bf_poly_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Eos/eos_bf_poly.C,v 1.22 2014/10/13 08:52:52 j_novak Exp $" ;
134 #include "eos_bifluid.h" 138 #include "utilitaires.h" 144 double puis(
double,
double) ;
145 double enthal1(
const double x,
const Param& parent) ;
146 double enthal23(
const double x,
const Param& parent) ;
147 double enthal(
const double x,
const Param& parent) ;
159 gam1(2), gam2(2),gam3(1),gam4(1),gam5(1),
160 gam6(1),kap1(kappa1), kap2(kappa2), kap3(kappa3),beta(bet),
161 relax(0.5), precis(1.e-9), ecart(1.e-8)
172 double gamma4,
double gamma5,
double gamma6,
173 double kappa1,
double kappa2,
double kappa3,
174 double bet,
double mass1,
double mass2,
175 double rel,
double prec,
double ec) :
176 Eos_bifluid(
"bi-fluid polytropic EOS", mass1, mass2),
250 cerr <<
"ERROR: could not read all required eos-paramters for Eos_bf_poly()" << endl;
256 if (get_typeos() == 4)
262 else if (get_typeos() == 0)
264 bool slowrot =
false;
265 read_variable (fname, const_cast<char*>(
"slow_rot_style"), slowrot);
273 cerr <<
"ERROR: could not read all required eos-paramters for Eos_bf_poly()" << endl;
328 cout <<
"WARNING!: Eos_bf_poly: the parameters are degenerate!" << endl ;
336 if ((
gam1 ==
double(2)) && (
gam2 ==
double(2)) && (
gam3 ==
double(1))
337 && (
gam4 ==
double(1)) && (
gam5 ==
double(1))
338 && (
gam6 ==
double(0))) {
340 cout <<
"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" <<endl ;
341 cout <<
"WARNING!! The entrainment factor does not depend on" <<endl ;
342 cout <<
" density 2! This may be incorrect and should only be used"<<endl ;
343 cout <<
" for testing purposes..." << endl ;
344 cout <<
"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" <<endl ;
346 else if ((
gam1 ==
double(2)) && (
gam2 ==
double(2)) && (
gam3 ==
double(1))
347 && (
gam4 ==
double(1)) && (
gam5 ==
double(0))
348 && (
gam6 ==
double(1))) {
350 cout <<
"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" <<endl ;
351 cout <<
"WARNING!! The entrainment factor does not depend on" << endl ;
352 cout <<
" density 1! This may be incorrect and should only be used"<<endl ;
353 cout <<
" for testing purposes..." << endl ;
354 cout <<
"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" <<endl ;
356 else if ((
gam1 ==
double(2)) && (
gam2 ==
double(2)) && (
gam3 ==
double(1))
357 && (
gam4 ==
double(1)) && (
gam5 ==
double(1))
358 && (
gam6 ==
double(1))) {
361 else if ((
gam3 ==
double(1)) && (
gam4 ==
double(1)) && (
gam5 ==
double(1))
362 && (
gam6 ==
double(1))) {
365 else if ((
gam3 ==
double(1)) && (
gam5 ==
double(1))) {
368 else if ((
gam4 ==
double(1)) && (
gam6 ==
double(1))) {
376 cout <<
"DEBUG: EOS-type was determined as typeos = " <<
typeos << endl;
391 cout <<
"The second EOS is not of type Eos_bf_poly !" << endl ;
401 <<
"The two Eos_bf_poly have different gammas : " <<
gam1 <<
" <-> " 402 << eos.
gam1 <<
", " <<
gam2 <<
" <-> " 403 << eos.
gam2 <<
", " <<
gam3 <<
" <-> " 404 << eos.
gam3 <<
", " <<
gam4 <<
" <-> " 405 << eos.
gam4 <<
", " <<
gam5 <<
" <-> " 406 << eos.
gam5 <<
", " <<
gam6 <<
" <-> " 407 << eos.
gam6 << endl ;
413 <<
"The two Eos_bf_poly have different kappas : " <<
kap1 <<
" <-> " 414 << eos.
kap1 <<
", " <<
kap2 <<
" <-> " 415 << eos.
kap2 <<
", " <<
kap3 <<
" <-> " 416 << eos.
kap3 << endl ;
422 <<
"The two Eos_bf_poly have different betas : " <<
beta <<
" <-> " 423 << eos.
beta << endl ;
429 <<
"The two Eos_bf_poly have different masses : " <<
m_1 <<
" <-> " 430 << eos.
m_1 <<
", " <<
m_2 <<
" <-> " 474 ost <<
"EOS of class Eos_bf_poly (relativistic polytrope) : " << endl ;
475 ost <<
" Adiabatic index gamma1 : " <<
gam1 << endl ;
476 ost <<
" Adiabatic index gamma2 : " <<
gam2 << endl ;
477 ost <<
" Adiabatic index gamma3 : " <<
gam3 << endl ;
478 ost <<
" Adiabatic index gamma4 : " <<
gam4 << endl ;
479 ost <<
" Adiabatic index gamma5 : " <<
gam5 << endl ;
480 ost <<
" Adiabatic index gamma6 : " <<
gam6 << endl ;
481 ost <<
" Pressure coefficient kappa1 : " <<
kap1 <<
482 " rho_nuc c^2 / n_nuc^gamma" << endl ;
483 ost <<
" Pressure coefficient kappa2 : " <<
kap2 <<
484 " rho_nuc c^2 / n_nuc^gamma" << endl ;
485 ost <<
" Pressure coefficient kappa3 : " <<
kap3 <<
486 " rho_nuc c^2 / n_nuc^gamma" << endl ;
487 ost <<
" Coefficient beta : " <<
beta <<
488 " rho_nuc c^2 / n_nuc^gamma" << endl ;
489 ost <<
" Type for inversion : " <<
typeos << endl ;
490 ost <<
" Parameters for inversion (used if typeos = 4) : " << endl ;
491 ost <<
" relaxation : " <<
relax << endl ;
492 ost <<
" precision for zerosec_b : " <<
precis << endl ;
493 ost <<
" final discrepancy in the result : " <<
ecart << endl ;
508 const double delta2,
double& nbar1,
509 double& nbar2)
const {
514 bool one_fluid =
false;
521 double determ =
kap1*
kap2 - kpd*kpd ;
525 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
529 double b1 =
exp(ent1) -
m_1 ;
530 double b2 =
exp(ent2) -
m_2 ;
532 if (fabs(denom) < 1.e-15) {
535 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
550 double f0 = enthal1(0.,parent) ;
551 if (fabs(f0)<1.e-15) {
555 assert (fabs(cas) > 1.e-15) ;
557 xmax = cas/fabs(cas) ;
560 if (fabs(xmax) > 1.e10) {
561 cout <<
"Problem in Eos_bf_poly::nbar_ent_p!" << endl ;
562 cout << f0 <<
", " << cas << endl ;
563 cout <<
"typeos = 1" << endl ;
566 }
while (enthal1(xmax,parent)*f0 > 0.) ;
567 double l_precis = 1.e-12 ;
570 nbar1 =
zerosec_b(enthal1, parent, xmin, xmax, l_precis,
573 nbar2 = (b1 - coef1*puis(nbar1,
gam1m1))/denom ;
574 double resu1 = coef1*puis(nbar1,
gam1m1) + denom*nbar2 ;
576 double erreur = fabs((
log(fabs(1+resu1))-ent1)/ent1) +
577 fabs((
log(fabs(1+resu2))-ent2)/ent2) ;
578 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)||(erreur > 1.e-4)) ;
583 double b1 =
exp(ent1) -
m_1 ;
584 double b2 =
exp(ent2) -
m_2 ;
586 if (fabs(denom) < 1.e-15) {
589 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
592 double coef0 =
beta*delta2 ;
595 double coef3 = 1./
gam1m1 ;
609 double f0 = enthal23(0.,parent) ;
610 if (fabs(f0)<1.e-15) {
613 double pmax = (fabs(
kap3) < 1.e-15 ? 0. :
gam4*(
gam4-1) ) ;
614 double ptemp = (fabs(
kap3*coef0) < 1.e-15 ? 0.
616 pmax = (pmax>ptemp ? pmax : ptemp) ;
617 ptemp = (fabs(
kap3*coef0) < 1.e-15 ? 0. :
gam4*(
gam6-1) ) ;
618 pmax = (pmax>ptemp ? pmax : ptemp) ;
619 ptemp = (fabs(coef0) < 1.e-15 ? 0. :
gam6*(
gam6-1) ) ;
620 pmax = (pmax>ptemp ? pmax : ptemp) ;
623 assert (fabs(cas) > 1.e-15) ;
625 xmax = cas/fabs(cas) ;
628 if (fabs(xmax) > 1.e10) {
629 cout <<
"Problem in Eos_bf_poly::nbar_ent_p!" << endl ;
630 cout <<
"typeos = 2" << endl ;
633 }
while (enthal23(xmax,parent)*f0 > 0.) ;
634 double l_precis = 1.e-12 ;
637 nbar2 =
zerosec_b(enthal23, parent, xmin, xmax, l_precis,
640 nbar1 = (b1 -
kap3*puis(nbar2,
gam4) - coef0*puis(nbar2,
gam6))
642 nbar1 = puis(nbar1,coef3) ;
644 + coef0*puis(nbar2,
gam6) ;
645 double resu2 = coef2*puis(nbar2,
gam2m1)
647 +
gam6*coef0*puis(nbar2,
gam6-1)*nbar1 ;
648 double erreur = fabs((
log(fabs(1+resu1))-ent1)/ent1) +
649 fabs((
log(fabs(1+resu2))-ent2)/ent2) ;
651 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)||(erreur > 1.e-4)) ;
656 double b1 =
exp(ent1) -
m_1 ;
657 double b2 =
exp(ent2) -
m_2 ;
659 if (fabs(denom) < 1.e-15) {
662 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
665 double coef0 =
beta*delta2 ;
668 double coef3 = 1./
gam2m1 ;
682 double f0 = enthal23(0.,parent) ;
683 if (fabs(f0)<1.e-15) {
686 double pmax = (fabs(
kap3) < 1.e-15 ? 0. :
gam3*(
gam3-1) ) ;
687 double ptemp = (fabs(
kap3*coef0) < 1.e-15 ? 0.
689 pmax = (pmax>ptemp ? pmax : ptemp) ;
690 ptemp = (fabs(
kap3*coef0) < 1.e-15 ? 0. :
gam3*(
gam5-1) ) ;
691 pmax = (pmax>ptemp ? pmax : ptemp) ;
692 ptemp = (fabs(coef0) < 1.e-15 ? 0. :
gam5*(
gam5-1) ) ;
693 pmax = (pmax>ptemp ? pmax : ptemp) ;
696 assert (fabs(cas) > 1.e-15) ;
698 xmax = cas/fabs(cas) ;
701 if (fabs(xmax) > 1.e10) {
702 cout <<
"Problem in Eos_bf_poly::nbar_ent_p!" << endl ;
703 cout <<
"typeos = 3" << endl ;
706 }
while (enthal23(xmax,parent)*f0 > 0.) ;
707 double l_precis = 1.e-12 ;
710 nbar1 =
zerosec_b(enthal23, parent, xmin, xmax, l_precis,
713 nbar2 = (b2 -
kap3*puis(nbar1,
gam3) - coef0*puis(nbar1,
gam5))
715 nbar2 = puis(nbar2,coef3) ;
716 double resu1 = coef1*puis(nbar1,
gam1m1)
718 + coef0*
gam5*puis(nbar1,
gam5-1)*nbar2 ;
719 double resu2 = coef2*puis(nbar2,
gam2m1)
721 double erreur = fabs((
log(fabs(1+resu1))-ent1)/ent1) +
722 fabs((
log(fabs(1+resu2))-ent2)/ent2) ;
723 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)||(erreur > 1.e-4)) ;
728 double b1 =
exp(ent1) -
m_1 ;
729 double b2 =
exp(ent2) -
m_2 ;
731 if (fabs(denom) < 1.e-15) {
734 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
749 double n1l, n2l, n1s, n2s ;
751 double delta = a11*a22 - (a13+a14)*(a23+a24) ;
752 n1l = (a22*b1 - (a13+a14)*b2)/delta ;
753 n2l = (a11*b2 - (a23+a24)*b1)/delta ;
754 n1s = puis(b1/a11, 1./(
gam1m1)) ;
755 n2s = puis(b2/a22, 1./(
gam2m1)) ;
757 double n1m = (n1l<0. ? n1s :
sqrt(n1l*n1s)) ;
758 double n2m = (n2l<0. ? n2s :
sqrt(n2l*n2s)) ;
762 double c1, c2, c3, c4, c5, c6, c7 ;
767 c5 = a13*puis(n2m,
gam4) ;
768 c6 = a14*puis(n2m,
gam6) ;
778 double d1, d2, d3, d4, d5, d6, d7 ;
783 d5 = a23*puis(n1m,
gam3) ;
784 d6 = a24*puis(n1m,
gam5) ;
794 double xmin = -3*(n1s>n2s ? n1s : n2s) ;
795 double xmax = 3*(n1s>n2s ? n1s : n2s) ;
810 sortie = (fabs((n1m-n1)/(n1m+n1)) + fabs((n2m-n2)/(n2m+n2)) >
ecart) ;
832 }
while (sortie&&(mer<nmermax)) ;
855 - kap3*(
exp(ent2) -
m_2)) / determ ;
858 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
865 - kap3*(
exp(ent2) -
m_2 -
beta*delta2)) / determ ;
867 - kap3*(
exp(ent1) -
m_1)) / determ ;
868 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
890 const double delta2)
const {
892 if (( nbar1 >
double(0) ) || ( nbar2 >
double(0))) {
894 double n1 = (nbar1>double(0) ? nbar1 : double(0)) ;
895 double n2 = (nbar2>double(0) ? nbar2 : double(0)) ;
896 double x2 = ((nbar1>double(0))&&(nbar2>
double(0))) ? delta2 : 0 ;
910 const double delta2)
const {
912 if (( nbar1 >
double(0) ) || ( nbar2 >
double(0))) {
914 double n1 = (nbar1>double(0) ? nbar1 : double(0)) ;
915 double n2 = (nbar2>double(0) ? nbar2 : double(0)) ;
916 double x2 = ((nbar1>double(0))&&(nbar2>
double(0))) ? delta2 : 0 ;
964 if ((n1 <= 0.) || (n2 <= 0.)) { xx = 0.; }
966 double gamma_delta3 =
pow(1-delta2,-1.5) ;
989 double puis(
double x,
double p) {
991 if (p==0.)
return (x>=0 ? 1 : -1) ;
993 if (x<0.)
return (-
pow(-x,p)) ;
994 else return pow(x,p) ;
1000 double enthal1(
const double x,
const Param& parent) {
1009 double enthal23(
const double x,
const Param& parent) {
1025 double enthal(
const double x,
const Param& parent) {
1036 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: ].
int get_n_double() const
Returns the number of stored double 's addresses.
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 get_n_double_mod() const
Returns the number of stored modifiable double 's addresses.
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.
const double & get_double(int position=0) const
Returns the reference of a double stored in the list.
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}.