LORENE
eos_bf_poly.C
1 /*
2  * Methods of the class Eos_bf_poly.
3  *
4  * (see file eos_bifluid.h for documentation).
5  */
6 
7 /*
8  * Copyright (c) 2001-2002 Jerome Novak
9  *
10  * This file is part of LORENE.
11  *
12  * LORENE is free software; you can redistribute it and/or modify
13  * it under the terms of the GNU General Public License as published by
14  * the Free Software Foundation; either version 2 of the License, or
15  * (at your option) any later version.
16  *
17  * LORENE is distributed in the hope that it will be useful,
18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20  * GNU General Public License for more details.
21  *
22  * You should have received a copy of the GNU General Public License
23  * along with LORENE; if not, write to the Free Software
24  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25  *
26  */
27 
28 
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 $" ;
30 
31 /*
32  * $Id: eos_bf_poly.C,v 1.22 2014/10/13 08:52:52 j_novak Exp $
33  * $Log: eos_bf_poly.C,v $
34  * Revision 1.22 2014/10/13 08:52:52 j_novak
35  * Lorene classes and functions now belong to the namespace Lorene.
36  *
37  * Revision 1.21 2014/10/06 15:13:06 j_novak
38  * Modified #include directives to use c++ syntax.
39  *
40  * Revision 1.20 2014/04/25 10:43:51 j_novak
41  * The member 'name' is of type string now. Correction of a few const-related issues.
42  *
43  * Revision 1.19 2008/08/19 06:42:00 j_novak
44  * Minor modifications to avoid warnings with gcc 4.3. Most of them concern
45  * cast-type operations, and constant strings that must be defined as const char*
46  *
47  * Revision 1.18 2004/05/13 15:27:42 r_prix
48  * fixed a little eos-bug also in the relativistic case (same as done already in Newt)
49  *
50  * Revision 1.17 2003/12/17 23:12:32 r_prix
51  * replaced use of C++ <string> by standard ANSI char* to be backwards compatible
52  * with broken compilers like MIPSpro Compiler 7.2 on SGI Origin200. ;-)
53  *
54  * Revision 1.16 2003/12/10 08:58:20 r_prix
55  * - added new Eos_bifluid paramter for eos-file: bool slow_rot_style
56  * to indicate if we want this particular kind of EOS-inversion (only works for
57  * the Newtonian 'analytic' polytropes) --> replaces former dirty hack with gamma1<0
58  *
59  * Revision 1.15 2003/12/05 15:09:47 r_prix
60  * adapted Eos_bifluid class and subclasses to use read_variable() for
61  * (formatted) file-reading.
62  *
63  * Revision 1.14 2003/12/04 14:24:41 r_prix
64  * a really dirty hack: gam1 < 0 signals to use 'slow-rot-style' EOS
65  * inversion (i.e. typeos=5). This only works for the 'analytic EOS'.
66  *
67  * Revision 1.13 2003/11/18 18:28:38 r_prix
68  * moved particle-masses m_1, m_2 of the two fluids into class eos_bifluid (from eos_bf_poly)
69  *
70  * Revision 1.12 2003/03/17 10:28:04 j_novak
71  * Removed too strong asserts on beta and kappa
72  *
73  * Revision 1.11 2003/02/07 10:47:43 j_novak
74  * The possibility of having gamma5 xor gamma6 =0 has been introduced for
75  * tests. The corresponding parameter files have been added.
76  *
77  * Revision 1.10 2002/10/16 14:36:34 j_novak
78  * Reorganization of #include instructions of standard C++, in order to
79  * use experimental version 3 of gcc.
80  *
81  * Revision 1.9 2002/05/31 16:13:36 j_novak
82  * better inversion for eos_bifluid
83  *
84  * Revision 1.8 2002/05/10 09:26:52 j_novak
85  * Added new class Et_rot_mag for magnetized rotating neutron stars (under development)
86  *
87  * Revision 1.7 2002/05/02 15:16:22 j_novak
88  * Added functions for more general bi-fluid EOS
89  *
90  * Revision 1.6 2002/04/05 09:09:36 j_novak
91  * The inversion of the EOS for 2-fluids polytrope has been modified.
92  * Some errors in the determination of the surface were corrected.
93  *
94  * Revision 1.5 2002/01/16 15:03:28 j_novak
95  * *** empty log message ***
96  *
97  * Revision 1.4 2002/01/11 14:09:34 j_novak
98  * Added newtonian version for 2-fluid stars
99  *
100  * Revision 1.3 2001/12/04 21:27:53 e_gourgoulhon
101  *
102  * All writing/reading to a binary file are now performed according to
103  * the big endian convention, whatever the system is big endian or
104  * small endian, thanks to the functions fwrite_be and fread_be
105  *
106  * Revision 1.2 2001/11/29 15:05:26 j_novak
107  * The entrainment term in 2-fluid eos is modified
108  *
109  * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
110  * LORENE
111  *
112  * Revision 1.4 2001/08/31 15:48:50 novak
113  * The flag tronc has been added to nbar_ent_p
114  *
115  * Revision 1.3 2001/08/27 09:52:21 novak
116  * New version of formulas
117  *
118  * Revision 1.2 2001/06/22 15:36:46 novak
119  * Modification de trans2Eos
120  *
121  * Revision 1.1 2001/06/21 15:24:46 novak
122  * Initial revision
123  *
124  *
125  * $Header: /cvsroot/Lorene/C++/Source/Eos/eos_bf_poly.C,v 1.22 2014/10/13 08:52:52 j_novak Exp $
126  *
127  */
128 
129 // Headers C
130 #include <cstdlib>
131 #include <cmath>
132 
133 // Headers Lorene
134 #include "eos_bifluid.h"
135 #include "param.h"
136 #include "eos.h"
137 #include "cmp.h"
138 #include "utilitaires.h"
139 
140 namespace Lorene {
141 
142 //****************************************************************************
143 // local prototypes
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) ;
148 //****************************************************************************
149 
150  //--------------//
151  // Constructors //
152  //--------------//
153 
154 // Standard constructor with gam1 = gam2 = 2,
155 // gam3 = gam4 = gam5 = gam6 = 1, m_1 = 1 and m_2 =1
156 // -------------------------------------------------
157 Eos_bf_poly::Eos_bf_poly(double kappa1, double kappa2, double kappa3, double bet) :
158  Eos_bifluid("bi-fluid polytropic EOS", 1, 1),
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)
162 {
163 
164  determine_type() ;
165  set_auxiliary() ;
166 
167 }
168 
169 // Standard constructor with everything specified
170 // -----------------------------------------------
171 Eos_bf_poly::Eos_bf_poly(double gamma1, double gamma2, double gamma3,
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)
180 {
181 
182  determine_type() ;
183  set_auxiliary() ;
184 
185 }
186 
187 // Copy constructor
188 // ----------------
190  Eos_bifluid(eosi),
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),
194  beta(eosi.beta),
195  typeos(eosi.typeos), relax(eosi.relax), precis(eosi.precis),
196  ecart(eosi.ecart) {
197 
198  set_auxiliary() ;
199 
200 }
201 
202 
203 // Constructor from binary file
204 // ----------------------------
206  Eos_bifluid(fich) {
207 
208  fread_be(&gam1, sizeof(double), 1, fich) ;
209  fread_be(&gam2, sizeof(double), 1, fich) ;
210  fread_be(&gam3, sizeof(double), 1, fich) ;
211  fread_be(&gam4, sizeof(double), 1, fich) ;
212  fread_be(&gam5, sizeof(double), 1, fich) ;
213  fread_be(&gam6, sizeof(double), 1, fich) ;
214  fread_be(&kap1, sizeof(double), 1, fich) ;
215  fread_be(&kap2, sizeof(double), 1, fich) ;
216  fread_be(&kap3, sizeof(double), 1, fich) ;
217  fread_be(&beta, sizeof(double), 1, fich) ;
218  fread_be(&relax, sizeof(double), 1, fich) ;
219  fread_be(&precis, sizeof(double), 1, fich) ;
220  fread_be(&ecart, sizeof(double), 1, fich) ;
221 
222  determine_type() ;
223  set_auxiliary() ;
224 
225 
226 }
227 
228 
229 // Constructor from a formatted file
230 // ---------------------------------
231 Eos_bf_poly::Eos_bf_poly(const char *fname ) :
232  Eos_bifluid(fname), relax(0.5), precis(1.e-8), ecart(1.e-7)
233 {
234  int res = 0;
235 
236  res += read_variable (fname, const_cast<char*>("gamma1"), gam1);
237  res += read_variable (fname, const_cast<char*>("gamma2"), gam2);
238  res += read_variable (fname, const_cast<char*>("gamma3"), gam3);
239  res += read_variable (fname, const_cast<char*>("gamma4"), gam4);
240  res += read_variable (fname, const_cast<char*>("gamma5"), gam5);
241  res += read_variable (fname, const_cast<char*>("gamma6"), gam6);
242  res += read_variable (fname, const_cast<char*>("kappa1"), kap1);
243  res += read_variable (fname, const_cast<char*>("kappa2"), kap2);
244  res += read_variable (fname, const_cast<char*>("kappa3"), kap3);
245  res += read_variable (fname, const_cast<char*>("beta"), beta);
246 
247 
248  if (res != 0)
249  {
250  cerr << "ERROR: could not read all required eos-paramters for Eos_bf_poly()" << endl;
251  exit (-1);
252  }
253 
254  determine_type() ;
255 
256  if (get_typeos() == 4)
257  {
258  res += read_variable (fname, const_cast<char*>("relax"), relax);
259  res += read_variable (fname, const_cast<char*>("precis"), precis);
260  res += read_variable (fname, const_cast<char*>("ecart"), ecart);
261  }
262  else if (get_typeos() == 0) // analytic EOS: check if we want slow-rot-style inversion
263  {
264  bool slowrot = false;
265  read_variable (fname, const_cast<char*>("slow_rot_style"), slowrot); // dont require this variable!
266  if (slowrot)
267  typeos = 5; // type=5 is reserved for (type0 + slow-rot-style)
268  }
269 
270 
271  if (res != 0)
272  {
273  cerr << "ERROR: could not read all required eos-paramters for Eos_bf_poly()" << endl;
274  exit (-1);
275  }
276 
277 
278  set_auxiliary() ;
279 
280 }
281  //--------------//
282  // Destructor //
283  //--------------//
284 
286 
287  //Does nothing ...
288 
289 }
290  //--------------//
291  // Assignment //
292  //--------------//
293 
295 
296  // Assignment of quantities common to all the derived classes of Eos_bifluid
297  Eos_bifluid::operator=(eosi) ;
298 
299  gam1 = eosi.gam1 ;
300  gam2 = eosi.gam2 ;
301  gam3 = eosi.gam3 ;
302  kap1 = eosi.kap1 ;
303  kap2 = eosi.kap2 ;
304  kap3 = eosi.kap3 ;
305  beta = eosi.beta ;
306  typeos = eosi.typeos ;
307  relax = eosi.relax ;
308  precis = eosi.precis ;
309  ecart = eosi.ecart ;
310 
311  set_auxiliary() ;
312 
313 }
314 
315 
316  //-----------------------//
317  // Miscellaneous //
318  //-----------------------//
319 
321 
322  gam1m1 = gam1 - double(1) ;
323  gam2m1 = gam2 - double(1) ;
324  gam34m1 = gam3 + gam4 - double(1) ;
325  gam56m1 = gam5 + gam6 - double(1) ;
326 
327  if (fabs(kap3*kap3-kap2*kap1) < 1.e-15) {
328  cout << "WARNING!: Eos_bf_poly: the parameters are degenerate!" << endl ;
329  abort() ;
330  }
331 
332 }
333 
335 
336  if ((gam1 == double(2)) && (gam2 == double(2)) && (gam3 == double(1))
337  && (gam4 == double(1)) && (gam5 == double(1))
338  && (gam6 == double(0))) {
339  typeos = -1 ;
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 ;
345  }
346  else if ((gam1 == double(2)) && (gam2 == double(2)) && (gam3 == double(1))
347  && (gam4 == double(1)) && (gam5 == double(0))
348  && (gam6 == double(1))) {
349  typeos = -2 ;
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 ;
355  }
356  else if ((gam1 == double(2)) && (gam2 == double(2)) && (gam3 == double(1))
357  && (gam4 == double(1)) && (gam5 == double(1))
358  && (gam6 == double(1))) {
359  typeos = 0 ;
360  }
361  else if ((gam3 == double(1)) && (gam4 == double(1)) && (gam5 == double(1))
362  && (gam6 == double(1))) {
363  typeos = 1 ;
364  }
365  else if ((gam3 == double(1)) && (gam5 == double(1))) {
366  typeos = 2 ;
367  }
368  else if ((gam4 == double(1)) && (gam6 == double(1))) {
369  typeos = 3 ;
370  return ;
371  }
372  else {
373  typeos = 4 ;
374  }
375 
376  cout << "DEBUG: EOS-type was determined as typeos = " << typeos << endl;
377 
378  return ;
379 }
380 
381  //------------------------//
382  // Comparison operators //
383  //------------------------//
384 
385 
386 bool Eos_bf_poly::operator==(const Eos_bifluid& eos_i) const {
387 
388  bool resu = true ;
389 
390  if ( eos_i.identify() != identify() ) {
391  cout << "The second EOS is not of type Eos_bf_poly !" << endl ;
392  resu = false ;
393  }
394  else{
395 
396  const Eos_bf_poly& eos = dynamic_cast<const Eos_bf_poly&>( eos_i ) ;
397 
398  if ((eos.gam1 != gam1)||(eos.gam2 != gam2)||(eos.gam3 != gam3)
399  ||(eos.gam4 != gam4)||(eos.gam5 != gam5)||(eos.gam6 != gam6)) {
400  cout
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 ;
408  resu = false ;
409  }
410 
411  if ((eos.kap1 != kap1)||(eos.kap2 != kap2)|| (eos.kap3 != kap3)){
412  cout
413  << "The two Eos_bf_poly have different kappas : " << kap1 << " <-> "
414  << eos.kap1 << ", " << kap2 << " <-> "
415  << eos.kap2 << ", " << kap3 << " <-> "
416  << eos.kap3 << endl ;
417  resu = false ;
418  }
419 
420  if (eos.beta != beta) {
421  cout
422  << "The two Eos_bf_poly have different betas : " << beta << " <-> "
423  << eos.beta << endl ;
424  resu = false ;
425  }
426 
427  if ((eos.m_1 != m_1)||(eos.m_2 != m_2)) {
428  cout
429  << "The two Eos_bf_poly have different masses : " << m_1 << " <-> "
430  << eos.m_1 << ", " << m_2 << " <-> "
431  << eos.m_2 << endl ;
432  resu = false ;
433  }
434 
435  }
436 
437  return resu ;
438 
439 }
440 
441 bool Eos_bf_poly::operator!=(const Eos_bifluid& eos_i) const {
442 
443  return !(operator==(eos_i)) ;
444 
445 }
446 
447 
448  //------------//
449  // Outputs //
450  //------------//
451 
452 void Eos_bf_poly::sauve(FILE* fich) const {
453 
454  Eos_bifluid::sauve(fich) ;
455 
456  fwrite_be(&gam1, sizeof(double), 1, fich) ;
457  fwrite_be(&gam2, sizeof(double), 1, fich) ;
458  fwrite_be(&gam3, sizeof(double), 1, fich) ;
459  fwrite_be(&gam4, sizeof(double), 1, fich) ;
460  fwrite_be(&gam5, sizeof(double), 1, fich) ;
461  fwrite_be(&gam6, sizeof(double), 1, fich) ;
462  fwrite_be(&kap1, sizeof(double), 1, fich) ;
463  fwrite_be(&kap2, sizeof(double), 1, fich) ;
464  fwrite_be(&kap3, sizeof(double), 1, fich) ;
465  fwrite_be(&beta, sizeof(double), 1, fich) ;
466  fwrite_be(&relax, sizeof(double), 1, fich) ;
467  fwrite_be(&precis, sizeof(double), 1, fich) ;
468  fwrite_be(&ecart, sizeof(double), 1, fich) ;
469 
470 }
471 
472 ostream& Eos_bf_poly::operator>>(ostream & ost) const {
473 
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 ;
494 
495  return ost ;
496 
497 }
498 
499 
500  //------------------------------//
501  // Computational routines //
502  //------------------------------//
503 
504 // Baryon densities from enthalpies
505 //---------------------------------
506 
507 bool Eos_bf_poly::nbar_ent_p(const double ent1, const double ent2,
508  const double delta2, double& nbar1,
509  double& nbar2) const {
510 
511  // RP: follow Newtonian case in this: the following is wrong, I think
512  // bool one_fluid = ((ent1<=0.)||(ent2<=0.)) ;
513 
514  bool one_fluid = false;
515 
516  if (!one_fluid) {
517  switch (typeos) {
518  case 5: // same as typeos=0 but with slow-rot-style inversion
519  case 0: {
520  double kpd = kap3+beta*delta2 ;
521  double determ = kap1*kap2 - kpd*kpd ;
522 
523  nbar1 = (kap2*(exp(ent1) - m_1) - kpd*(exp(ent2) - m_2)) / determ ;
524  nbar2 = (kap1*(exp(ent2) - m_2) - kpd*(exp(ent1) - m_1)) / determ ;
525  one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
526  break ;
527  }
528  case 1: {
529  double b1 = exp(ent1) - m_1 ;
530  double b2 = exp(ent2) - m_2 ;
531  double denom = kap3 + beta*delta2 ;
532  if (fabs(denom) < 1.e-15) {
533  nbar1 = puis(2*b1/(gam1*kap1), 1./gam1m1) ;
534  nbar2 = puis(2*b2/(gam2*kap2), 1./gam2m1) ;
535  one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
536  }
537  else {
538  double coef0 = 0.5*gam2*kap2/pow(denom, gam2m1) ;
539  double coef1 = 0.5*kap1*gam1 ;
540  Param parent ;
541  parent.add_double(coef0, 0) ;
542  parent.add_double(b1, 1) ;
543  parent.add_double(coef1, 2) ;
544  parent.add_double(gam1m1,3) ;
545  parent.add_double(gam2m1,4) ;
546  parent.add_double(denom, 5) ;
547  parent.add_double(b2, 6) ;
548 
549  double xmin, xmax ;
550  double f0 = enthal1(0.,parent) ;
551  if (fabs(f0)<1.e-15) {
552  nbar1 = 0. ;}
553  else {
554  double cas = (gam1m1*gam2m1 - 1.)*f0;
555  assert (fabs(cas) > 1.e-15) ;
556  xmin = 0. ;
557  xmax = cas/fabs(cas) ;
558  do {
559  xmax *= 1.1 ;
560  if (fabs(xmax) > 1.e10) {
561  cout << "Problem in Eos_bf_poly::nbar_ent_p!" << endl ;
562  cout << f0 << ", " << cas << endl ; // to be removed!!
563  cout << "typeos = 1" << endl ;
564  abort() ;
565  }
566  } while (enthal1(xmax,parent)*f0 > 0.) ;
567  double l_precis = 1.e-12 ;
568  int nitermax = 400 ;
569  int niter = 0 ;
570  nbar1 = zerosec_b(enthal1, parent, xmin, xmax, l_precis,
571  nitermax, niter) ;
572  }
573  nbar2 = (b1 - coef1*puis(nbar1, gam1m1))/denom ;
574  double resu1 = coef1*puis(nbar1,gam1m1) + denom*nbar2 ;
575  double resu2 = 0.5*gam2*kap2*puis(nbar2,gam2m1) + denom*nbar1 ;
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)) ;
579  }
580  break ;
581  }
582  case 2: {
583  double b1 = exp(ent1) - m_1 ;
584  double b2 = exp(ent2) - m_2 ;
585  double denom = kap3 + beta*delta2 ;
586  if (fabs(denom) < 1.e-15) {
587  nbar1 = puis(2*b1/(gam1*kap1), 1./gam1m1) ;
588  nbar2 = puis(2*b2/(gam2*kap2), 1./gam2m1) ;
589  one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
590  }
591  else {
592  double coef0 = beta*delta2 ;
593  double coef1 = 0.5*kap1*gam1 ;
594  double coef2 = 0.5*kap2*gam2 ;
595  double coef3 = 1./gam1m1 ;
596  Param parent ;
597  parent.add_double(b1, 0) ;
598  parent.add_double(kap3, 1) ;
599  parent.add_double(gam4, 2) ;
600  parent.add_double(coef0, 3) ;
601  parent.add_double(gam6,4) ;
602  parent.add_double(coef1, 5) ;
603  parent.add_double(coef3, 6) ;
604  parent.add_double(coef2, 7) ;
605  parent.add_double(gam2m1, 8) ;
606  parent.add_double(b2, 9) ;
607 
608  double xmin, xmax ;
609  double f0 = enthal23(0.,parent) ;
610  if (fabs(f0)<1.e-15) {
611  nbar2 = 0. ;}
612  else {
613  double pmax = (fabs(kap3) < 1.e-15 ? 0. : gam4*(gam4-1) ) ;
614  double ptemp = (fabs(kap3*coef0) < 1.e-15 ? 0.
615  : gam6*(gam4-1) ) ;
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) ;
621  double cas = (pmax - gam1m1*gam2m1)*f0;
622  // cout << "Pmax, cas: " << pmax << ", " << cas << endl ;
623  assert (fabs(cas) > 1.e-15) ;
624  xmin = 0. ;
625  xmax = cas/fabs(cas) ;
626  do {
627  xmax *= 1.1 ;
628  if (fabs(xmax) > 1.e10) {
629  cout << "Problem in Eos_bf_poly::nbar_ent_p!" << endl ;
630  cout << "typeos = 2" << endl ;
631  abort() ;
632  }
633  } while (enthal23(xmax,parent)*f0 > 0.) ;
634  double l_precis = 1.e-12 ;
635  int nitermax = 400 ;
636  int niter = 0 ;
637  nbar2 = zerosec_b(enthal23, parent, xmin, xmax, l_precis,
638  nitermax, niter) ;
639  }
640  nbar1 = (b1 - kap3*puis(nbar2,gam4) - coef0*puis(nbar2,gam6))
641  / coef1 ;
642  nbar1 = puis(nbar1,coef3) ;
643  double resu1 = coef1*puis(nbar1,gam1m1) + kap3*puis(nbar2,gam4)
644  + coef0*puis(nbar2, gam6) ;
645  double resu2 = coef2*puis(nbar2,gam2m1)
646  + gam4*kap3*puis(nbar2, gam4-1)*nbar1
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) ;
650  //cout << "Erreur d'inversion: " << erreur << endl ;
651  one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)||(erreur > 1.e-4)) ;
652  }
653  break ;
654  }
655  case 3: {
656  double b1 = exp(ent1) - m_1 ;
657  double b2 = exp(ent2) - m_2 ;
658  double denom = kap3 + beta*delta2 ;
659  if (fabs(denom) < 1.e-15) {
660  nbar1 = puis(2*b1/(gam1*kap1), 1./gam1m1) ;
661  nbar2 = puis(2*b2/(gam2*kap2), 1./gam2m1) ;
662  one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
663  }
664  else {
665  double coef0 = beta*delta2 ;
666  double coef1 = 0.5*kap1*gam1 ;
667  double coef2 = 0.5*kap2*gam2 ;
668  double coef3 = 1./gam2m1 ;
669  Param parent ;
670  parent.add_double(b2, 0) ;
671  parent.add_double(kap3, 1) ;
672  parent.add_double(gam3, 2) ;
673  parent.add_double(coef0, 3) ;
674  parent.add_double(gam5, 4) ;
675  parent.add_double(coef2, 5) ;
676  parent.add_double(coef3, 6) ;
677  parent.add_double(coef1, 7) ;
678  parent.add_double(gam1m1, 8) ;
679  parent.add_double(b1, 9) ;
680 
681  double xmin, xmax ;
682  double f0 = enthal23(0.,parent) ;
683  if (fabs(f0)<1.e-15) {
684  nbar1 = 0. ;}
685  else {
686  double pmax = (fabs(kap3) < 1.e-15 ? 0. : gam3*(gam3-1) ) ;
687  double ptemp = (fabs(kap3*coef0) < 1.e-15 ? 0.
688  : gam5*(gam3-1) ) ;
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) ;
694  double cas = (pmax - gam1m1*gam2m1)*f0;
695  // cout << "Pmax, cas: " << pmax << ", " << cas << endl ;
696  assert (fabs(cas) > 1.e-15) ;
697  xmin = 0. ;
698  xmax = cas/fabs(cas) ;
699  do {
700  xmax *= 1.1 ;
701  if (fabs(xmax) > 1.e10) {
702  cout << "Problem in Eos_bf_poly::nbar_ent_p!" << endl ;
703  cout << "typeos = 3" << endl ;
704  abort() ;
705  }
706  } while (enthal23(xmax,parent)*f0 > 0.) ;
707  double l_precis = 1.e-12 ;
708  int nitermax = 400 ;
709  int niter = 0 ;
710  nbar1 = zerosec_b(enthal23, parent, xmin, xmax, l_precis,
711  nitermax, niter) ;
712  }
713  nbar2 = (b2 - kap3*puis(nbar1,gam3) - coef0*puis(nbar1,gam5))
714  / coef2 ;
715  nbar2 = puis(nbar2,coef3) ;
716  double resu1 = coef1*puis(nbar1,gam1m1)
717  + gam3*kap3*puis(nbar1,gam3-1)*nbar2
718  + coef0*gam5*puis(nbar1, gam5-1)*nbar2 ;
719  double resu2 = coef2*puis(nbar2,gam2m1)
720  + kap3*puis(nbar1, gam3) + coef0*puis(nbar1, gam5);
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)) ;
724  }
725  break ;
726  }
727  case 4:{
728  double b1 = exp(ent1) - m_1 ;
729  double b2 = exp(ent2) - m_2 ;
730  double denom = kap3 + beta*delta2 ;
731  if (fabs(denom) < 1.e-15) {
732  nbar1 = puis(2*b1/(gam1*kap1), 1./gam1m1) ;
733  nbar2 = puis(2*b2/(gam2*kap2), 1./gam2m1) ;
734  one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
735  }
736  else {
737  int nitermax = 200 ;
738  int niter = 0 ;
739  int nmermax = 800 ;
740 
741  double a11 = 0.5*gam1*kap1 ;
742  double a13 = gam3*kap3 ;
743  double a14 = beta*gam5*delta2 ;
744 
745  double a22 = 0.5*gam2*kap2 ;
746  double a23 = gam4*kap3 ;
747  double a24 = beta*gam6*delta2 ;
748 
749  double n1l, n2l, n1s, n2s ;
750 
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)) ;
756 
757  double n1m = (n1l<0. ? n1s : sqrt(n1l*n1s)) ;
758  double n2m = (n2l<0. ? n2s : sqrt(n2l*n2s)) ;
759 
760  Param parf1 ;
761  Param parf2 ;
762  double c1, c2, c3, c4, c5, c6, c7 ;
763  c1 = gam1m1 ;
764  c2 = gam3 - 1. ;
765  c3 = gam5 - 1. ;
766  c4 = a11 ;
767  c5 = a13*puis(n2m,gam4) ;
768  c6 = a14*puis(n2m,gam6) ;
769  c7 = b1 ;
770  parf1.add_double_mod(c1,0) ;
771  parf1.add_double_mod(c2,1) ;
772  parf1.add_double_mod(c3,2) ;
773  parf1.add_double_mod(c4,3) ;
774  parf1.add_double_mod(c5,4) ;
775  parf1.add_double_mod(c6,5) ;
776  parf1.add_double_mod(c7,6) ;
777 
778  double d1, d2, d3, d4, d5, d6, d7 ;
779  d1 = gam2m1 ;
780  d2 = gam4 - 1. ;
781  d3 = gam6 - 1. ;
782  d4 = a22 ;
783  d5 = a23*puis(n1m,gam3) ;
784  d6 = a24*puis(n1m,gam5) ;
785  d7 = b2 ;
786  parf2.add_double_mod(d1,0) ;
787  parf2.add_double_mod(d2,1) ;
788  parf2.add_double_mod(d3,2) ;
789  parf2.add_double_mod(d4,3) ;
790  parf2.add_double_mod(d5,4) ;
791  parf2.add_double_mod(d6,5) ;
792  parf2.add_double_mod(d7,6) ;
793 
794  double xmin = -3*(n1s>n2s ? n1s : n2s) ;
795  double xmax = 3*(n1s>n2s ? n1s : n2s) ;
796 
797  double n1 = n1m ;
798  double n2 = n2m ;
799  bool sortie = true ;
800  int mer = 0 ;
801 
802  //cout << "Initial guess: " << n1 << ", " << n2 << endl ;
803  n1s *= 0.1 ;
804  n2s *= 0.1 ;
805  do {
806  //cout << "n1, n2: " << n1 << ", " << n2 << endl ;
807  n1 = zerosec_b(&enthal, parf1, xmin, xmax, precis, nitermax, niter) ;
808  n2 = zerosec_b(&enthal, parf2, xmin, xmax, precis, nitermax, niter) ;
809 
810  sortie = (fabs((n1m-n1)/(n1m+n1)) + fabs((n2m-n2)/(n2m+n2)) > ecart) ;
811  n1m = relax*n1 + (1.-relax)*n1m ;
812  n2m = relax*n2 + (1.-relax)*n2m ;
813  if (n2m>-n2s) {
814  parf1.get_double_mod(4) = a13*puis(n2m,gam4) ;
815  parf1.get_double_mod(5) = a14*puis(n2m,gam6) ;
816  }
817  else {
818  parf1.get_double_mod(4) = a13*puis(-n2s,gam4) ;
819  parf1.get_double_mod(5) = a14*puis(-n2s,gam6) ;
820  }
821 
822  if (n1m>-n1s) {
823  parf2.get_double_mod(4) = a23*puis(n1m,gam3) ;
824  parf2.get_double_mod(5) = a24*puis(n1m,gam5) ;
825  }
826  else {
827  parf2.get_double_mod(4) = a23*puis(-n1s,gam3) ;
828  parf2.get_double_mod(5) = a24*puis(-n1s,gam5) ;
829  }
830 
831  mer++ ;
832  } while (sortie&&(mer<nmermax)) ;
833  nbar1 = n1m ;
834  nbar2 = n2m ;
835 
836 // double resu1 = a11*puis(n1,gam1m1) + a13*puis(n1,gam3-1.)*puis(n2,gam4)
837 // +a14*puis(n1,gam5-1.)*puis(n2,gam6) ;
838 // double resu2 = a22*puis(n2,gam2m1) + a23*puis(n1,gam3)*puis(n2,gam4-1.)
839 // +a24*puis(n1,gam5)*puis(n2,gam6-1.) ;
840  //cout << "Nbre d'iterations: " << mer << endl ;
841  //cout << "Resus: " << n1m << ", " << n2m << endl ;
842  //cout << "Verification: " << log(fabs(1+resu1)) << ", "
843  // << log(fabs(1+resu2)) << endl ;
844  // cout << "Erreur: " << fabs(enthal(n1, parf1)/b1) +
845  // fabs(enthal(n2, parf2)/b2) << endl ;
846  //cout << "Erreur: " << fabs((log(fabs(1+resu1))-ent1)/ent1) +
847  //fabs((log(fabs(1+resu2))-ent2)/ent2) << endl ;
848  }
849  break ;
850  }
851  case -1: {
852  double determ = kap1*kap2 - kap3*kap3 ;
853 
854  nbar1 = (kap2*(exp(ent1) - m_1 - beta*delta2)
855  - kap3*(exp(ent2) - m_2)) / determ ;
856  nbar2 = (kap1*(exp(ent2) - m_2) - kap3*(exp(ent1) - m_1 - beta*delta2))
857  / determ ;
858  one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
859  break ;
860  }
861  case -2: {
862  double determ = kap1*kap2 - kap3*kap3 ;
863 
864  nbar1 = (kap2*(exp(ent1) - m_1 )
865  - kap3*(exp(ent2) - m_2 - beta*delta2)) / determ ;
866  nbar2 = (kap1*(exp(ent2) - m_2 - beta*delta2)
867  - kap3*(exp(ent1) - m_1)) / determ ;
868  one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
869  break ;
870  }
871  }
872  }
873  return one_fluid ;
874 }
875 // One fluid sub-EOSs
876 //-------------------
877 
878 double Eos_bf_poly::nbar_ent_p1(const double ent1) const {
879  return puis(2*(exp(ent1) - m_1)/(gam1*kap1), 1./gam1m1) ;
880 }
881 
882 double Eos_bf_poly::nbar_ent_p2(const double ent2) const {
883  return puis(2*(exp(ent2) - m_2)/(gam2*kap2), 1./gam2m1) ;
884 }
885 
886 // Energy density from baryonic densities
887 //---------------------------------------
888 
889 double Eos_bf_poly::ener_nbar_p(const double nbar1, const double nbar2,
890  const double delta2) const {
891 
892  if (( nbar1 > double(0) ) || ( nbar2 > double(0))) {
893 
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 ;
897 
898  double resu = 0.5*kap1*pow(n1, gam1) + 0.5*kap2*pow(n2,gam2)
899  + kap3*pow(n1,gam3)*pow(n2,gam4) + m_1*n1 + m_2*n2
900  + x2*beta*pow(n1,gam5)*pow(n2,gam6) ;
901  return resu ;
902  }
903  else return 0 ;
904 }
905 
906 // Pressure from baryonic densities
907 //---------------------------------
908 
909 double Eos_bf_poly::press_nbar_p(const double nbar1, const double nbar2,
910  const double delta2) const {
911 
912  if (( nbar1 > double(0) ) || ( nbar2 > double(0))) {
913 
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 ;
917 
918  double resu = 0.5*gam1m1*kap1*pow(n1,gam1) + 0.5*gam2m1*kap2*pow(n2,gam2)
919  + gam34m1*kap3*pow(n1,gam3)*pow(n2,gam4) +
920  x2*gam56m1*beta*pow(n1,gam5)*pow(n2,gam6) ;
921 
922  return resu ;
923  }
924  else return 0 ;
925 }
926 
927 // Derivatives of energy
928 //----------------------
929 
930 double Eos_bf_poly::get_K11(const double n1, const double n2, const
931  double delta2) const
932 {
933  double xx ;
934  if (n1 <= 0.) {
935  xx = 0. ;
936  }
937  else {
938  xx = 0.5*gam1*kap1 * pow(n1,gam1 - 2) + m_1/n1 +
939  gam3*kap3 * pow(n1,gam3 - 2) * pow(n2,gam4) +
940  (delta2*(gam5 + 2) - 2)*beta * pow(n1,gam5 - 2)*pow(n2, gam6) ;
941  }
942  return xx ;
943 }
944 
945 double Eos_bf_poly::get_K22(const double n1, const double n2, const
946  double delta2) const
947 {
948  double xx ;
949  if (n2 <= 0.) {
950  xx = 0. ;
951  }
952  else {
953  xx = 0.5*gam2*kap2 * pow(n2,gam2 - 2) + m_2/n2 +
954  gam4*kap3 * pow(n2,gam4 - 2) * pow(n1,gam3) +
955  (delta2*(gam6 + 2) - 2)*beta * pow(n1, gam5) * pow(n2,gam6 - 2) ;
956  }
957  return xx ;
958 }
959 
960 double Eos_bf_poly::get_K12(const double n1, const double n2, const
961  double delta2) const
962 {
963  double xx ;
964  if ((n1 <= 0.) || (n2 <= 0.)) { xx = 0.; }
965  else {
966  double gamma_delta3 = pow(1-delta2,-1.5) ;
967  xx = 2*beta*pow(n1,gam5-1)*pow(n2,gam6-1) / gamma_delta3 ;
968  }
969  return xx ;
970 }
971 
972 // Conversion functions
973 // ---------------------
974 
976 
977  Eos_poly* eos_simple = new Eos_poly(gam1, kap1, m_1) ;
978  return eos_simple ;
979 }
980 /***************************************************************************
981  *
982  * Non-class functions
983  *
984  ***************************************************************************/
985 
986 // New "pow"
987 //-----------
988 
989  double puis(double x, double p) {
990  assert(p>=0.) ;
991  if (p==0.) return (x>=0 ? 1 : -1) ;
992  //if (p==0.) return 1 ;
993  if (x<0.) return (-pow(-x,p)) ;
994  else return pow(x,p) ;
995 }
996 
997 // Auxilliary functions for nbar_ent_p
998 //------------------------------------
999 
1000 double enthal1(const double x, const Param& parent) {
1001  assert(parent.get_n_double() == 7) ;
1002 
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) ;
1006 
1007 }
1008 
1009 double enthal23(const double x, const Param& parent) {
1010  assert(parent.get_n_double() == 10) ;
1011 
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) ;
1022 
1023 }
1024 
1025 double enthal(const double x, const Param& parent) {
1026  assert(parent.get_n_double_mod() == 7) ;
1027 
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) ;
1035 
1036  return (cc1*puis(x,alp1) + cc2*puis(x,alp2) + cc3*puis(x,alp3) - cc4) ;
1037 
1038 }
1039 
1040 }
1041 
Analytic equation of state for two fluids (relativistic case).
Definition: eos_bifluid.h:814
double gam5
Adiabatic indexes , see Eq.~eeosbfpolye}.
Definition: eos_bifluid.h:833
virtual bool operator!=(const Eos_bifluid &) const
Comparison operator (difference)
Definition: eos_bf_poly.C:441
double kap1
Pressure coefficient , see Eq.
Definition: eos_bifluid.h:843
void determine_type()
Determines the type of the analytical EOS (see typeos )
Definition: eos_bf_poly.C:334
double gam1
Adiabatic indexes , see Eq.~eeosbfpolye}.
Definition: eos_bifluid.h:821
double kap2
Pressure coefficient , see Eq.
Definition: eos_bifluid.h:850
void operator=(const Eos_bf_poly &)
Assignment to another Eos_bf_poly.
Definition: eos_bf_poly.C:294
double gam4
Adiabatic indexes , see Eq.~eeosbfpolye}.
Definition: eos_bifluid.h:830
double gam3
Adiabatic indexes , see Eq.~eeosbfpolye}.
Definition: eos_bifluid.h:827
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) .
Definition: eos_bf_poly.C:930
virtual double get_K12(const double n1, const double n2, const double delta2) const
Computes the derivative of the energy with respect to .
Definition: eos_bf_poly.C:960
virtual bool operator==(const Eos_bifluid &) const
Comparison operator (egality)
Definition: eos_bf_poly.C:386
virtual ostream & operator>>(ostream &) const
Operator >>
Definition: eos_bf_poly.C:472
virtual Eos * trans2Eos() const
Makes a translation from Eos_bifluid to Eos .
Definition: eos_bf_poly.C:975
virtual double get_K22(const double n1, const double n2, const double delta2) const
Computes the derivative of the energy/(baryonic density 2) .
Definition: eos_bf_poly.C:945
double gam6
Adiabatic indexes , see Eq.~eeosbfpolye}.
Definition: eos_bifluid.h:836
double gam2
Adiabatic indexes , see Eq.~eeosbfpolye}.
Definition: eos_bifluid.h:824
double precis
contains the precision required in zerosec_b
Definition: eos_bifluid.h:901
virtual int identify() const
Returns a number to identify the sub-classe of Eos_bifluid the object belongs to.
Definition: eos_bf_file.C:95
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.
Definition: eos_bf_poly.C:909
virtual void sauve(FILE *) const
Save in a file.
Definition: eos_bf_poly.C:452
void set_auxiliary()
Computes the auxiliary quantities gam1m1 , gam2m1 and gam3m1.
Definition: eos_bf_poly.C:320
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.
Definition: eos_bf_poly.C:889
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.
Definition: eos_bf_poly.C:507
int typeos
The bi-fluid analytical EOS type:
Definition: eos_bifluid.h:893
double beta
Coefficient , see Eq.
Definition: eos_bifluid.h:864
double kap3
Pressure coefficient , see Eq.
Definition: eos_bifluid.h:857
virtual double nbar_ent_p2(const double ent2) const
Computes baryon density out of the log-enthalpy assuming that only fluid 2 is present.
Definition: eos_bf_poly.C:882
virtual double nbar_ent_p1(const double ent1) const
Computes baryon density out of the log-enthalpy asuming that only fluid 1 is present.
Definition: eos_bf_poly.C:878
double ecart
contains the precision required in the relaxation nbar_ent_p
Definition: eos_bifluid.h:904
Eos_bf_poly(double kappa1, double kappa2, double kappa3, double beta)
Standard constructor.
Definition: eos_bf_poly.C:157
double relax
Parameters needed for some inversions of the EOS.
Definition: eos_bifluid.h:899
virtual ~Eos_bf_poly()
Destructor.
Definition: eos_bf_poly.C:285
2-fluids equation of state base class.
Definition: eos_bifluid.h:173
virtual void sauve(FILE *) const
Save in a file.
Definition: eos_bifluid.C:245
double m_1
Individual particle mass [unit: ].
Definition: eos_bifluid.h:184
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.
Definition: eos_bifluid.C:232
double m_2
Individual particle mass [unit: ].
Definition: eos_bifluid.h:189
Polytropic equation of state (relativistic case).
Definition: eos.h:757
Equation of state base class.
Definition: eos.h:190
Parameter storage.
Definition: param.h:125
void add_double(const double &x, int position=0)
Adds the the address of a new double to the list.
Definition: param.C:315
void add_double_mod(double &x, int position=0)
Adds the address of a new modifiable double to the list.
Definition: param.C:453
double & get_double_mod(int position=0) const
Returns the reference of a stored modifiable double .
Definition: param.C:498
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:220
Cmp exp(const Cmp &)
Exponential.
Definition: cmp_math.C:270
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
Cmp log(const Cmp &)
Neperian logarithm.
Definition: cmp_math.C:296
int fread_be(int *aa, int size, int nb, FILE *fich)
Reads integer(s) from a binary file according to the big endian convention.
Definition: fread_be.C:69
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.
Definition: zerosec_borne.C:68
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.
Definition: fwrite_be.C:70
Lorene prototypes.
Definition: app_hor.h:64