LORENE
eos_bf_poly_newt.C
1 /*
2  * Methods of the class Eos_bf_poly_newt.
3  *
4  * (see file eos_bifluid.h for documentation).
5  */
6 
7 /*
8  * Copyright (c) 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_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 $" ;
30 
31 /*
32  * $Id: eos_bf_poly_newt.C,v 1.15 2014/10/13 08:52:52 j_novak Exp $
33  * $Log: eos_bf_poly_newt.C,v $
34  * Revision 1.15 2014/10/13 08:52:52 j_novak
35  * Lorene classes and functions now belong to the namespace Lorene.
36  *
37  * Revision 1.14 2014/10/06 15:13:06 j_novak
38  * Modified #include directives to use c++ syntax.
39  *
40  * Revision 1.13 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.12 2003/12/17 23:12:32 r_prix
44  * replaced use of C++ <string> by standard ANSI char* to be backwards compatible
45  * with broken compilers like MIPSpro Compiler 7.2 on SGI Origin200. ;-)
46  *
47  * Revision 1.11 2003/12/05 15:09:47 r_prix
48  * adapted Eos_bifluid class and subclasses to use read_variable() for
49  * (formatted) file-reading.
50  *
51  * Revision 1.10 2003/12/04 14:22:33 r_prix
52  * removed enthalpy restrictions in eos-inversion
53  *
54  * Revision 1.9 2003/11/18 18:28:38 r_prix
55  * moved particle-masses m_1, m_2 of the two fluids into class eos_bifluid (from eos_bf_poly)
56  *
57  * Revision 1.8 2003/02/07 10:47:43 j_novak
58  * The possibility of having gamma5 xor gamma6 =0 has been introduced for
59  * tests. The corresponding parameter files have been added.
60  *
61  * Revision 1.7 2003/02/06 16:05:56 j_novak
62  * Corrected an error in the inversion of the EOS for typeos =1,2 and 3.
63  * Added new parameter files for sfstar.
64  *
65  * Revision 1.6 2002/10/16 14:36:34 j_novak
66  * Reorganization of #include instructions of standard C++, in order to
67  * use experimental version 3 of gcc.
68  *
69  * Revision 1.5 2002/05/31 16:13:36 j_novak
70  * better inversion for eos_bifluid
71  *
72  * Revision 1.4 2002/05/02 15:16:22 j_novak
73  * Added functions for more general bi-fluid EOS
74  *
75  * Revision 1.3 2002/04/05 09:09:36 j_novak
76  * The inversion of the EOS for 2-fluids polytrope has been modified.
77  * Some errors in the determination of the surface were corrected.
78  *
79  * Revision 1.2 2002/01/16 15:03:28 j_novak
80  * *** empty log message ***
81  *
82  * Revision 1.1 2002/01/11 14:09:34 j_novak
83  * Added newtonian version for 2-fluid stars
84  *
85  *
86  * $Header: /cvsroot/Lorene/C++/Source/Eos/eos_bf_poly_newt.C,v 1.15 2014/10/13 08:52:52 j_novak Exp $
87  *
88  */
89 
90 // Headers C
91 #include <cstdlib>
92 #include <cmath>
93 
94 // Headers Lorene
95 #include "eos_bifluid.h"
96 #include "eos.h"
97 #include "cmp.h"
98 #include "utilitaires.h"
99 
100 namespace Lorene {
101 
102 //****************************************************************************
103 // local prototypes
104 double puis(double, double) ;
105 double enthal1(const double x, const Param& parent) ;
106 double enthal23(const double x, const Param& parent) ;
107 double enthal(const double x, const Param& parent) ;
108 //****************************************************************************
109 
110  //--------------//
111  // Constructors //
112  //--------------//
113 
114 // Standard constructor with gam1 = gam2 = 2,
115 // gam3 = gam4 = gam5 = gam6 = 1, m_1 = 1 and m_2 =1
116 // -------------------------------------------------
117 Eos_bf_poly_newt::Eos_bf_poly_newt(double kappa1, double kappa2, double kappa3,
118  double bet) :
119  Eos_bf_poly(kappa1, kappa2, kappa3, bet) {
120  name = "bi-fluid polytropic non-relativistic EOS" ;
121 }
122 
123 // Standard constructor with everything specified
124 // -----------------------------------------------
125 Eos_bf_poly_newt::Eos_bf_poly_newt(double gamma1, double gamma2, double gamma3,
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,
132  l_ecart) {
133  name = "bi-fluid polytropic non-relativistic EOS" ;
134 }
135 
136 // Copy constructor
137 // ----------------
139  Eos_bf_poly(eosi) {}
140 
141 
142 // Constructor from binary file
143 // ----------------------------
145  Eos_bf_poly(fich) {}
146 
147 // Constructor from a formatted file
148 // ---------------------------------
150  Eos_bf_poly(fname) {}
151 
152  //--------------//
153  // Destructor //
154  //--------------//
155 
157 
158  // does nothing
159 
160 }
161  //--------------//
162  // Assignment //
163  //--------------//
164 
166 
167  Eos_bf_poly::operator=(eosi) ;
168 
169 }
170 
171  //------------------------//
172  // Comparison operators //
173  //------------------------//
174 
175 
176 bool Eos_bf_poly_newt::operator==(const Eos_bifluid& eos_i) const {
177 
178  bool resu = true ;
179 
180  if ( eos_i.identify() != identify() ) {
181  cout << "The second EOS is not of type Eos_bf_poly_newt !" << endl ;
182  resu = false ;
183  }
184  else{
185 
186  const Eos_bf_poly_newt& eos =
187  dynamic_cast<const Eos_bf_poly_newt&>( eos_i ) ;
188 
189  if ((eos.gam1 != gam1)||(eos.gam2 != gam2)||(eos.gam3 != gam3)
190  ||(eos.gam4 != gam4)||(eos.gam5 != gam5)||(eos.gam6 != gam6)) {
191  cout
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 ;
199  resu = false ;
200  }
201 
202  if ((eos.kap1 != kap1)||(eos.kap2 != kap2)|| (eos.kap3 != kap3)){
203  cout
204  << "The two Eos_bf_poly_newt have different kappas : " << kap1 << " <-> "
205  << eos.kap1 << ", " << kap2 << " <-> "
206  << eos.kap2 << ", " << kap3 << " <-> "
207  << eos.kap3 << endl ;
208  resu = false ;
209  }
210 
211  if (eos.beta != beta) {
212  cout
213  << "The two Eos_bf_poly_newt have different betas : " << beta << " <-> "
214  << eos.beta << endl ;
215  resu = false ;
216  }
217 
218  if ((eos.m_1 != m_1)||(eos.m_2 != m_2)) {
219  cout
220  << "The two Eos_bf_poly_newt have different masses : " << m_1 << " <-> "
221  << eos.m_1 << ", " << m_2 << " <-> "
222  << eos.m_2 << endl ;
223  resu = false ;
224  }
225 
226  }
227 
228  return resu ;
229 
230 }
231 
232 bool Eos_bf_poly_newt::operator!=(const Eos_bifluid& eos_i) const {
233 
234  return !(operator==(eos_i)) ;
235 
236 }
237 
238 
239  //------------//
240  // Outputs //
241  //------------//
242 
243 void Eos_bf_poly_newt::sauve(FILE* fich) const {
244 
245  Eos_bf_poly::sauve(fich) ;
246 }
247 
248 ostream& Eos_bf_poly_newt::operator>>(ostream & ost) const {
249 
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 ;
271 
272  return ost ;
273 
274 }
275 
276 
277  //------------------------------//
278  // Computational routines //
279  //------------------------------//
280 
281 // Baryon densities from enthalpies
282 //---------------------------------
283 // RETURN: bool true = if in a one-fluid region, false if two-fluids
284 
285 bool Eos_bf_poly_newt::nbar_ent_p(const double ent1, const double ent2,
286  const double delta2, double& nbar1,
287  double& nbar2) const {
288 
289  //RP: I think this is wrong!!
290  // bool one_fluid = ((ent1<=0.)||(ent2<=0.)) ;
291 
292  bool one_fluid = false;
293 
294  if (!one_fluid) {
295  switch (typeos) {
296  case 5: // same as typeos==0, just with slow-rot-style inversion
297  case 0: {//gamma1=gamma2=2 all others = 1 easy case of a linear system
298  double kpd = kap3+beta*delta2 ;
299  double determ = kap1*kap2 - kpd*kpd ;
300 
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.)) ;
304  break ;
305  }
306  case 1: { //gamma1 or gamma 2 not= 2; all others = 1
307  double b1 = ent1*m_1 ;
308  double b2 = ent2*m_2 ;
309  double denom = kap3 + beta*delta2 ;
310  if (fabs(denom) < 1.e-15) {
311  nbar1 = puis(2*b1/(gam1*kap1), 1./gam1m1) ;
312  nbar2 = puis(2*b2/(gam2*kap2), 1./gam2m1) ;
313  one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
314  }
315  else {
316  double coef0 = 0.5*gam2*kap2/pow(denom, gam2m1) ;
317  double coef1 = 0.5*kap1*gam1 ;
318  Param parent ;
319  parent.add_double(coef0, 0) ;
320  parent.add_double(b1, 1) ;
321  parent.add_double(coef1, 2) ;
322  parent.add_double(gam1m1,3) ;
323  parent.add_double(gam2m1,4) ;
324  parent.add_double(denom, 5) ;
325  parent.add_double(b2, 6) ;
326 
327  double xmin, xmax ;
328  double f0 = enthal1(0.,parent) ;
329  if (fabs(f0)<1.e-15) {
330  nbar1 = 0. ;}
331  else {
332  double cas = (gam1m1*gam2m1 - 1.)*f0;
333  assert (fabs(cas) > 1.e-15) ;
334  xmin = 0. ;
335  xmax = cas/fabs(cas) ;
336  do {
337  xmax *= 1.1 ;
338  if (fabs(xmax) > 1.e10) {
339  cout << "Problem in Eos_bf_poly::nbar_ent_p!" << endl ;
340  cout << f0 << ", " << cas << endl ; // to be removed!!
341  cout << "typeos = 1" << endl ;
342  abort() ;
343  }
344  } while (enthal1(xmax,parent)*f0 > 0.) ;
345  double l_precis = 1.e-12 ;
346  int nitermax = 400 ;
347  int niter = 0 ;
348  nbar1 = zerosec_b(enthal1, parent, xmin, xmax, l_precis,
349  nitermax, niter) ;
350  }
351  nbar2 = (b1 - coef1*puis(nbar1, gam1m1))/denom ;
352  double resu1 = coef1*puis(nbar1,gam1m1) + denom*nbar2 ;
353  double resu2 = 0.5*gam2*kap2*puis(nbar2,gam2m1) + denom*nbar1 ;
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)) ;
357  }
358  break ;
359  }
360  case 2: { // gamma3 = gamma5 = 1 at least
361  double b1 = ent1*m_1 ;
362  double b2 = ent2*m_2 ;
363  double denom = kap3 + beta*delta2 ;
364  if (fabs(denom) < 1.e-15) {
365  nbar1 = puis(2*b1/(gam1*kap1), 1./gam1m1) ;
366  nbar2 = puis(2*b2/(gam2*kap2), 1./gam2m1) ;
367  one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
368  }
369  else {
370  double coef0 = beta*delta2 ;
371  double coef1 = 0.5*kap1*gam1 ;
372  double coef2 = 0.5*kap2*gam2 ;
373  double coef3 = 1./gam1m1 ;
374  Param parent ;
375  parent.add_double(b1, 0) ;
376  parent.add_double(kap3, 1) ;
377  parent.add_double(gam4, 2) ;
378  parent.add_double(coef0, 3) ;
379  parent.add_double(gam6,4) ;
380  parent.add_double(coef1, 5) ;
381  parent.add_double(coef3, 6) ;
382  parent.add_double(coef2, 7) ;
383  parent.add_double(gam2m1, 8) ;
384  parent.add_double(b2, 9) ;
385 
386  double xmin, xmax ;
387  double f0 = enthal23(0.,parent) ;
388  if (fabs(f0)<1.e-15) {
389  nbar2 = 0. ;}
390  else {
391  double pmax = (fabs(kap3) < 1.e-15 ? 0. : gam4*(gam4-1) ) ;
392  double ptemp = (fabs(kap3*coef0) < 1.e-15 ? 0.
393  : gam6*(gam4-1) ) ;
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) ;
399  double cas = (pmax - gam1m1*gam2m1)*f0;
400  // cout << "Pmax, cas: " << pmax << ", " << cas << endl ;
401  assert (fabs(cas) > 1.e-15) ;
402  xmin = 0. ;
403  xmax = cas/fabs(cas) ;
404  do {
405  xmax *= 1.1 ;
406  if (fabs(xmax) > 1.e10) {
407  cout << "Problem in Eos_bf_poly::nbar_ent_p!" << endl ;
408  cout << "typeos = 2" << endl ;
409  abort() ;
410  }
411  } while (enthal23(xmax,parent)*f0 > 0.) ;
412  double l_precis = 1.e-12 ;
413  int nitermax = 400 ;
414  int niter = 0 ;
415  nbar2 = zerosec_b(enthal23, parent, xmin, xmax, l_precis,
416  nitermax, niter) ;
417  }
418  nbar1 = (b1 - kap3*puis(nbar2,gam4) - coef0*puis(nbar2,gam6))
419  / coef1 ;
420  nbar1 = puis(nbar1,coef3) ;
421  double resu1 = coef1*puis(nbar1,gam1m1) + kap3*puis(nbar2,gam4)
422  + coef0*puis(nbar2, gam6) ;
423  double resu2 = coef2*puis(nbar2,gam2m1)
424  + gam4*kap3*puis(nbar2, gam4-1)*nbar1
425  + gam6*coef0*puis(nbar2, gam6-1)*nbar1 ;
426  double erreur = fabs((resu1/m_1-ent1)/ent1) +
427  fabs((resu2/m_2-ent2)/ent2) ;
428  //cout << "Erreur d'inversion: " << erreur << endl ;
429  one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)||(erreur > 1.e-4)) ;
430  }
431  break ;
432  }
433  case 3: { //gamma4 = gamm6 = 1 at least
434  double b1 = ent1*m_1 ;
435  double b2 = ent2*m_2 ;
436  double denom = kap3 + beta*delta2 ;
437  if (fabs(denom) < 1.e-15) {
438  nbar1 = puis(2*b1/(gam1*kap1), 1./gam1m1) ;
439  nbar2 = puis(2*b2/(gam2*kap2), 1./gam2m1) ;
440  one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
441  }
442  else {
443  double coef0 = beta*delta2 ;
444  double coef1 = 0.5*kap1*gam1 ;
445  double coef2 = 0.5*kap2*gam2 ;
446  double coef3 = 1./gam2m1 ;
447  Param parent ;
448  parent.add_double(b2, 0) ;
449  parent.add_double(kap3, 1) ;
450  parent.add_double(gam3, 2) ;
451  parent.add_double(coef0, 3) ;
452  parent.add_double(gam5, 4) ;
453  parent.add_double(coef2, 5) ;
454  parent.add_double(coef3, 6) ;
455  parent.add_double(coef1, 7) ;
456  parent.add_double(gam1m1, 8) ;
457  parent.add_double(b1, 9) ;
458 
459  double xmin, xmax ;
460  double f0 = enthal23(0.,parent) ;
461  if (fabs(f0)<1.e-15) {
462  nbar1 = 0. ;}
463  else {
464  double pmax = (fabs(kap3) < 1.e-15 ? 0. : gam3*(gam3-1) ) ;
465  double ptemp = (fabs(kap3*coef0) < 1.e-15 ? 0.
466  : gam5*(gam3-1) ) ;
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) ;
472  double cas = (pmax - gam1m1*gam2m1)*f0;
473  // cout << "Pmax, cas: " << pmax << ", " << cas << endl ;
474  assert (fabs(cas) > 1.e-15) ;
475  xmin = 0. ;
476  xmax = cas/fabs(cas) ;
477  do {
478  xmax *= 1.1 ;
479  if (fabs(xmax) > 1.e10) {
480  cout << "Problem in Eos_bf_poly::nbar_ent_p!" << endl ;
481  cout << "typeos = 3" << endl ;
482  abort() ;
483  }
484  } while (enthal23(xmax,parent)*f0 > 0.) ;
485  double l_precis = 1.e-12 ;
486  int nitermax = 400 ;
487  int niter = 0 ;
488  nbar1 = zerosec_b(enthal23, parent, xmin, xmax, l_precis,
489  nitermax, niter) ;
490  }
491  nbar2 = (b2 - kap3*puis(nbar1,gam3) - coef0*puis(nbar1,gam5))
492  / coef2 ;
493  nbar2 = puis(nbar2,coef3) ;
494  double resu1 = coef1*puis(nbar1,gam1m1)
495  + gam3*kap3*puis(nbar1,gam3-1)*nbar2
496  + coef0*gam5*puis(nbar1, gam5-1)*nbar2 ;
497  double resu2 = coef2*puis(nbar2,gam2m1)
498  + kap3*puis(nbar1, gam3) + coef0*puis(nbar1, gam5);
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)) ;
502  }
503  break ;
504  }
505  case 4:{ // most general case
506  double b1 = ent1*m_1 ;
507  double b2 = ent2*m_2 ;
508  double denom = kap3 + beta*delta2 ;
509  if (fabs(denom) < 1.e-15) {
510  nbar1 = puis(2*b1/(gam1*kap1), 1./gam1m1) ;
511  nbar2 = puis(2*b2/(gam2*kap2), 1./gam2m1) ;
512  one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
513  }
514  else {
515  int nitermax = 200 ;
516  int niter = 0 ;
517  int nmermax = 800 ;
518 
519  double a11 = 0.5*gam1*kap1 ;
520  double a13 = gam3*kap3 ;
521  double a14 = beta*gam5*delta2 ;
522 
523  double a22 = 0.5*gam2*kap2 ;
524  double a23 = gam4*kap3 ;
525  double a24 = beta*gam6*delta2 ;
526 
527  double n1l, n2l, n1s, n2s ;
528 
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)) ;
534 
535  double n1m = (n1l<0. ? n1s : sqrt(n1l*n1s)) ;
536  double n2m = (n2l<0. ? n2s : sqrt(n2l*n2s)) ;
537 
538  Param parf1 ;
539  Param parf2 ;
540  double c1, c2, c3, c4, c5, c6, c7 ;
541  c1 = gam1m1 ;
542  c2 = gam3 - 1. ;
543  c3 = gam5 - 1. ;
544  c4 = a11 ;
545  c5 = a13*puis(n2m,gam4) ;
546  c6 = a14*puis(n2m,gam6) ;
547  c7 = b1 ;
548  parf1.add_double_mod(c1,0) ;
549  parf1.add_double_mod(c2,1) ;
550  parf1.add_double_mod(c3,2) ;
551  parf1.add_double_mod(c4,3) ;
552  parf1.add_double_mod(c5,4) ;
553  parf1.add_double_mod(c6,5) ;
554  parf1.add_double_mod(c7,6) ;
555 
556  double d1, d2, d3, d4, d5, d6, d7 ;
557  d1 = gam2m1 ;
558  d2 = gam4 - 1. ;
559  d3 = gam6 - 1. ;
560  d4 = a22 ;
561  d5 = a23*puis(n1m,gam3) ;
562  d6 = a24*puis(n1m,gam5) ;
563  d7 = b2 ;
564  parf2.add_double_mod(d1,0) ;
565  parf2.add_double_mod(d2,1) ;
566  parf2.add_double_mod(d3,2) ;
567  parf2.add_double_mod(d4,3) ;
568  parf2.add_double_mod(d5,4) ;
569  parf2.add_double_mod(d6,5) ;
570  parf2.add_double_mod(d7,6) ;
571 
572  double xmin = -3*(n1s>n2s ? n1s : n2s) ;
573  double xmax = 3*(n1s>n2s ? n1s : n2s) ;
574 
575  double n1 = n1m ;
576  double n2 = n2m ;
577  bool sortie = true ;
578  int mer = 0 ;
579 
580  //cout << "Initial guess: " << n1 << ", " << n2 << endl ;
581  n1s *= 0.1 ;
582  n2s *= 0.1 ;
583  do {
584  //cout << "n1, n2: " << n1 << ", " << n2 << endl ;
585  n1 = zerosec_b(&enthal, parf1, xmin, xmax, precis, nitermax, niter) ;
586  n2 = zerosec_b(&enthal, parf2, xmin, xmax, precis, nitermax, niter) ;
587 
588  sortie = (fabs((n1m-n1)/(n1m+n1)) + fabs((n2m-n2)/(n2m+n2)) > ecart) ;
589  n1m = relax*n1 + (1.-relax)*n1m ;
590  n2m = relax*n2 + (1.-relax)*n2m ;
591  if (n2m>-n2s) {
592  parf1.get_double_mod(4) = a13*puis(n2m,gam4) ;
593  parf1.get_double_mod(5) = a14*puis(n2m,gam6) ;
594  }
595  else {
596  parf1.get_double_mod(4) = a13*puis(-n2s,gam4) ;
597  parf1.get_double_mod(5) = a14*puis(-n2s,gam6) ;
598  }
599 
600  if (n1m>-n1s) {
601  parf2.get_double_mod(4) = a23*puis(n1m,gam3) ;
602  parf2.get_double_mod(5) = a24*puis(n1m,gam5) ;
603  }
604  else {
605  parf2.get_double_mod(4) = a23*puis(-n1s,gam3) ;
606  parf2.get_double_mod(5) = a24*puis(-n1s,gam5) ;
607  }
608 
609  mer++ ;
610  } while (sortie&&(mer<nmermax)) ;
611  nbar1 = n1m ;
612  nbar2 = n2m ;
613 
614 // double resu1 = a11*puis(n1,gam1m1) + a13*puis(n1,gam3-1.)*puis(n2,gam4)
615 // +a14*puis(n1,gam5-1.)*puis(n2,gam6) ;
616 // double resu2 = a22*puis(n2,gam2m1) + a23*puis(n1,gam3)*puis(n2,gam4-1.)
617 // +a24*puis(n1,gam5)*puis(n2,gam6-1.) ;
618  //cout << "Nbre d'iterations: " << mer << endl ;
619  //cout << "Resus: " << n1m << ", " << n2m << endl ;
620  //cout << "Verification: " << log(fabs(1+resu1)) << ", "
621  // << log(fabs(1+resu2)) << endl ;
622  // cout << "Erreur: " << fabs(enthal(n1, parf1)/b1) +
623  // fabs(enthal(n2, parf2)/b2) << endl ;
624  //cout << "Erreur: " << fabs((log(fabs(1+resu1))-ent1)/ent1) +
625  //fabs((log(fabs(1+resu2))-ent2)/ent2) << endl ;
626  }
627  break ;
628  }
629  case -1: {
630  double determ = kap1*kap2 - kap3*kap3 ;
631 
632  nbar1 = (kap2*(ent1*m_1 - beta*delta2) - kap3*ent2*m_2) / determ ;
633  nbar2 = (kap1*ent2*m_2 - kap3*(ent1*m_1 - beta*delta2)) / determ ;
634  one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
635  break ;
636  }
637  case -2: {
638  double determ = kap1*kap2 - kap3*kap3 ;
639 
640  nbar1 = (kap2*ent1*m_1 - kap3*(ent2*m_2 - beta*delta2)) / determ ;
641  nbar2 = (kap1*(ent2*m_2 - beta*delta2) - kap3*ent1*m_1) / determ ;
642  one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
643  break ;
644  }
645  }
646  }
647 
648  return one_fluid ;
649 }
650 // One fluid sub-EOSs
651 //-------------------
652 
653 double Eos_bf_poly_newt::nbar_ent_p1(const double ent1) const {
654  return puis(2*ent1*m_1/(gam1*kap1), 1./gam1m1) ;
655 }
656 
657 double Eos_bf_poly_newt::nbar_ent_p2(const double ent2) const {
658  return puis(2*ent2*m_2/(gam2*kap2), 1./gam2m1) ;
659 }
660 
661 // Energy density from baryonic densities
662 //---------------------------------------
663 
664 double Eos_bf_poly_newt::ener_nbar_p(const double nbar1, const double nbar2,
665  const double delta2) const {
666 
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 ;
670 
671  double resu = 0.5*kap1*pow(n1, gam1) + 0.5*kap2*pow(n2,gam2)
672  + kap3*pow(n1,gam3)*pow(n2,gam4)
673  + x2*beta*pow(n1,gam5)*pow(n2,gam6) ;
674 
675  return resu ;
676 
677 }
678 
679 // Pressure from baryonic densities
680 //---------------------------------
681 
682 double Eos_bf_poly_newt::press_nbar_p(const double nbar1, const double nbar2,
683  const double delta2) const {
684 
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 ;
688 
689  double resu = 0.5*gam1m1*kap1*pow(n1,gam1) + 0.5*gam2m1*kap2*pow(n2,gam2)
690  + gam34m1*kap3*pow(n1,gam3)*pow(n2,gam4) +
691  x2*gam56m1*beta*pow(n1,gam5)*pow(n2,gam6) ;
692 
693  return resu ;
694 
695 }
696 
697 // Derivatives of energy
698 //----------------------
699 
700 double Eos_bf_poly_newt::get_K11(const double n1, const double n2, const
701  double) const
702 {
703  double xx ;
704  if (n1 <= 0.) xx = 0. ;
705  else xx = m_1/n1 -2*beta*pow(n1,gam5-2)*pow(n2,gam6) ;
706 
707  return xx ;
708 }
709 
710 double Eos_bf_poly_newt::get_K22(const double n1, const double n2, const
711  double ) const
712 {
713  double xx ;
714  if (n2 <= 0.) xx = 0. ;
715  else xx = m_2/n2 - 2*beta*pow(n1,gam5)*pow(n2,gam6-2) ;
716 
717  return xx ;
718 }
719 
720 double Eos_bf_poly_newt::get_K12(const double n1, const double n2, const
721  double) const
722 {
723  double xx ;
724  if ((n1 <= 0.) || (n2 <= 0.)) xx = 0.;
725  else xx = 2*beta*pow(n1,gam5-1)*pow(n2,gam6-1) ;
726 
727  return xx ;
728 }
729 
730 // Conversion functions
731 // ---------------------
732 
734 
735  Eos_poly_newt* eos_simple = new Eos_poly_newt(gam1, kap1) ;
736  return eos_simple ;
737 
738 }
739 
740 }
Analytic equation of state for two fluids (Newtonian case).
Definition: eos_bifluid.h:1258
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.
Definition: eos_bf_file.C:97
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).
Definition: eos_bifluid.h:814
double gam5
Adiabatic indexes , see Eq.~eeosbfpolye}.
Definition: eos_bifluid.h:833
double kap1
Pressure coefficient , see Eq.
Definition: eos_bifluid.h:843
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
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 void sauve(FILE *) const
Save in a file.
Definition: eos_bf_poly.C:452
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
double ecart
contains the precision required in the relaxation nbar_ent_p
Definition: eos_bifluid.h:904
double relax
Parameters needed for some inversions of the EOS.
Definition: eos_bifluid.h:899
2-fluids equation of state base class.
Definition: eos_bifluid.h:173
double m_1
Individual particle mass [unit: ].
Definition: eos_bifluid.h:184
string name
EOS name.
Definition: eos_bifluid.h:179
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: ].
Definition: eos_bifluid.h:189
Polytropic equation of state (Newtonian case).
Definition: eos.h:1044
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 pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
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
Lorene prototypes.
Definition: app_hor.h:64