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
29char 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// -----------------------------------------------
125namespace Lorene {
126Eos_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// ----------------------------------
136Eos_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// ----------------------------------
146Eos_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// ----------------------------
167Eos_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// ---------------------------------
190Eos_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
224void 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
256double Eos_poly::get_gam() const {
257 return gam ;
258}
259
260double Eos_poly::get_kap() const {
261 return kap ;
262}
263
264double Eos_poly::get_m_0() const {
265 return m_0 ;
266}
267
268double Eos_poly::get_mu_0() const {
269 return mu_0 ;
270}
271
272
273 //------------------------//
274 // Comparison operators //
275 //------------------------//
276
277
278bool 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
324bool Eos_poly::operator!=(const Eos& eos_i) const {
325
326 return !(operator==(eos_i)) ;
327
328}
329
330
331 //------------//
332 // Outputs //
333 //------------//
334
335void 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
347ostream& 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
368double 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
382double 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
400double 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
418double 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
438double 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
469double 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}
Polytropic equation of state (relativistic case).
Definition eos.h:757
virtual int identify() const
Returns a number to identify the sub-classe of Eos the object belongs to.
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
double ent_0
Enthalpy at zero pressure ( )
Definition eos.h:790
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
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
double kap
Pressure coefficient (cf.
Definition eos.h:771
virtual void sauve(FILE *) const
Save in a file.
Definition eos_poly.C:335
Eos_poly(double gamma, double kappa)
Standard constructor (sets both m_0 and mu_0 to 1).
Definition eos_poly.C:126
double get_mu_0() const
Return the relativistic chemical potential at zero pressure [unit: , with ].
Definition eos_poly.C:268
double get_gam() const
Returns the adiabatic index (cf. Eq. (3))
Definition eos_poly.C:256
double gam1sgamkap
Definition eos.h:788
double unsgam1
Definition eos.h:787
double gam
Adiabatic index (cf. Eq. (3))
Definition eos.h:764
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
void operator=(const Eos_poly &)
Assignment to another Eos_poly.
Definition eos_poly.C:224
virtual double press_ent_p(double ent, const Param *par=0x0) const
Computes the pressure from the log-enthalpy.
Definition eos_poly.C:400
double m_0
Individual particule mass (cf.
Definition eos.h:776
double mu_0
Relativistic chemical potential at zero pressure [unit: , with ].
Definition eos.h:782
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
virtual bool operator!=(const Eos &) const
Comparison operator (difference)
Definition eos_poly.C:324
virtual bool operator==(const Eos &) const
Comparison operator (egality)
Definition eos_poly.C:278
virtual ~Eos_poly()
Destructor.
Definition eos_poly.C:215
double get_m_0() const
Return the individual particule mass (cf.
Definition eos_poly.C:264
double get_kap() const
Returns the pressure coefficient (cf.
Definition eos_poly.C:260
double rel_mu_0
Definition eos.h:789
double gam1
Definition eos.h:786
virtual ostream & operator>>(ostream &) const
Operator >>
Definition eos_poly.C:347
void set_auxiliary()
Computes the auxiliary quantities gam1 , unsgam1 , gam1sgamkap from the values of gam and kap.
Definition eos_poly.C:242
Equation of state base class.
Definition eos.h:190
virtual int identify() const =0
Returns a number to identify the sub-classe of Eos the object belongs to.
virtual void sauve(FILE *) const
Save in a file.
Definition eos.C:179
char name[100]
EOS name.
Definition eos.h:196
void set_name(const char *name_i)
Sets the EOS name.
Definition eos.C:163
Parameter storage.
Definition param.h:125
Cmp exp(const Cmp &)
Exponential.
Definition cmp_math.C:270
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
Cmp log(const Cmp &)
Neperian logarithm.
Definition cmp_math.C:296
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
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 prototypes.
Definition app_hor.h:64