29 char eos_bf_poly_newt_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Eos/eos_bf_poly_newt.C,v 1.15 2014/10/13 08:52:52 j_novak Exp $" ;
95 #include "eos_bifluid.h" 98 #include "utilitaires.h" 104 double puis(
double,
double) ;
105 double enthal1(
const double x,
const Param& parent) ;
106 double enthal23(
const double x,
const Param& parent) ;
107 double enthal(
const double x,
const Param& parent) ;
120 name =
"bi-fluid polytropic non-relativistic EOS" ;
126 double gamma4,
double gamma5,
double gamma6,
127 double kappa1,
double kappa2,
double kappa3,
128 double bet,
double mass1,
double mass2,
129 double l_relax,
double l_precis,
double l_ecart) :
130 Eos_bf_poly(gamma1, gamma2, gamma3, gamma4, gamma5, gamma6,
131 kappa1, kappa2, kappa3, bet, mass1, mass2, l_relax, l_precis,
133 name =
"bi-fluid polytropic non-relativistic EOS" ;
181 cout <<
"The second EOS is not of type Eos_bf_poly_newt !" << endl ;
192 <<
"The two Eos_bf_poly_newt have different gammas : " <<
gam1 <<
" <-> " 193 << eos.
gam1 <<
", " <<
gam2 <<
" <-> " 194 << eos.
gam2 <<
", " <<
gam3 <<
" <-> " 195 << eos.
gam3 <<
", " <<
gam4 <<
" <-> " 196 << eos.
gam4 <<
", " <<
gam5 <<
" <-> " 197 << eos.
gam5 <<
", " <<
gam6 <<
" <-> " 198 << eos.
gam6 << endl ;
204 <<
"The two Eos_bf_poly_newt have different kappas : " <<
kap1 <<
" <-> " 205 << eos.
kap1 <<
", " <<
kap2 <<
" <-> " 206 << eos.
kap2 <<
", " <<
kap3 <<
" <-> " 207 << eos.
kap3 << endl ;
213 <<
"The two Eos_bf_poly_newt have different betas : " <<
beta <<
" <-> " 214 << eos.
beta << endl ;
220 <<
"The two Eos_bf_poly_newt have different masses : " <<
m_1 <<
" <-> " 221 << eos.
m_1 <<
", " <<
m_2 <<
" <-> " 250 ost <<
"EOS of class Eos_bf_poly_newt (non-relativistic polytrope) : " << endl ;
251 ost <<
" Adiabatic index gamma1 : " <<
gam1 << endl ;
252 ost <<
" Adiabatic index gamma2 : " <<
gam2 << endl ;
253 ost <<
" Adiabatic index gamma3 : " <<
gam3 << endl ;
254 ost <<
" Adiabatic index gamma4 : " <<
gam4 << endl ;
255 ost <<
" Adiabatic index gamma5 : " <<
gam5 << endl ;
256 ost <<
" Adiabatic index gamma6 : " <<
gam6 << endl ;
257 ost <<
" Pressure coefficient kappa1 : " <<
kap1 <<
258 " rho_nuc c^2 / n_nuc^gamma" << endl ;
259 ost <<
" Pressure coefficient kappa2 : " <<
kap2 <<
260 " rho_nuc c^2 / n_nuc^gamma" << endl ;
261 ost <<
" Pressure coefficient kappa3 : " <<
kap3 <<
262 " rho_nuc c^2 / n_nuc^gamma" << endl ;
263 ost <<
" Coefficient beta : " <<
beta <<
264 " rho_nuc c^2 / n_nuc^gamma" << endl ;
265 ost <<
" Mean particle 1 mass : " <<
m_1 <<
" m_B" << endl ;
266 ost <<
" Mean particle 2 mass : " <<
m_2 <<
" m_B" << endl ;
267 ost <<
" Parameters for inversion (used if typeos = 4 : " << endl ;
268 ost <<
" relaxation : " <<
relax << endl ;
269 ost <<
" precision for zerosec_b : " <<
precis << endl ;
270 ost <<
" final discrepancy in the result : " <<
ecart << endl ;
286 const double delta2,
double& nbar1,
287 double& nbar2)
const {
292 bool one_fluid =
false;
299 double determ =
kap1*
kap2 - kpd*kpd ;
301 nbar1 = (
kap2*ent1*
m_1 - kpd*ent2*
m_2) / determ ;
302 nbar2 = (
kap1*ent2*
m_2 - kpd*ent1*
m_1) / determ ;
303 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
307 double b1 = ent1*
m_1 ;
308 double b2 = ent2*
m_2 ;
310 if (fabs(denom) < 1.e-15) {
313 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
328 double f0 = enthal1(0.,parent) ;
329 if (fabs(f0)<1.e-15) {
333 assert (fabs(cas) > 1.e-15) ;
335 xmax = cas/fabs(cas) ;
338 if (fabs(xmax) > 1.e10) {
339 cout <<
"Problem in Eos_bf_poly::nbar_ent_p!" << endl ;
340 cout << f0 <<
", " << cas << endl ;
341 cout <<
"typeos = 1" << endl ;
344 }
while (enthal1(xmax,parent)*f0 > 0.) ;
345 double l_precis = 1.e-12 ;
348 nbar1 =
zerosec_b(enthal1, parent, xmin, xmax, l_precis,
351 nbar2 = (b1 - coef1*puis(nbar1,
gam1m1))/denom ;
352 double resu1 = coef1*puis(nbar1,
gam1m1) + denom*nbar2 ;
354 double erreur = fabs((resu1/m_1-ent1)/ent1) +
355 fabs((resu2/m_2-ent2)/ent2) ;
356 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)||(erreur > 1.e-4)) ;
361 double b1 = ent1*
m_1 ;
362 double b2 = ent2*
m_2 ;
364 if (fabs(denom) < 1.e-15) {
367 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
370 double coef0 =
beta*delta2 ;
373 double coef3 = 1./
gam1m1 ;
387 double f0 = enthal23(0.,parent) ;
388 if (fabs(f0)<1.e-15) {
391 double pmax = (fabs(
kap3) < 1.e-15 ? 0. :
gam4*(
gam4-1) ) ;
392 double ptemp = (fabs(
kap3*coef0) < 1.e-15 ? 0.
394 pmax = (pmax>ptemp ? pmax : ptemp) ;
395 ptemp = (fabs(
kap3*coef0) < 1.e-15 ? 0. :
gam4*(
gam6-1) ) ;
396 pmax = (pmax>ptemp ? pmax : ptemp) ;
397 ptemp = (fabs(coef0) < 1.e-15 ? 0. :
gam6*(
gam6-1) ) ;
398 pmax = (pmax>ptemp ? pmax : ptemp) ;
401 assert (fabs(cas) > 1.e-15) ;
403 xmax = cas/fabs(cas) ;
406 if (fabs(xmax) > 1.e10) {
407 cout <<
"Problem in Eos_bf_poly::nbar_ent_p!" << endl ;
408 cout <<
"typeos = 2" << endl ;
411 }
while (enthal23(xmax,parent)*f0 > 0.) ;
412 double l_precis = 1.e-12 ;
415 nbar2 =
zerosec_b(enthal23, parent, xmin, xmax, l_precis,
418 nbar1 = (b1 -
kap3*puis(nbar2,
gam4) - coef0*puis(nbar2,
gam6))
420 nbar1 = puis(nbar1,coef3) ;
422 + coef0*puis(nbar2,
gam6) ;
423 double resu2 = coef2*puis(nbar2,
gam2m1)
425 +
gam6*coef0*puis(nbar2,
gam6-1)*nbar1 ;
426 double erreur = fabs((resu1/m_1-ent1)/ent1) +
427 fabs((resu2/m_2-ent2)/ent2) ;
429 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)||(erreur > 1.e-4)) ;
434 double b1 = ent1*
m_1 ;
435 double b2 = ent2*
m_2 ;
437 if (fabs(denom) < 1.e-15) {
440 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
443 double coef0 =
beta*delta2 ;
446 double coef3 = 1./
gam2m1 ;
460 double f0 = enthal23(0.,parent) ;
461 if (fabs(f0)<1.e-15) {
464 double pmax = (fabs(
kap3) < 1.e-15 ? 0. :
gam3*(
gam3-1) ) ;
465 double ptemp = (fabs(
kap3*coef0) < 1.e-15 ? 0.
467 pmax = (pmax>ptemp ? pmax : ptemp) ;
468 ptemp = (fabs(
kap3*coef0) < 1.e-15 ? 0. :
gam3*(
gam5-1) ) ;
469 pmax = (pmax>ptemp ? pmax : ptemp) ;
470 ptemp = (fabs(coef0) < 1.e-15 ? 0. :
gam5*(
gam5-1) ) ;
471 pmax = (pmax>ptemp ? pmax : ptemp) ;
474 assert (fabs(cas) > 1.e-15) ;
476 xmax = cas/fabs(cas) ;
479 if (fabs(xmax) > 1.e10) {
480 cout <<
"Problem in Eos_bf_poly::nbar_ent_p!" << endl ;
481 cout <<
"typeos = 3" << endl ;
484 }
while (enthal23(xmax,parent)*f0 > 0.) ;
485 double l_precis = 1.e-12 ;
488 nbar1 =
zerosec_b(enthal23, parent, xmin, xmax, l_precis,
491 nbar2 = (b2 -
kap3*puis(nbar1,
gam3) - coef0*puis(nbar1,
gam5))
493 nbar2 = puis(nbar2,coef3) ;
494 double resu1 = coef1*puis(nbar1,
gam1m1)
496 + coef0*
gam5*puis(nbar1,
gam5-1)*nbar2 ;
497 double resu2 = coef2*puis(nbar2,
gam2m1)
499 double erreur = fabs((resu1/m_1-ent1)/ent1) +
500 fabs((resu2/m_2-ent2)/ent2) ;
501 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)||(erreur > 1.e-4)) ;
506 double b1 = ent1*
m_1 ;
507 double b2 = ent2*
m_2 ;
509 if (fabs(denom) < 1.e-15) {
512 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
527 double n1l, n2l, n1s, n2s ;
529 double delta = a11*a22 - (a13+a14)*(a23+a24) ;
530 n1l = (a22*b1 - (a13+a14)*b2)/delta ;
531 n2l = (a11*b2 - (a23+a24)*b1)/delta ;
532 n1s = puis(b1/a11, 1./(
gam1m1)) ;
533 n2s = puis(b2/a22, 1./(
gam2m1)) ;
535 double n1m = (n1l<0. ? n1s :
sqrt(n1l*n1s)) ;
536 double n2m = (n2l<0. ? n2s :
sqrt(n2l*n2s)) ;
540 double c1, c2, c3, c4, c5, c6, c7 ;
545 c5 = a13*puis(n2m,
gam4) ;
546 c6 = a14*puis(n2m,
gam6) ;
556 double d1, d2, d3, d4, d5, d6, d7 ;
561 d5 = a23*puis(n1m,
gam3) ;
562 d6 = a24*puis(n1m,
gam5) ;
572 double xmin = -3*(n1s>n2s ? n1s : n2s) ;
573 double xmax = 3*(n1s>n2s ? n1s : n2s) ;
588 sortie = (fabs((n1m-n1)/(n1m+n1)) + fabs((n2m-n2)/(n2m+n2)) >
ecart) ;
610 }
while (sortie&&(mer<nmermax)) ;
632 nbar1 = (
kap2*(ent1*
m_1 -
beta*delta2) - kap3*ent2*
m_2) / determ ;
633 nbar2 = (
kap1*ent2*m_2 - kap3*(ent1*
m_1 -
beta*delta2)) / determ ;
634 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
640 nbar1 = (
kap2*ent1*
m_1 - kap3*(ent2*
m_2 -
beta*delta2)) / determ ;
641 nbar2 = (
kap1*(ent2*
m_2 -
beta*delta2) - kap3*ent1*
m_1) / determ ;
642 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
665 const double delta2)
const {
667 double n1 = (nbar1>double(0) ? nbar1 : 0) ;
668 double n2 = (nbar2>double(0) ? nbar2 : 0) ;
669 double x2 = ((nbar1>double(0))&&(nbar2>
double(0))) ? delta2 : 0 ;
683 const double delta2)
const {
685 double n1 = (nbar1>double(0) ? nbar1 : 0) ;
686 double n2 = (nbar2>double(0) ? nbar2 : 0) ;
687 double x2 = ((nbar1>double(0))&&(nbar2>
double(0))) ? delta2 : 0 ;
704 if (n1 <= 0.) xx = 0. ;
714 if (n2 <= 0.) xx = 0. ;
724 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}.