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 }
virtual bool operator!=(const Eos_bifluid &) const
Comparison operator (difference)
double m_1
Individual particle mass [unit: ].
Definition: eos_bifluid.h:184
double & get_double_mod(int position=0) const
Returns the reference of a stored modifiable double .
Definition: param.C:498
double kap2
Pressure coefficient , see Eq.
Definition: eos_bifluid.h:850
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:220
Lorene prototypes.
Definition: app_hor.h:64
double relax
Parameters needed for some inversions of the EOS.
Definition: eos_bifluid.h:899
void operator=(const Eos_bf_poly_newt &)
Assignment to another Eos_bf_poly_newt.
string name
EOS name.
Definition: eos_bifluid.h:179
virtual bool operator==(const Eos_bifluid &) const
Comparison operator (egality)
Equation of state base class.
Definition: eos.h:190
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:294
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:814
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
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:904
void add_double_mod(double &x, int position=0)
Adds the address of a new modifiable double to the list.
Definition: param.C:453
2-fluids equation of state base class.
Definition: eos_bifluid.h:173
int typeos
The bi-fluid analytical EOS type:
Definition: eos_bifluid.h:893
double beta
Coefficient , see Eq.
Definition: eos_bifluid.h:864
double m_2
Individual particle mass [unit: ].
Definition: eos_bifluid.h:189
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:857
Parameter storage.
Definition: param.h:125
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 .
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
virtual void sauve(FILE *) const
Save in a file.
Polytropic equation of state (Newtonian case).
Definition: eos.h:1044
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:452
double precis
contains the precision required in zerosec_b
Definition: eos_bifluid.h:901
double gam5
Adiabatic indexes , see Eq.~eeosbfpolye}.
Definition: eos_bifluid.h:833
void add_double(const double &x, int position=0)
Adds the the address of a new double to the list.
Definition: param.C:315
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
double kap1
Pressure coefficient , see Eq.
Definition: eos_bifluid.h:843
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:821
double gam6
Adiabatic indexes , see Eq.~eeosbfpolye}.
Definition: eos_bifluid.h:836
double gam3
Adiabatic indexes , see Eq.~eeosbfpolye}.
Definition: eos_bifluid.h:827
double gam4
Adiabatic indexes , see Eq.~eeosbfpolye}.
Definition: eos_bifluid.h:830
double gam2
Adiabatic indexes , see Eq.~eeosbfpolye}.
Definition: eos_bifluid.h:824