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