LORENE
eos_poly.C
1 /*
2  * Methods of the class Eos_poly.
3  *
4  * (see file eos.h for documentation).
5  */
6 
7 /*
8  * Copyright (c) 2000-2001 Eric Gourgoulhon
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_poly.C,v 1.14 2023/03/13 15:48:46 j_novak Exp $
33  * $Log: eos_poly.C,v $
34  * Revision 1.14 2023/03/13 15:48:46 j_novak
35  * Use of function expm1.
36  *
37  * Revision 1.13 2023/01/11 15:30:00 j_novak
38  * Change the computation at low enthalpies for better accuracy.
39  *
40  * Revision 1.12 2023/01/10 15:46:11 j_novak
41  * Improvement computation of derivatives, including sound speed.
42  *
43  * Revision 1.11 2021/05/14 15:39:22 g_servignat
44  * Added sound speed computation from enthalpy to Eos class and tabulated+polytropic derived classes
45  *
46  * Revision 1.10 2016/12/05 16:17:51 j_novak
47  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
48  *
49  * Revision 1.9 2014/10/13 08:52:53 j_novak
50  * Lorene classes and functions now belong to the namespace Lorene.
51  *
52  * Revision 1.8 2014/10/06 15:13:06 j_novak
53  * Modified #include directives to use c++ syntax.
54  *
55  * Revision 1.7 2009/05/25 06:52:27 k_taniguchi
56  * Allowed the case of mu_0 != 1 for der_ener_ent_p and der_press_ent_p.
57  *
58  * Revision 1.6 2003/12/10 08:58:20 r_prix
59  * - added new Eos_bifluid paramter for eos-file: bool slow_rot_style
60  * to indicate if we want this particular kind of EOS-inversion (only works for
61  * the Newtonian 'analytic' polytropes) --> replaces former dirty hack with gamma1<0
62  *
63  * Revision 1.5 2002/10/16 14:36:35 j_novak
64  * Reorganization of #include instructions of standard C++, in order to
65  * use experimental version 3 of gcc.
66  *
67  * Revision 1.4 2002/04/11 13:28:40 e_gourgoulhon
68  * Added the parameter mu_0
69  *
70  * Revision 1.3 2002/04/09 14:32:15 e_gourgoulhon
71  * 1/ Added extra parameters in EOS computational functions (argument par)
72  * 2/ New class MEos for multi-domain EOS
73  *
74  * Revision 1.2 2001/12/04 21:27:53 e_gourgoulhon
75  *
76  * All writing/reading to a binary file are now performed according to
77  * the big endian convention, whatever the system is big endian or
78  * small endian, thanks to the functions fwrite_be and fread_be
79  *
80  * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
81  * LORENE
82  *
83  * Revision 2.9 2001/02/23 15:17:30 eric
84  * Methodes der_nbar_ent_p, der_ener_ent_p, der_press_ent_p :
85  * traitement du cas ent<1.e-13 par un DL
86  * continuite des quantites pour ent<=0.
87  *
88  * Revision 2.8 2001/02/07 09:50:30 eric
89  * Suppression de la fonction derent_ent_p.
90  * Ajout des fonctions donnant les derivees de l'EOS:
91  * der_nbar_ent_p
92  * der_ener_ent_p
93  * der_press_ent_p
94  *
95  * Revision 2.7 2000/06/20 08:34:46 eric
96  * Ajout des fonctions get_gam(), etc...
97  *
98  * Revision 2.6 2000/02/14 14:49:34 eric
99  * Modif affichage.
100  *
101  * Revision 2.5 2000/02/14 14:33:15 eric
102  * Ajout du constructeur par lecture de fichier formate.
103  *
104  * Revision 2.4 2000/01/21 16:05:47 eric
105  * Corrige erreur dans set_auxiliary: calcul de gam1.
106  *
107  * Revision 2.3 2000/01/21 15:18:45 eric
108  * Ajout des operateurs de comparaison == et !=
109  *
110  * Revision 2.2 2000/01/18 14:26:37 eric
111  * *** empty log message ***
112  *
113  * Revision 2.1 2000/01/18 13:47:17 eric
114  * Premiere version operationnelle
115  *
116  * Revision 2.0 2000/01/18 10:46:28 eric
117  * *** empty log message ***
118  *
119  *
120  * $Header: /cvsroot/Lorene/C++/Source/Eos/eos_poly.C,v 1.14 2023/03/13 15:48:46 j_novak Exp $
121  *
122  */
123 
124 // Headers C
125 #include <cstdlib>
126 #include <cstring>
127 #include <cmath>
128 
129 // Headers Lorene
130 #include "eos.h"
131 #include "cmp.h"
132 #include "utilitaires.h"
133 
134  //--------------//
135  // Constructors //
136  //--------------//
137 
138 // Standard constructor with m_0 = 1 and mu_0 = 1
139 // -----------------------------------------------
140 namespace Lorene {
141 Eos_poly::Eos_poly(double gam0, double kappa) :
142  Eos("Relativistic polytropic EOS"),
143  gam(gam0), kap(kappa), m_0(double(1)), mu_0(double(1)) {
144 
145  set_auxiliary() ;
146 
147 }
148 
149 // Standard constructor with mu_0 = 1
150 // ----------------------------------
151 Eos_poly::Eos_poly(double gam0, double kappa, double mass) :
152  Eos("Relativistic polytropic EOS"),
153  gam(gam0), kap(kappa), m_0(mass), mu_0(double(1)) {
154 
155  set_auxiliary() ;
156 
157 }
158 
159 // Standard constructor with mu_0 = 1
160 // ----------------------------------
161 Eos_poly::Eos_poly(double gam0, double kappa, double mass, double mu_zero) :
162  Eos("Relativistic polytropic EOS"),
163  gam(gam0), kap(kappa), m_0(mass), mu_0(mu_zero) {
164 
165  set_auxiliary() ;
166 
167 }
168 
169 // Copy constructor
170 // ----------------
172  Eos(eosi),
173  gam(eosi.gam), kap(eosi.kap), m_0(eosi.m_0), mu_0(eosi.mu_0) {
174 
175  set_auxiliary() ;
176 
177 }
178 
179 
180 // Constructor from binary file
181 // ----------------------------
182 Eos_poly::Eos_poly(FILE* fich) :
183  Eos(fich) {
184 
185  fread_be(&gam, sizeof(double), 1, fich) ;
186  fread_be(&kap, sizeof(double), 1, fich) ;
187  fread_be(&m_0, sizeof(double), 1, fich) ;
188 
189  if (m_0 < 0) { // to ensure compatibility with previous version (revision <= 1.2)
190  // of Eos_poly
191  m_0 = fabs( m_0 ) ;
192  fread_be(&mu_0, sizeof(double), 1, fich) ;
193  }
194  else {
195  mu_0 = double(1) ;
196  }
197 
198  set_auxiliary() ;
199 
200 }
201 
202 
203 // Constructor from a formatted file
204 // ---------------------------------
205 Eos_poly::Eos_poly(ifstream& fich) :
206  Eos(fich) {
207 
208  char blabla[80] ;
209 
210  fich >> gam ; fich.getline(blabla, 80) ;
211  fich >> kap ; fich.getline(blabla, 80) ;
212  fich >> m_0 ; fich.getline(blabla, 80) ;
213 
214  if (m_0 < 0) { // to ensure compatibility with previous version (revision <= 1.2)
215  // of Eos_poly
216  m_0 = fabs( m_0 ) ;
217  fich >> mu_0 ; fich.getline(blabla, 80) ;
218  }
219  else {
220  mu_0 = double(1) ;
221  }
222 
223  set_auxiliary() ;
224 
225 }
226  //--------------//
227  // Destructor //
228  //--------------//
229 
231 
232  // does nothing
233 
234 }
235  //--------------//
236  // Assignment //
237  //--------------//
238 
239 void Eos_poly::operator=(const Eos_poly& eosi) {
240 
241  set_name(eosi.name) ;
242 
243  gam = eosi.gam ;
244  kap = eosi.kap ;
245  m_0 = eosi.m_0 ;
246  mu_0 = eosi.mu_0 ;
247 
248  set_auxiliary() ;
249 
250 }
251 
252 
253  //-----------------------//
254  // Miscellaneous //
255  //-----------------------//
256 
258 
259  gam1 = gam - double(1) ;
260 
261  unsgam1 = double(1) / gam1 ;
262 
263  gam1sgamkap = m_0 * gam1 / (gam * kap) ;
264 
265  rel_mu_0 = mu_0 / m_0 ;
266 
267  ent_0 = log( rel_mu_0 ) ;
268 
269 }
270 
271 double Eos_poly::get_gam() const {
272  return gam ;
273 }
274 
275 double Eos_poly::get_kap() const {
276  return kap ;
277 }
278 
279 double Eos_poly::get_m_0() const {
280  return m_0 ;
281 }
282 
283 double Eos_poly::get_mu_0() const {
284  return mu_0 ;
285 }
286 
287 
288  //------------------------//
289  // Comparison operators //
290  //------------------------//
291 
292 
293 bool Eos_poly::operator==(const Eos& eos_i) const {
294 
295  bool resu = true ;
296 
297  if ( eos_i.identify() != identify() ) {
298  cout << "The second EOS is not of type Eos_poly !" << endl ;
299  resu = false ;
300  }
301  else{
302 
303  const Eos_poly& eos = dynamic_cast<const Eos_poly&>( eos_i ) ;
304 
305  if (eos.gam != gam) {
306  cout
307  << "The two Eos_poly have different gamma : " << gam << " <-> "
308  << eos.gam << endl ;
309  resu = false ;
310  }
311 
312  if (eos.kap != kap) {
313  cout
314  << "The two Eos_poly have different kappa : " << kap << " <-> "
315  << eos.kap << endl ;
316  resu = false ;
317  }
318 
319  if (eos.m_0 != m_0) {
320  cout
321  << "The two Eos_poly have different m_0 : " << m_0 << " <-> "
322  << eos.m_0 << endl ;
323  resu = false ;
324  }
325 
326  if (eos.mu_0 != mu_0) {
327  cout
328  << "The two Eos_poly have different mu_0 : " << mu_0 << " <-> "
329  << eos.mu_0 << endl ;
330  resu = false ;
331  }
332 
333  }
334 
335  return resu ;
336 
337 }
338 
339 bool Eos_poly::operator!=(const Eos& eos_i) const {
340 
341  return !(operator==(eos_i)) ;
342 
343 }
344 
345 
346  //------------//
347  // Outputs //
348  //------------//
349 
350 void Eos_poly::sauve(FILE* fich) const {
351 
352  Eos::sauve(fich) ;
353 
354  fwrite_be(&gam, sizeof(double), 1, fich) ;
355  fwrite_be(&kap, sizeof(double), 1, fich) ;
356  double tempo = - m_0 ;
357  fwrite_be(&tempo, sizeof(double), 1, fich) ;
358  fwrite_be(&mu_0, sizeof(double), 1, fich) ;
359 
360 }
361 
362 ostream& Eos_poly::operator>>(ostream & ost) const {
363 
364  ost << "EOS of class Eos_poly (relativistic polytrope) : " << endl ;
365  ost << " Adiabatic index gamma : " << gam << endl ;
366  ost << " Pressure coefficient kappa : " << kap <<
367  " rho_nuc c^2 / n_nuc^gamma" << endl ;
368  ost << " Mean particle mass : " << m_0 << " m_B" << endl ;
369  ost << " Relativistic chemical potential at zero pressure : " << mu_0 << " m_B c^2" << endl ;
370 
371  return ost ;
372 
373 }
374 
375 
376  //------------------------------//
377  // Computational routines //
378  //------------------------------//
379 
380 // Baryon density from enthalpy
381 //------------------------------
382 
383 double Eos_poly::nbar_ent_p(double ent, const Param* ) const {
384 
385  if ( ent > ent_0 ) {
386 
387  double xx = ent - ent_0 ;
388  return pow( gam1sgamkap * rel_mu_0 * expm1(xx), unsgam1 ) ;
389  }
390  else{
391  return 0 ;
392  }
393 }
394 
395 // Energy density from enthalpy
396 //------------------------------
397 
398 double Eos_poly::ener_ent_p(double ent, const Param* ) const {
399 
400  double nn = nbar_ent_p(ent) ;
401  double pp = kap * pow( nn, gam ) ;
402 
403  return unsgam1 * pp + mu_0 * nn ;
404 
405 }
406 
407 // Pressure from enthalpy
408 //------------------------
409 
410 double Eos_poly::press_ent_p(double ent, const Param* ) const {
411 
412  double nn = nbar_ent_p(ent) ;
413 
414  return kap * pow( nn, gam ) ;
415 
416 }
417 
418 // dln(n)/ln(H) from enthalpy
419 //---------------------------
420 
421 double Eos_poly::der_nbar_ent_p(double ent, const Param* ) const {
422 
423  if ( ent > ent_0 )
424  {
425  double xx = ent - ent_0 ;
426  return ent / (- expm1(-xx)) / gam1 ;
427  }
428  else{
429  return double(1) / gam1 ; // to ensure continuity at ent=0
430  }
431 }
432 
433 // dln(e)/ln(H) from enthalpy
434 //---------------------------
435 
436 double Eos_poly::der_ener_ent_p(double ent, const Param* ) const {
437 
438  if ( ent > ent_0 ) {
439 
440  double nn = nbar_ent_p(ent) ;
441  double pp = kap * pow( nn, gam ) ;
442  double ee = unsgam1 * pp + mu_0 * nn ;
443 
444  double xx = ent - ent_0 ;
445  return ent / (- expm1(-xx)) / gam1 * ( double(1) + pp / ee) ;
446  }
447  else{
448  return double(1) / gam1 ; // to ensure continuity at ent=0
449  }
450 }
451 
452 // dln(p)/ln(H) from enthalpy
453 //---------------------------
454 
455 double Eos_poly::der_press_ent_p(double ent, const Param* ) const {
456 
457  if ( ent > ent_0 )
458  {
459  double xx = ent - ent_0 ;
460  return gam * ent / (- expm1(-xx)) / gam1 ;
461  }
462  else{
463  return gam / gam1 ; // to ensure continuity at ent=0
464  }
465 }
466 
467 // Sound speed from enthalpy
468 //---------------------------------
469 double Eos_poly::csound_square_ent_p(double ent, const Param* ) const {
470 
471  if ( ent > ent_0 ) {
472  double xx = ent - ent_0 ;
473  return gam1*( - expm1(-xx) ) ;
474  }
475  else
476  return 0. ;
477 }
478 
479 }
virtual double ener_ent_p(double ent, const Param *par=0x0) const
Computes the total energy density from the log-enthalpy.
Definition: eos_poly.C:398
virtual double der_press_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
Definition: eos_poly.C:455
double unsgam1
Definition: eos.h:842
Cmp log(const Cmp &)
Neperian logarithm.
Definition: cmp_math.C:299
Lorene prototypes.
Definition: app_hor.h:67
Equation of state base class.
Definition: eos.h:209
virtual int identify() const =0
Returns a number to identify the sub-classe of Eos the object belongs to.
virtual ~Eos_poly()
Destructor.
Definition: eos_poly.C:230
double rel_mu_0
Definition: eos.h:844
virtual double der_nbar_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
Definition: eos_poly.C:421
double ent_0
Enthalpy at zero pressure ( )
Definition: eos.h:845
double get_gam() const
Returns the adiabatic index (cf. Eq. (3))
Definition: eos_poly.C:271
double kap
Pressure coefficient (cf.
Definition: eos.h:826
virtual void sauve(FILE *) const
Save in a file.
Definition: eos.C:189
double gam1sgamkap
Definition: eos.h:843
double gam1
Definition: eos.h:841
virtual ostream & operator>>(ostream &) const
Operator >>
Definition: eos_poly.C:362
double m_0
Individual particule mass (cf.
Definition: eos.h:831
void set_auxiliary()
Computes the auxiliary quantities gam1 , unsgam1 , gam1sgamkap from the values of gam and kap...
Definition: eos_poly.C:257
double get_kap() const
Returns the pressure coefficient (cf.
Definition: eos_poly.C:275
double get_m_0() const
Return the individual particule mass (cf.
Definition: eos_poly.C:279
void operator=(const Eos_poly &)
Assignment to another Eos_poly.
Definition: eos_poly.C:239
virtual double nbar_ent_p(double ent, const Param *par=0x0) const
Computes the baryon density from the log-enthalpy.
Definition: eos_poly.C:383
Parameter storage.
Definition: param.h:125
Polytropic equation of state (relativistic case).
Definition: eos.h:812
virtual void sauve(FILE *) const
Save in a file.
Definition: eos_poly.C:350
virtual double csound_square_ent_p(double ent, const Param *par=0x0) const
Computes the sound speed squared from the enthapy with extra parameters (virtual function implemente...
Definition: eos_poly.C:469
Eos_poly(double gamma, double kappa)
Standard constructor (sets both m_0 and mu_0 to 1).
Definition: eos_poly.C:141
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
virtual double press_ent_p(double ent, const Param *par=0x0) const
Computes the pressure from the log-enthalpy.
Definition: eos_poly.C:410
double get_mu_0() const
Return the relativistic chemical potential at zero pressure [unit: , with ].
Definition: eos_poly.C:283
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
virtual int identify() const
Returns a number to identify the sub-classe of Eos the object belongs to.
double gam
Adiabatic index (cf. Eq. (3))
Definition: eos.h:819
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
virtual bool operator==(const Eos &) const
Comparison operator (egality)
Definition: eos_poly.C:293
void set_name(const char *name_i)
Sets the EOS name.
Definition: eos.C:173
virtual bool operator!=(const Eos &) const
Comparison operator (difference)
Definition: eos_poly.C:339
char name[100]
EOS name.
Definition: eos.h:215
double mu_0
Relativistic chemical potential at zero pressure [unit: , with ].
Definition: eos.h:837
virtual double der_ener_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
Definition: eos_poly.C:436