29char 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"
104double puis(
double,
double) ;
105double enthal1(
const double x,
const Param& parent) ;
106double enthal23(
const double x,
const Param& parent) ;
107double 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)) ;
634 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
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.;
Analytic equation of state for two fluids (Newtonian case).
virtual Eos * trans2Eos() const
Makes a translation from Eos_bifluid to Eos .
virtual bool operator!=(const Eos_bifluid &) const
Comparison operator (difference)
void operator=(const Eos_bf_poly_newt &)
Assignment to another Eos_bf_poly_newt.
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.
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 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 ostream & operator>>(ostream &) const
Operator >>
virtual bool operator==(const Eos_bifluid &) const
Comparison operator (egality)
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.
Eos_bf_poly_newt(double kappa1, double kappa2, double kappa3, double beta)
Standard constructor.
virtual void sauve(FILE *) const
Save in a file.
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 ~Eos_bf_poly_newt()
Destructor.
virtual double get_K12(const double n1, const double n2, const double delta2) const
Computes the derivative of the energy with respect to .
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 press_nbar_p(const double nbar1, const double nbar2, const double delta2) const
Computes the pressure from the baryonic densities and the relative velocity.
Analytic equation of state for two fluids (relativistic case).
double gam5
Adiabatic indexes , see Eq.~eeosbfpolye}.
double kap1
Pressure coefficient , see Eq.
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}.
double gam6
Adiabatic indexes , see Eq.~eeosbfpolye}.
double gam2
Adiabatic indexes , see Eq.~eeosbfpolye}.
double precis
contains the precision required in zerosec_b
virtual void sauve(FILE *) const
Save in a file.
int typeos
The bi-fluid analytical EOS type:
double beta
Coefficient , see Eq.
double kap3
Pressure coefficient , see Eq.
double ecart
contains the precision required in the relaxation nbar_ent_p
double relax
Parameters needed for some inversions of the EOS.
2-fluids equation of state base class.
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.
double m_2
Individual particle mass [unit: ].
Polytropic equation of state (Newtonian 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 pow(const Cmp &, int)
Power .
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.