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
29char 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
100namespace Lorene {
101
102//****************************************************************************
103// local prototypes
104double puis(double, double) ;
105double enthal1(const double x, const Param& parent) ;
106double enthal23(const double x, const Param& parent) ;
107double 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// -------------------------------------------------
117Eos_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// -----------------------------------------------
125Eos_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// ----------------
140
141
142// Constructor from binary file
143// ----------------------------
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
168
169}
170
171 //------------------------//
172 // Comparison operators //
173 //------------------------//
174
175
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
233
234 return !(operator==(eos_i)) ;
235
236}
237
238
239 //------------//
240 // Outputs //
241 //------------//
242
243void Eos_bf_poly_newt::sauve(FILE* fich) const {
244
245 Eos_bf_poly::sauve(fich) ;
246}
247
248ostream& 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
285bool 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
653double Eos_bf_poly_newt::nbar_ent_p1(const double ent1) const {
654 return puis(2*ent1*m_1/(gam1*kap1), 1./gam1m1) ;
655}
656
657double 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
664double 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
682double 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
700double 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
710double 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
720double 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).
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).
double gam5
Adiabatic indexes , see Eq.~eeosbfpolye}.
double kap1
Pressure coefficient , see Eq.
double gam1
Adiabatic indexes , see Eq.~eeosbfpolye}.
double kap2
Pressure coefficient , see Eq.
void operator=(const Eos_bf_poly &)
Assignment to another Eos_bf_poly.
double gam4
Adiabatic indexes , see Eq.~eeosbfpolye}.
double gam3
Adiabatic indexes , see Eq.~eeosbfpolye}.
double gam6
Adiabatic indexes , see Eq.~eeosbfpolye}.
double gam2
Adiabatic indexes , see Eq.~eeosbfpolye}.
double precis
contains the precision required in zerosec_b
virtual void sauve(FILE *) const
Save in a file.
int typeos
The bi-fluid analytical EOS type:
double beta
Coefficient , see Eq.
double kap3
Pressure coefficient , see Eq.
double ecart
contains the precision required in the relaxation nbar_ent_p
double relax
Parameters needed for some inversions of the EOS.
2-fluids equation of state base class.
double m_1
Individual particle mass [unit: ].
string name
EOS name.
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: ].
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.
Lorene prototypes.
Definition app_hor.h:64