LORENE
piecewise_polytrope_1D.C
1 /*
2  * Methods of the class Piecewise_polytrope_1D.
3  *
4  * (see file eos.h for documentation).
5  */
6 
7 /*
8  * Copyright (c) 2023 GaĆ«l Servignat
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 
33  */
34 
35 // Headers C
36 #include <cstdlib>
37 #include <cstring>
38 #include <cmath>
39 
40 // Headers Lorene
41 #include "eos.h"
42 #include "cmp.h"
43 #include "unites.h"
44 
45  //--------------//
46  // Constructors //
47  //--------------//
48 
49 // Standard constructor
50 // --------------------
51 namespace Lorene {
52 void huntm(const Tbl& xx, double& x, int& i_low) ;
53 
54 Piecewise_polytrope_1D::Piecewise_polytrope_1D(const Tbl& Gamma, const Tbl& Kappa, const Tbl& Lambda,
55  const Tbl& AA, const Tbl& N_lim, const Tbl& Ent_lim, double gamma, double kappa, double n_lim0):
56  gamma_high(0x0), kappa_high(0x0), Lambda_high(0x0), a_high(0x0), n_lim_high(0x0), ent_lim_high(0x0), gamma_low(gamma), kappa_low(kappa), n_lim(n_lim0) {
57 
58  set_name("Pseudo-polytropic fit of cold, beta-equilibrated EoS") ;
59 
60  assert(Gamma.get_ndim() == 1) ;
61  assert(Kappa.get_ndim() == 1) ;
62  assert(Lambda.get_ndim() == 1) ;
63  assert(AA.get_ndim() == 1) ;
64  assert(N_lim.get_ndim() == 1) ;
65  assert(Ent_lim.get_ndim() == 1) ;
66  gamma_high = new Tbl(Gamma) ;
67  n_param_high = gamma_high->get_taille() ;
68  assert(Kappa.get_taille() == n_param_high) ;
69  kappa_high = new Tbl(Kappa) ;
70  assert(Lambda.get_taille() == n_param_high) ;
71  Lambda_high = new Tbl(Lambda) ;
72  assert(AA.get_taille() == n_param_high) ;
73  a_high = new Tbl(AA) ;
74  assert(N_lim.get_taille() == n_param_high) ;
75  n_lim_high = new Tbl(N_lim) ;
76  assert(Ent_lim.get_taille() == n_param_high) ;
77  ent_lim_high = new Tbl(Ent_lim) ; // maybe ent_lim should be computed here ? Anyway this constructor is very likely not to be used a lot so maybe it is not a problem.
78 
79  eos_low = new Eos_poly(gamma_low, kappa_low) ;
80 }
81 
82 
83 // Copy constructor
84 // ----------------
85 Piecewise_polytrope_1D::Piecewise_polytrope_1D(const Piecewise_polytrope_1D& eosi): Eos(eosi), gamma_high(eosi.gamma_high), kappa_high(eosi.kappa_high),
86  Lambda_high(eosi.Lambda_high), a_high(eosi.a_high), n_lim_high(eosi.n_lim_high), ent_lim_high(eosi.ent_lim_high), gamma_low(eosi.gamma_low),
87  kappa_low(eosi.kappa_low), n_lim(eosi.n_lim), ent_lim(eosi.ent_lim), n_param_high(eosi.n_param_high) {}
88 
89 
90 // Constructor from a binary file
91 // ------------------------------
93 
94  fread_be(&n_param_high, sizeof(int), 1, fich) ;
95  gamma_high = new Tbl(fich) ;
96  kappa_high = new Tbl(fich) ;
97  Lambda_high = new Tbl(fich) ;
98  a_high = new Tbl(fich) ;
99  n_lim_high = new Tbl(fich) ;
100  ent_lim_high = new Tbl(fich) ;
101  fread_be(&kappa_low, sizeof(double), 1, fich) ;
102  fread_be(&gamma_low, sizeof(double), 1, fich) ;
103  fread_be(&n_lim, sizeof(double), 1, fich) ;
104  fread_be(&ent_lim, sizeof(double), 1, fich) ;
105 
106  eos_low = new Eos_poly(gamma_low, kappa_low) ;
107 }
108 
109 // Constructor from a formatted file
110 // ---------------------------------
112 
113  fich >> n_param_high ; fich.ignore(1000, '\n') ;
114 
115  gamma_high = new Tbl(n_param_high) ; gamma_high->set_etat_qcq() ;
116  for (int i=0; i < n_param_high; i++)
117  fich >> gamma_high->set(i) ;
118  fich.ignore(1000, '\n') ;
119 
120  kappa_high = new Tbl(n_param_high) ; kappa_high->set_etat_qcq() ;
121  for (int i=0; i < n_param_high; i++)
122  fich >> kappa_high->set(i) ;
123  fich.ignore(1000, '\n') ;
124 
125  Lambda_high = new Tbl(n_param_high) ; Lambda_high->set_etat_qcq() ;
126  for (int i=0; i < n_param_high; i++)
127  fich >> Lambda_high->set(i) ;
128  fich.ignore(1000, '\n') ;
129 
130  a_high = new Tbl(n_param_high) ; a_high->set_etat_qcq() ;
131  for (int i=0; i < n_param_high; i++)
132  fich >> a_high->set(i) ;
133  fich.ignore(1000, '\n') ;
134 
135  n_lim_high = new Tbl(n_param_high) ; n_lim_high->set_etat_qcq() ;
136  for (int i=0; i < n_param_high; i++)
137  fich >> n_lim_high->set(i) ;
138  fich.ignore(1000, '\n') ;
139 
140  ent_lim_high = new Tbl(n_param_high) ; ent_lim_high->set_etat_qcq() ;
141  for (int i=0; i < n_param_high; i++)
142  fich >> ent_lim_high->set(i) ;
143  fich.ignore(1000, '\n') ;
144 
145  n_lim = (*n_lim_high)(0) ;
146  ent_lim = (*ent_lim_high)(0) ;
147 
148  fich >> kappa_low ; fich.ignore(1000, '\n') ;
149 
150  fich >> gamma_low ; fich.ignore(1000, '\n') ;
151 
152  // double Gamma1 = (*gamma_high)(0) ;
153  // double Kappa1 = (*kappa_high)(0) ;
154  // double a1 = (*a_high)(0) ;
155  // ent_lim = log(1. + a1 + Kappa1*Gamma1/(Gamma1-1.) * pow(0.1*n_lim,Gamma1-1.)) ;
156 
157  eos_low = new Eos_poly(gamma_low, kappa_low) ;
158 
159 }
160 
161 
162  //--------------//
163  // Destructor //
164  //--------------//
165 
167 
168  if (gamma_high != 0x0) delete gamma_high ;
169  if (kappa_high != 0x0) delete kappa_high ;
170  if (Lambda_high != 0x0) delete Lambda_high ;
171  if (a_high != 0x0) delete a_high ;
172  if (n_lim_high != 0x0) delete n_lim_high ;
173  if (ent_lim_high != 0x0) delete ent_lim_high ;
174  if (eos_low != 0x0) delete eos_low ;
175 
176 }
177  //--------------//
178  // Assignment //
179  //--------------//
180 
182 
183  set_name(eosi.name) ;
184 
185  gamma_high = eosi.gamma_high ;
186  kappa_high = eosi.kappa_high ;
187  Lambda_high = eosi.Lambda_high ;
188  a_high = eosi.a_high ;
189  n_lim_high = eosi.n_lim_high ;
190  ent_lim_high = eosi.ent_lim_high ;
191  n_lim = eosi.n_lim ;
192  n_param_high = eosi.n_param_high ;
193  gamma_low = eosi.gamma_low ;
194  kappa_low = eosi.kappa_low ;
195  ent_lim = eosi.ent_lim ;
196  eos_low = eosi.eos_low ;
197 
198 }
199  //------------------------//
200  // Comparison operators //
201  //------------------------//
202 
203 bool Piecewise_polytrope_1D::operator==(const Eos& eos_i) const {
204 
205  bool resu = true ;
206 
207  if ( eos_i.identify() != identify() ) {
208  cout << "The second EOS is not of type Piecewise_polytrope_1D !" << endl ;
209  resu = false ;
210  }
211  else{
212 
213  const Piecewise_polytrope_1D& eos = dynamic_cast<const Piecewise_polytrope_1D&>( eos_i ) ;
214 
215  for (int i=0; i < n_param_high; i++){
216  if ((*eos.gamma_high)(i) != (*gamma_high)(i)) {
217  cout
218  << "The two Piecewise_polytrope_1D have different coefficients: " << (*gamma_high)(i) << " <-> "
219  << (*eos.gamma_high)(i) << endl ;
220  resu = false ;
221  }
222 
223  if ((*eos.kappa_high)(i) != (*kappa_high)(i)) {
224  cout
225  << "The two Piecewise_polytrope_1D have different coefficients: " << (*kappa_high)(i) << " <-> "
226  << (*eos.kappa_high)(i) << endl ;
227  resu = false ;
228  }
229 
230  if ((*eos.Lambda_high)(i) != (*Lambda_high)(i)) {
231  cout
232  << "The two Piecewise_polytrope_1D have different coefficients: " << (*Lambda_high)(i) << " <-> "
233  << (*eos.Lambda_high)(i) << endl ;
234  resu = false ;
235  }
236 
237  if ((*eos.a_high)(i) != (*a_high)(i)) {
238  cout
239  << "The two Piecewise_polytrope_1D have different coefficients: " << (*a_high)(i) << " <-> "
240  << (*eos.a_high)(i) << endl ;
241  resu = false ;
242  }
243 
244  if ((*eos.n_lim_high)(i) != (*n_lim_high)(i)) {
245  cout
246  << "The two Piecewise_polytrope_1D have different coefficients: " << (*n_lim_high)(i) << " <-> "
247  << (*eos.n_lim_high)(i) << endl ;
248  resu = false ;
249  }
250 
251  if ((*eos.ent_lim_high)(i) != (*ent_lim_high)(i)) {
252  cout
253  << "The two Piecewise_polytrope_1D have different coefficients: " << (*ent_lim_high)(i) << " <-> "
254  << (*eos.ent_lim_high)(i) << endl ;
255  resu = false ;
256  }
257  }
258 
259  if (eos.n_lim != n_lim) {
260  cout
261  << "The two Piecewise_polytrope_1D have different limiting densities n_lim: " << n_lim << " <-> "
262  << eos.n_lim << endl ;
263  resu = false ;
264  }
265 
266  if (eos.ent_lim != ent_lim) {
267  cout
268  << "The two Piecewise_polytrope_1D have different limiting enthalpies ent_lim: " << ent_lim << " <-> "
269  << eos.ent_lim << endl ;
270  resu = false ;
271  }
272 
273  if (eos.n_param_high != n_param_high) {
274  cout
275  << "The two Piecewise_polytrope_1D have different number of coefficients: " << n_param_high << " <-> "
276  << eos.n_param_high << endl ;
277  resu = false ;
278  }
279 
280  if (eos.kappa_low != kappa_low) {
281  cout
282  << "The two Piecewise_polytrope_1D have different kappa in low polytrope: " << kappa_low << " <-> "
283  << eos.kappa_low << endl ;
284  resu = false ;
285  }
286 
287  if (eos.gamma_low != gamma_low) {
288  cout
289  << "The two Piecewise_polytrope_1D have different kappa in low polytrope: " << gamma_low << " <-> "
290  << eos.gamma_low << endl ;
291  resu = false ;
292  }
293 
294  }
295 
296  return resu ;
297 
298 }
299 
300 bool Piecewise_polytrope_1D::operator!=(const Eos& eos_i) const {
301 
302  return !(operator==(eos_i)) ;
303 
304 }
305 
306 
307  //------------//
308  // Outputs //
309  //------------//
310 
311 void Piecewise_polytrope_1D::sauve(FILE* fich) const {
312 
313  Eos::sauve(fich) ;
314 
315  fwrite_be(&n_param_high, sizeof(int), 1, fich) ;
316  gamma_high->sauve(fich) ;
317  kappa_high->sauve(fich) ;
318  Lambda_high->sauve(fich) ;
319  a_high->sauve(fich) ;
320  n_lim_high->sauve(fich) ;
321  ent_lim_high->sauve(fich) ;
322  fwrite_be(&kappa_low, sizeof(double), 1, fich) ;
323  fwrite_be(&gamma_low, sizeof(double), 1, fich) ;
324  fwrite_be(&n_lim, sizeof(double), 1, fich) ;
325  fwrite_be(&ent_lim, sizeof(double), 1, fich) ;
326 
327 }
328 
329 ostream& Piecewise_polytrope_1D::operator>>(ostream & ost) const {
330 
331  ost << setprecision(16) << "EOS of class Piecewise_polytrope_1D (analytical fit of cold EoS): " << '\n' ;
332  ost << " Gamma_high coefficients: " << *gamma_high << '\n' ;
333  ost << " Kappa_high coefficients: " << *kappa_high << '\n' ;
334  ost << " Lambda_high coefficients: " << *Lambda_high << '\n' ;
335  ost << " a_high coefficients: " << *a_high << '\n' ;
336  ost << " n_lim_high coefficients: " << *n_lim_high << '\n' ;
337  ost << " ent_lim_high coefficients: " << *ent_lim_high << '\n' ;
338  ost << "Low densities polytrope :" << '\n' ;
339  cout << *eos_low << '\n' ;
340  ost << setprecision(16) << " Limiting density n_lim: " << 0.1*n_lim << " [fm^-3]" << '\n' ;
341  ost << setprecision(16) << " Limiting enthalpy h_lim: " << ent_lim << " [c^2]" << endl;
342 
343  return ost ;
344 
345 }
346 
347 
348  //------------------------------//
349  // Computational routines //
350  //------------------------------//
351 
352 // Baryon density from enthalpy
353 //------------------------------
354 
355 double Piecewise_polytrope_1D::nbar_ent_p(double ent, const Param* ) const {
356 
357  if (ent > ent_lim ) {
358 
359  int i_low = 0;
360  huntm(*ent_lim_high, ent, i_low) ;
361  Tbl& Gamma = *gamma_high ;
362  Tbl& Kappa = *kappa_high ;
363  Tbl& AA = *a_high ;
364  double nn = 1./(Gamma(i_low)-1.) ;
365 
366  return pow((exp(ent)-1.-AA(i_low))/Kappa(i_low)/(nn+1.), nn) ;
367 
368  }
369  else if ( (0 < ent) && (ent < ent_lim)){
370 
371  return eos_low->nbar_ent_p(ent) ;
372  }
373  else{
374  return 0 ;
375  }
376 }
377 
378 // Energy density from enthalpy
379 //------------------------------
380 
381 double Piecewise_polytrope_1D::ener_ent_p(double ent, const Param* ) const {
382  if (ent > ent_lim ) {
383 
384  int i_low = 0;
385  huntm(*ent_lim_high, ent, i_low) ;
386  Tbl& Gamma = *gamma_high ;
387  Tbl& Kappa = *kappa_high ;
388  Tbl& AA = *a_high ;
389  Tbl& Lambda = *Lambda_high ;
390  double nn = 1./(Gamma(i_low)-1.) ;
391  double nbar = pow((exp(ent)-1.-AA(i_low))/Kappa(i_low)/(nn+1.), nn) ;
392 
393 
394  return Kappa(i_low)*nn * pow(nbar, Gamma(i_low)) + (1.+AA(i_low))*nbar - Lambda(i_low) ;
395  }
396  else if ( (0 < ent) && (ent < ent_lim)){
397 
398  return eos_low->ener_ent_p(ent) ;
399  }
400  else{
401  return 0 ;
402  }
403 }
404 
405 // Pressure from enthalpy
406 //------------------------
407 
408 double Piecewise_polytrope_1D::press_ent_p(double ent, const Param* ) const {
409  if (ent > ent_lim ) {
410 
411  int i_low = 0;
412  huntm(*ent_lim_high, ent, i_low) ;
413  Tbl& Gamma = *gamma_high ;
414  Tbl& Kappa = *kappa_high ;
415  Tbl& AA = *a_high ;
416  Tbl& Lambda = *Lambda_high ;
417  double nn = 1./(Gamma(i_low)-1.) ;
418  double nbar = pow((exp(ent)-1.-AA(i_low))/Kappa(i_low)/(nn+1.), nn) ;
419 
420  return Kappa(i_low) * pow(nbar, Gamma(i_low)) + Lambda(i_low) ;
421  }
422  else if ( (0 < ent) && (ent < ent_lim)){
423 
424  return eos_low->press_ent_p(ent) ;
425  }
426  else{
427  return 0 ;
428  }
429 }
430 
431 // dln(n)/ln(h) from enthalpy
432 //---------------------------
433 
434 double Piecewise_polytrope_1D::der_nbar_ent_p(double , const Param* ) const {
435  c_est_pas_fait("Piecewise_polytrope_1D::der_nbar_ent_p") ;
436  return 0. ;
437 }
438 
439 // dln(e)/ln(h) from enthalpy
440 //---------------------------
441 
442 double Piecewise_polytrope_1D::der_ener_ent_p(double, const Param* ) const {
443  c_est_pas_fait("Piecewise_polytrope_1D::der_ener_ent_p") ;
444  return 0.;
445 }
446 
447 // dln(p)/ln(h) from enthalpy
448 //---------------------------
449 
450 double Piecewise_polytrope_1D::der_press_ent_p(double, const Param* ) const {
451  c_est_pas_fait("Piecewise_polytrope_1D::der_press_ent_p") ;
452  return 0. ;
453 }
454 
455 double Piecewise_polytrope_1D::csound_square_ent_p(double ent, const Param*) const {
456  if (ent > ent_lim ) {
457 
458  int i_low = 0;
459  huntm(*ent_lim_high, ent, i_low) ;
460  Tbl& AA = *a_high ;
461  Tbl& Gamma = *gamma_high ;
462 
463  return (Gamma(i_low)-1.)*(1.-exp(-ent)*(1.+AA(i_low))) ;
464  }
465  else if ( (0 < ent) && (ent < ent_lim)){
466 
467  return eos_low->csound_square_ent_p(ent) ;
468  }
469  else{
470  return 0 ;
471  }
472 }
473 
474 }
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 nbar_ent_p(double ent, const Param *par=0x0) const
Computes the baryon density from the specific enthalpy.
virtual double ener_ent_p(double ent, const Param *par=0x0) const
Computes the total energy density from the specific enthalpy.
virtual bool operator!=(const Eos &) const
Comparison operator (difference)
Cmp exp(const Cmp &)
Exponential.
Definition: cmp_math.C:273
Lorene prototypes.
Definition: app_hor.h:67
Equation of state base class.
Definition: eos.h:209
void operator=(const Piecewise_polytrope_1D &)
Assignment to another Piecewise_polytrope_1D.
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:301
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.
virtual double der_ener_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the specific enthalpy.
Piecewise_polytrope_1D(const Tbl &, const Tbl &, const Tbl &, const Tbl &, const Tbl &, const Tbl &, double, double, double)
Standard constructor.
virtual void sauve(FILE *) const
Save in a file.
Definition: eos.C:189
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:364
int get_ndim() const
Gives the number of dimensions (ie dim.ndim)
Definition: tbl.h:420
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 double press_ent_p(double ent, const Param *par=0x0) const
Computes the pressure from the specific enthalpy.
virtual double csound_square_ent_p(double, const Param *) const
Computes the sound speed squared from the enthapy with extra parameters (virtual function implemente...
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
virtual double der_nbar_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the specific enthalpy.
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 bool operator==(const Eos &) const
Comparison operator (egality)
virtual double press_ent_p(double ent, const Param *par=0x0) const
Computes the pressure from the log-enthalpy.
Definition: eos_poly.C:410
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
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 double der_press_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the specific enthalpy.
void c_est_pas_fait(const char *)
Helpful function to say something is not implemented yet.
void set_name(const char *name_i)
Sets the EOS name.
Definition: eos.C:173
int get_taille() const
Gives the total size (ie dim.taille)
Definition: tbl.h:417
char name[100]
EOS name.
Definition: eos.h:215
void sauve(FILE *) const
Save in a file.
Definition: tbl.C:329
Basic array class.
Definition: tbl.h:164
virtual ~Piecewise_polytrope_1D()
Destructor.
virtual ostream & operator>>(ostream &) const
Operator >>
virtual int identify() const
Returns a number to identify the sub-classe of Eos the object belongs to.