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