LORENE
dyneos_poly.C
1 /*
2  * Methods of the class Dyn_eos_poly.
3  *
4  * (see file dyneos.h for documentation).
5  */
6 
7 /*
8  * Copyright (c) 2019 Jerome Novak
9  * (c) 2000 Eric Gourgoulhon for Eos classes
10  *
11  * This file is part of LORENE.
12  *
13  * LORENE is free software; you can redistribute it and/or modify
14  * it under the terms of the GNU General Public License as published by
15  * the Free Software Foundation; either version 2 of the License, or
16  * (at your option) any later version.
17  *
18  * LORENE is distributed in the hope that it will be useful,
19  * but WITHOUT ANY WARRANTY; without even the implied warranty of
20  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21  * GNU General Public License for more details.
22  *
23  * You should have received a copy of the GNU General Public License
24  * along with LORENE; if not, write to the Free Software
25  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26  *
27  */
28 
29 
30 
31 
32 /*
33  * $Id: dyneos_poly.C,v 1.3 2020/12/17 17:00:27 j_novak Exp $
34  * $Log: dyneos_poly.C,v $
35  * Revision 1.3 2020/12/17 17:00:27 j_novak
36  * Output of sound speed squared, instead of sound speed.
37  *
38  * Revision 1.2 2019/12/19 13:38:37 e_declerck
39  * Correction pour la formaule de la vitesse du son dans dyneos_poly.C
40  *
41  * Revision 1.1 2019/12/06 14:30:50 j_novak
42  * New classes Dyn_eos... for cold Eos's with baryon density as input.
43  *
44  *
45  *
46  * $Header: /cvsroot/Lorene/C++/Source/Eos/dyneos_poly.C,v 1.3 2020/12/17 17:00:27 j_novak Exp $
47  *
48  */
49 
50 // Headers C
51 #include <cstdlib>
52 #include <cstring>
53 #include <cmath>
54 
55 // Headers Lorene
56 #include "dyneos.h"
57 #include "utilitaires.h"
58 
59 namespace Lorene {
60 
61  //--------------//
62  // Constructors //
63  //--------------//
64 
65  // Standard constructor with m_0 = 1 and mu_0 = 1
66  // -----------------------------------------------
67  Dyn_eos_poly::Dyn_eos_poly(double gam0, double kappa) :
68  Dyn_eos("Relativistic polytropic EOS"),
69  gam(gam0), kap(kappa), m_0(double(1)), mu_0(double(1))
70  {
71  set_auxiliary() ;
72  }
73 
74  // Standard constructor with mu_0 = 1
75  // ----------------------------------
76  Dyn_eos_poly::Dyn_eos_poly(double gam0, double kappa, double mass) :
77  Dyn_eos("Relativistic polytropic EOS"),
78  gam(gam0), kap(kappa), m_0(mass), mu_0(double(1))
79  {
80  set_auxiliary() ;
81  }
82 
83  // Standard constructor with mu_0 = 1
84  // ----------------------------------
85  Dyn_eos_poly::Dyn_eos_poly(double gam0, double kappa, double mass, double mu_zero) :
86  Dyn_eos("Relativistic polytropic EOS"),
87  gam(gam0), kap(kappa), m_0(mass), mu_0(mu_zero)
88  {
89  set_auxiliary() ;
90  }
91 
92  // Copy constructor
93  // ----------------
95  Dyn_eos(eosi),
96  gam(eosi.gam), kap(eosi.kap), m_0(eosi.m_0), mu_0(eosi.mu_0)
97  {
98  set_auxiliary() ;
99  }
100 
101 
102  // Constructor from binary file
103  // ----------------------------
105  {
106  fread_be(&gam, sizeof(double), 1, fich) ;
107  fread_be(&kap, sizeof(double), 1, fich) ;
108  fread_be(&m_0, sizeof(double), 1, fich) ;
109 
110  if (m_0 < 0) // to ensure compatibility with previous version (revision <= 1.2)
111  { // of Eos_poly
112  m_0 = fabs( m_0 ) ;
113  fread_be(&mu_0, sizeof(double), 1, fich) ;
114  }
115  else
116  {
117  mu_0 = double(1) ;
118  }
119  set_auxiliary() ;
120  }
121 
122 
123  // Constructor from a formatted file
124  // ---------------------------------
125  Dyn_eos_poly::Dyn_eos_poly(ifstream& fich) : Dyn_eos(fich) {
126 
127  fich >> gam ; fich.ignore(1000, '\n') ;
128  fich >> kap ; fich.ignore(1000, '\n') ;
129  fich >> m_0 ; fich.ignore(1000, '\n') ;
130 
131  if (m_0 < 0) // to ensure compatibility with previous version (revision <= 1.2)
132  { // of Eos_poly
133  m_0 = fabs( m_0 ) ;
134  fich >> mu_0 ; fich.ignore(1000, '\n') ;
135  }
136  else
137  {
138  mu_0 = double(1) ;
139  }
140  set_auxiliary() ;
141  }
142  //--------------//
143  // Destructor //
144  //--------------//
145 
146  Dyn_eos_poly::~Dyn_eos_poly(){} // does nothing
147 
148  //--------------//
149  // Assignment //
150  //--------------//
151 
153 
154  set_name(eosi.name) ;
155 
156  gam = eosi.gam ;
157  kap = eosi.kap ;
158  m_0 = eosi.m_0 ;
159  mu_0 = eosi.mu_0 ;
160 
161  set_auxiliary() ;
162 
163 }
164 
165 
166  //-----------------------//
167  // Miscellaneous //
168  //-----------------------//
169 
171  {
172  gam1 = gam - double(1) ;
173  kapsgam1 = kap / gam1 ;
174  gamkapsgam1 = gam * kap / (gam1*m_0) ;
175  rel_mu_0 = mu_0 / m_0 ;
176  }
177 
178  double Dyn_eos_poly::get_gam() const
179  {
180  return gam ;
181  }
182 
183  double Dyn_eos_poly::get_kap() const
184  {
185  return kap ;
186  }
187 
188  double Dyn_eos_poly::get_m_0() const
189  {
190  return m_0 ;
191  }
192 
193  double Dyn_eos_poly::get_mu_0() const
194  {
195  return mu_0 ;
196  }
197 
198 
199  //------------------------//
200  // Comparison operators //
201  //------------------------//
202 
203 
204  bool Dyn_eos_poly::operator==(const Dyn_eos& eos_i) const
205  {
206  bool resu = true ;
207 
208  if ( eos_i.identify() != identify() )
209  {
210  cout << "The second EOS is not of type Dyn_eos_poly !" << endl ;
211  resu = false ;
212  }
213  else
214  {
215  const Dyn_eos_poly& eos = dynamic_cast<const Dyn_eos_poly&>( eos_i ) ;
216  if (eos.gam != gam)
217  {
218  cout
219  << "The two Dyn_eos_poly have different gamma : " << gam << " <-> "
220  << eos.gam << endl ;
221  resu = false ;
222  }
223  if (eos.kap != kap)
224  {
225  cout
226  << "The two Dyn_eos_poly have different kappa : " << kap << " <-> "
227  << eos.kap << endl ;
228  resu = false ;
229  }
230  if (eos.m_0 != m_0)
231  {
232  cout
233  << "The two Dyn_eos_poly have different m_0 : " << m_0 << " <-> "
234  << eos.m_0 << endl ;
235  resu = false ;
236  }
237  if (eos.mu_0 != mu_0)
238  {
239  cout
240  << "The two Dyn_eos_poly have different mu_0 : " << mu_0 << " <-> "
241  << eos.mu_0 << endl ;
242  resu = false ;
243  }
244  }
245  return resu ;
246  }
247 
248  bool Dyn_eos_poly::operator!=(const Dyn_eos& eos_i) const
249  {
250  return !(operator==(eos_i)) ;
251  }
252 
253 
254  //------------//
255  // Outputs //
256  //------------//
257 
258  void Dyn_eos_poly::sauve(FILE* fich) const
259  {
260  Dyn_eos::sauve(fich) ;
261 
262  fwrite_be(&gam, sizeof(double), 1, fich) ;
263  fwrite_be(&kap, sizeof(double), 1, fich) ;
264  double tempo = - m_0 ;
265  fwrite_be(&tempo, sizeof(double), 1, fich) ;
266  fwrite_be(&mu_0, sizeof(double), 1, fich) ;
267  }
268 
269  ostream& Dyn_eos_poly::operator>>(ostream & ost) const
270  {
271  ost << "EOS of class Dyn_eos_poly (relativistic polytrope) : " << endl ;
272  ost << " Adiabatic index gamma : " << gam << endl ;
273  ost << " Pressure coefficient kappa : " << kap <<
274  " rho_nuc c^2 / n_nuc^gamma" << endl ;
275  ost << " Mean particle mass : " << m_0 << " m_B" << endl ;
276  ost << " Relativistic chemical potential at zero pressure : "
277  << mu_0 << " m_B c^2" << endl ;
278  return ost ;
279  }
280 
281 
282  //------------------------------//
283  // Computational routines //
284  //------------------------------//
285 
286  // Baryon log-enthalpy from baryon density
287  //-----------------------------------------
288  double Dyn_eos_poly::ent_nbar_p(double nbar, const Param* ) const
289  {
290  if ( nbar > 0. )
291  return log( gamkapsgam1*pow(nbar, gam1) + rel_mu_0 ) ;
292  else
293  return 0. ;
294  }
295 
296  // Energy density from baryon density
297  //------------------------------------
298  double Dyn_eos_poly::ener_nbar_p(double nbar, const Param* ) const
299  {
300  if ( nbar > 0. )
301  return kapsgam1*pow(nbar, gam) + mu_0*nbar ;
302  else
303  return 0. ;
304  }
305 
306  // Pressure from baryon density
307  //------------------------------
308  double Dyn_eos_poly::press_nbar_p(double nbar, const Param* ) const
309  {
310  if ( nbar > 0. )
311  return kap * pow( nbar, gam ) ;
312  else
313  return 0 ;
314  }
315 
316  // Sound speed from baryon density
317  //---------------------------------
318  double Dyn_eos_poly::csound_square_nbar_p(double nbar, const Param* ) const
319  {
320  if ( nbar > 0. )
321  {
322  double ngam = pow(nbar, gam1) ;
323  return kap*gam*ngam / ( gam*kapsgam1*ngam + mu_0 ) ;
324  }
325  else
326  return 0. ;
327  }
328 
329 
330 }
virtual double ent_nbar_p(double nbar, const Param *par=0x0) const
Computes the log-enthalpy from the baryon density and extra parameters (virtual function implemented ...
Definition: dyneos_poly.C:288
void set_name(const string &)
Sets the EOS name.
Definition: dyneos.C:105
virtual ~Dyn_eos_poly()
Destructor.
Definition: dyneos_poly.C:146
Cmp log(const Cmp &)
Neperian logarithm.
Definition: cmp_math.C:299
double mu_0
Relativistic chemical potential at zero pressure [unit: , with ].
Definition: dyneos.h:403
virtual int identify() const
Returns a number to identify the sub-classe of Dyn_eos the object belongs to.
Definition: dyneos.C:313
virtual double csound_square_nbar_p(double nbar, const Param *par=0x0) const
Computes the sound speed squared from the baryon density with extra parameters.
Definition: dyneos_poly.C:318
virtual void sauve(FILE *) const
Save in a file.
Definition: dyneos.C:178
virtual bool operator==(const Dyn_eos &) const
Comparison operator (egality)
Definition: dyneos_poly.C:204
virtual int identify() const =0
Returns a number to identify the sub-classe of Dyn_eos the object belongs to.
Lorene prototypes.
Definition: app_hor.h:67
Dyn_eos_poly(double gamma, double kappa)
Standard constructor (sets both m_0 and mu_0 to 1).
Definition: dyneos_poly.C:67
virtual double press_nbar_p(double nbar, const Param *par=0x0) const
Computes the pressure from the baryon density and extra parameters (virtual function implemented in t...
Definition: dyneos_poly.C:308
Polytropic equation of state (relativistic case) for use in dynamical code.
Definition: dyneos.h:378
void set_auxiliary()
Computes the auxiliary quantities gam1 , unsgam1 , gam1sgamkap from the values of gam and kap...
Definition: dyneos_poly.C:170
double get_gam() const
Returns the adiabatic index (cf. Eq. (3))
Definition: dyneos_poly.C:178
double kapsgam1
Definition: dyneos.h:408
virtual double ener_nbar_p(double nbar, const Param *par=0x0) const
Computes the total energy density from the baryon density and extra parameters (virtual function impl...
Definition: dyneos_poly.C:298
void operator=(const Dyn_eos_poly &)
Assignment to another Dyn_eos_poly.
Definition: dyneos_poly.C:152
Equation of state for use in dynamical code base class.
Definition: dyneos.h:75
virtual ostream & operator>>(ostream &) const
Operator >>
Definition: dyneos_poly.C:269
double get_mu_0() const
Return the relativistic chemical potential at zero pressure [unit: , with ].
Definition: dyneos_poly.C:193
double gam
Adiabatic index (cf. Eq. (3))
Definition: dyneos.h:385
double gamkapsgam1
Definition: dyneos.h:409
double get_kap() const
Returns the pressure coefficient (cf.
Definition: dyneos_poly.C:183
virtual bool operator!=(const Dyn_eos &) const
Comparison operator (difference)
Definition: dyneos_poly.C:248
Parameter storage.
Definition: param.h:125
double get_m_0() const
Return the individual particule mass (cf.
Definition: dyneos_poly.C:188
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
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
string name
EOS name.
Definition: dyneos.h:81
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
double kap
Pressure coefficient (cf.
Definition: dyneos.h:392
double rel_mu_0
Definition: dyneos.h:410
virtual void sauve(FILE *) const
Save in a file.
Definition: dyneos_poly.C:258
double m_0
Individual particule mass (cf.
Definition: dyneos.h:397