29char 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"
144double puis(
double,
double) ;
145double enthal1(
const double x,
const Param& parent) ;
146double enthal23(
const double x,
const Param& parent) ;
147double 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),
177 gam1(gamma1),gam2(gamma2),gam3(gamma3),gam4(gamma4),gam5(gamma5),
178 gam6(gamma6),kap1(kappa1),kap2(kappa2),kap3(kappa3),beta(bet),
179 relax(rel), precis(prec), ecart(ec)
191 gam1(eosi.gam1), gam2(eosi.gam2), gam3(eosi.gam3),
192 gam4(eosi.gam4), gam5(eosi.gam5), gam6(eosi.gam6),
193 kap1(eosi.kap1), kap2(eosi.kap2), kap3(eosi.kap3),
195 typeos(eosi.typeos), relax(eosi.relax), precis(eosi.precis),
232 Eos_bifluid(fname), relax(0.5), precis(1.e-8), ecart(1.e-7)
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)) ;
858 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
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) ;
1000double enthal1(
const double x,
const Param& parent) {
1001 assert(parent.get_n_double() == 7) ;
1003 return parent.get_double(0)*puis(parent.get_double(1) - parent.get_double(2)
1004 *puis(x,parent.get_double(3)), parent.get_double(4))
1005 + parent.get_double(5)*x - parent.get_double(6) ;
1009double enthal23(
const double x,
const Param& parent) {
1010 assert(parent.get_n_double() == 10) ;
1012 double nx = (parent.get_double(0) - parent.get_double(1)*
1013 puis(x,parent.get_double(2)) - parent.get_double(3)*
1014 puis(x,parent.get_double(4)) )/parent.get_double(5) ;
1015 nx = puis(nx,parent.get_double(6)) ;
1016 return parent.get_double(7)*puis(x,parent.get_double(8))
1017 + parent.get_double(1)*parent.get_double(2)*nx*
1018 puis(x,parent.get_double(2) - 1)
1019 + parent.get_double(3)*parent.get_double(4)*nx*
1020 puis(x,parent.get_double(4) - 1)
1021 - parent.get_double(9) ;
1025double enthal(
const double x,
const Param& parent) {
1026 assert(parent.get_n_double_mod() == 7) ;
1028 double alp1 = parent.get_double_mod(0) ;
1029 double alp2 = parent.get_double_mod(1) ;
1030 double alp3 = parent.get_double_mod(2) ;
1031 double cc1 = parent.get_double_mod(3) ;
1032 double cc2 = parent.get_double_mod(4) ;
1033 double cc3 = parent.get_double_mod(5) ;
1034 double cc4 = parent.get_double_mod(6) ;
1036 return (cc1*puis(x,alp1) + cc2*puis(x,alp2) + cc3*puis(x,alp3) - cc4) ;
Analytic equation of state for two fluids (relativistic case).
double gam5
Adiabatic indexes , see Eq.~eeosbfpolye}.
virtual bool operator!=(const Eos_bifluid &) const
Comparison operator (difference)
double kap1
Pressure coefficient , see Eq.
void determine_type()
Determines the type of the analytical EOS (see typeos )
double gam1
Adiabatic indexes , see Eq.~eeosbfpolye}.
double kap2
Pressure coefficient , see Eq.
void operator=(const Eos_bf_poly &)
Assignment to another Eos_bf_poly.
double gam4
Adiabatic indexes , see Eq.~eeosbfpolye}.
double gam3
Adiabatic indexes , see Eq.~eeosbfpolye}.
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 double get_K12(const double n1, const double n2, const double delta2) const
Computes the derivative of the energy with respect to .
virtual bool operator==(const Eos_bifluid &) const
Comparison operator (egality)
virtual ostream & operator>>(ostream &) const
Operator >>
virtual Eos * trans2Eos() const
Makes a translation from Eos_bifluid to Eos .
virtual double get_K22(const double n1, const double n2, const double delta2) const
Computes the derivative of the energy/(baryonic density 2) .
double gam6
Adiabatic indexes , see Eq.~eeosbfpolye}.
double gam2
Adiabatic indexes , see Eq.~eeosbfpolye}.
double precis
contains the precision required in zerosec_b
virtual int identify() const
Returns a number to identify the sub-classe of Eos_bifluid the object belongs to.
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 void sauve(FILE *) const
Save in a file.
void set_auxiliary()
Computes the auxiliary quantities gam1m1 , gam2m1 and gam3m1.
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 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.
int typeos
The bi-fluid analytical EOS type:
double beta
Coefficient , see Eq.
double kap3
Pressure coefficient , see Eq.
virtual double nbar_ent_p2(const double ent2) const
Computes baryon density out of the log-enthalpy assuming that only fluid 2 is present.
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 ecart
contains the precision required in the relaxation nbar_ent_p
Eos_bf_poly(double kappa1, double kappa2, double kappa3, double beta)
Standard constructor.
double relax
Parameters needed for some inversions of the EOS.
virtual ~Eos_bf_poly()
Destructor.
2-fluids equation of state base class.
virtual void sauve(FILE *) const
Save in a file.
double m_1
Individual particle mass [unit: ].
virtual int identify() const =0
Returns a number to identify the sub-classe of Eos_bifluid the object belongs to.
void operator=(const Eos_bifluid &)
Assignment to another Eos_bifluid.
double m_2
Individual particle mass [unit: ].
Polytropic equation of state (relativistic case).
Equation of state base class.
void add_double(const double &x, int position=0)
Adds the the address of a new double to the list.
void add_double_mod(double &x, int position=0)
Adds the address of a new modifiable double to the list.
double & get_double_mod(int position=0) const
Returns the reference of a stored modifiable double .
Cmp sqrt(const Cmp &)
Square root.
Cmp exp(const Cmp &)
Exponential.
Cmp pow(const Cmp &, int)
Power .
Cmp log(const Cmp &)
Neperian logarithm.
int fread_be(int *aa, int size, int nb, FILE *fich)
Reads integer(s) from a binary file according to the big endian convention.
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.
int read_variable(const char *fname, const char *var_name, char *fmt, void *varp)
Reads a variable from file.
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.