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