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