LORENE
eos_tabul.C
1 /*
2  * Methods of class Eos_tabul
3  *
4  * (see file eos.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2000-2001 Eric Gourgoulhon
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 char eos_tabul_C[] = "$Header: /cvsroot/Lorene/C++/Source/Eos/eos_tabul.C,v 1.16 2015/08/04 14:41:29 j_novak Exp $" ;
31 
32 /*
33  * $Id: eos_tabul.C,v 1.16 2015/08/04 14:41:29 j_novak Exp $
34  * $Log: eos_tabul.C,v $
35  * Revision 1.16 2015/08/04 14:41:29 j_novak
36  * Back to previous version for Eos_CompOSE. Enthalpy-consistent EoS can be accessed using Eos_consistent class (derived from Eos_CompOSE).
37  *
38  * Revision 1.15 2015/01/27 14:22:38 j_novak
39  * New methods in Eos_tabul to correct for EoS themro consistency (optional).
40  *
41  * Revision 1.14 2014/10/13 08:52:54 j_novak
42  * Lorene classes and functions now belong to the namespace Lorene.
43  *
44  * Revision 1.13 2014/06/30 16:13:18 j_novak
45  * New methods for reading directly from CompOSE files.
46  *
47  * Revision 1.12 2014/03/06 15:53:35 j_novak
48  * Eos_compstar is now Eos_compOSE. Eos_tabul uses strings and contains informations about authors.
49  *
50  * Revision 1.11 2012/11/09 10:32:44 m_bejger
51  * Implementing der_ener_ent_p
52  *
53  * Revision 1.10 2010/02/16 11:14:50 j_novak
54  * More verbose opeining of the file.
55  *
56  * Revision 1.9 2010/02/02 13:22:16 j_novak
57  * New class Eos_Compstar.
58  *
59  * Revision 1.8 2004/03/25 10:29:02 j_novak
60  * All LORENE's units are now defined in the namespace Unites (in file unites.h).
61  *
62  * Revision 1.7 2003/11/25 13:42:50 m_bejger
63  * read_table written in more ordered way
64  *
65  * Revision 1.6 2003/11/21 16:11:16 m_bejger
66  * added the computation dlnp/dlnn_b, dlnn/dlnH
67  *
68  * Revision 1.5 2003/05/30 07:50:06 e_gourgoulhon
69  *
70  * Reformulate the computation of der_nbar_ent
71  * Added computation of der_press_ent_p.
72  *
73  * Revision 1.4 2003/05/15 09:42:47 e_gourgoulhon
74  *
75  * Now computes d ln / dln H.
76  *
77  * Revision 1.3 2002/10/16 14:36:35 j_novak
78  * Reorganization of #include instructions of standard C++, in order to
79  * use experimental version 3 of gcc.
80  *
81  * Revision 1.2 2002/04/09 14:32:15 e_gourgoulhon
82  * 1/ Added extra parameters in EOS computational functions (argument par)
83  * 2/ New class MEos for multi-domain EOS
84  *
85  * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
86  * LORENE
87  *
88  * Revision 2.3 2001/09/13 13:35:49 eric
89  * Suppression des affichages dans read_table().
90  *
91  * Revision 2.2 2001/02/07 09:48:05 eric
92  * Suppression de la fonction derent_ent_p.
93  * Ajout des fonctions donnant les derivees de l'EOS:
94  * der_nbar_ent_p
95  * der_ener_ent_p
96  * der_press_ent_p
97  *
98  * Revision 2.1 2000/11/23 00:10:20 eric
99  * Enthalpie minimale fixee a 1e-14.
100  *
101  * Revision 2.0 2000/11/22 19:31:30 eric
102  * *** empty log message ***
103  *
104  *
105  * $Header: /cvsroot/Lorene/C++/Source/Eos/eos_tabul.C,v 1.16 2015/08/04 14:41:29 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 "tbl.h"
117 #include "utilitaires.h"
118 #include "unites.h"
119 
120 
121 namespace Lorene {
122 void interpol_herm(const Tbl& , const Tbl&, const Tbl&, double, int&,
123  double&, double& ) ;
124 
125 void interpol_linear(const Tbl&, const Tbl&, double, int&, double&) ;
126 
127  //----------------------------//
128  // Constructors //
129  //----------------------------//
130 
131 // Standard constructor
132 // --------------------
133 Eos_tabul::Eos_tabul(const char* name_i, const char* table,
134  const char* path) : Eos(name_i)
135 {
136  tablename = path ;
137  tablename += "/" ;
138  tablename += table ;
139 
140  read_table() ;
141 }
142 
143 // Standard constructor with full filename
144 // ---------------------------------------
145  Eos_tabul::Eos_tabul(const char* name_i, const char* file_name)
146  : Eos(name_i) {
147 
148  tablename = file_name ;
149 
150  read_table() ;
151 
152 }
153 
154 
155 // Constructor from binary file
156 // ----------------------------
157  Eos_tabul::Eos_tabul(FILE* fich) : Eos(fich) {
158 
159  char tmp_string[160] ;
160  fread(tmp_string, sizeof(char), 160, fich) ;
161  tablename = tmp_string ;
162 
163  read_table() ;
164 
165 }
166 
167 
168 
169 // Constructor from a formatted file
170 // ---------------------------------
171  Eos_tabul::Eos_tabul(ifstream& fich, const char* table) :
172  Eos(fich) {
173 
174  fich >> tablename ;
175  tablename += "/" ;
176  tablename += table ;
177 
178  read_table() ;
179 
180 }
181 
182  Eos_tabul::Eos_tabul(ifstream& fich) : Eos(fich)
183 {
184  fich >> tablename ;
185 
186  read_table() ;
187 }
188 
189 // Standard constructor with a name
190 // ---------------------------------
191  Eos_tabul::Eos_tabul(const char* name_i) : Eos(name_i),
192  logh(0x0), logp(0x0), dlpsdlh(0x0),
193  lognb(0x0), dlpsdlnb(0x0)
194 {}
195 
196 
197  //--------------//
198  // Destructor //
199  //--------------//
200 
202  if (logh != 0x0) delete logh ;
203  if (logp != 0x0) delete logp ;
204  if (dlpsdlh != 0x0) delete dlpsdlh ;
205  if (lognb != 0x0) delete lognb ;
206  if (dlpsdlnb != 0x0) delete dlpsdlnb ;
207 }
208 
209  //------------//
210  // Outputs //
211  //------------//
212 
213 void Eos_tabul::sauve(FILE* fich) const {
214 
215  Eos::sauve(fich) ;
216 
217  char tmp_string[160] ;
218  strcpy(tmp_string, tablename.c_str()) ;
219  fwrite(tmp_string, sizeof(char), 160, fich) ;
220 }
221 
222  //------------------------//
223  // Reading of the table //
224  //------------------------//
225 
227 
228  using namespace Unites ;
229 
230  ifstream fich(tablename.data()) ;
231 
232  if (!fich) {
233  cerr << "Eos_tabul::read_table(): " << endl ;
234  cerr << "Problem in opening the EOS file!" << endl ;
235  cerr << "While trying to open " << tablename << endl ;
236  cerr << "Aborting..." << endl ;
237  abort() ;
238  }
239 
240  fich.ignore(1000, '\n') ; // Jump over the first header
241  fich.ignore(1) ;
242  getline(fich, authors, '\n') ;
243  for (int i=0; i<3; i++) { // jump over the file
244  fich.ignore(1000, '\n') ; // header
245  } //
246 
247  int nbp ;
248  fich >> nbp ; fich.ignore(1000, '\n') ; // number of data
249  if (nbp<=0) {
250  cerr << "Eos_tabul::read_table(): " << endl ;
251  cerr << "Wrong value for the number of lines!" << endl ;
252  cerr << "nbp = " << nbp << endl ;
253  cerr << "Aborting..." << endl ;
254  abort() ;
255  }
256 
257  for (int i=0; i<3; i++) { // jump over the table
258  fich.ignore(1000, '\n') ; // description
259  }
260 
261  press = new double[nbp] ;
262  nb = new double[nbp] ;
263  ro = new double[nbp] ;
264 
265  logh = new Tbl(nbp) ;
266  logp = new Tbl(nbp) ;
267  dlpsdlh = new Tbl(nbp) ;
268  lognb = new Tbl(nbp) ;
269  dlpsdlnb = new Tbl(nbp) ;
270 
271  logh->set_etat_qcq() ;
272  logp->set_etat_qcq() ;
273  dlpsdlh->set_etat_qcq() ;
274  lognb->set_etat_qcq() ;
275  dlpsdlnb->set_etat_qcq() ;
276 
277  double rhonuc_cgs = rhonuc_si * 1e-3 ;
278  double c2_cgs = c_si * c_si * 1e4 ;
279 
280  int no ;
281  double nb_fm3, rho_cgs, p_cgs ;
282 
283  for (int i=0; i<nbp; i++) {
284 
285  fich >> no ;
286  fich >> nb_fm3 ;
287  fich >> rho_cgs ;
288  fich >> p_cgs ; fich.ignore(1000,'\n') ;
289  if ( (nb_fm3<0) || (rho_cgs<0) || (p_cgs < 0) ){
290  cout << "Eos_tabul::read_table(): " << endl ;
291  cout << "Negative value in table!" << endl ;
292  cout << "nb = " << nb_fm3 << ", rho = " << rho_cgs <<
293  ", p = " << p_cgs << endl ;
294  cout << "Aborting..." << endl ;
295  abort() ;
296  }
297 
298  press[i] = p_cgs / c2_cgs ;
299  nb[i] = nb_fm3 ;
300  ro[i] = rho_cgs ;
301  }
302 
303  double ww = 0. ;
304  for (int i=0; i<nbp; i++) {
305  double h = log( (ro[i] + press[i]) /
306  (10 * nb[i] * rhonuc_cgs) ) ;
307 
308  if (i==0) { ww = h ; }
309  h = h - ww + 1.e-14 ;
310 
311  logh->set(i) = log10( h ) ;
312  logp->set(i) = log10( press[i] / rhonuc_cgs ) ;
313  dlpsdlh->set(i) = h * (ro[i] + press[i]) / press[i] ;
314  lognb->set(i) = log10(nb[i]) ;
315  }
316 
317  double p0, p1, p2, n0, n1, n2, dpdnb;
318 
319  // special case: i=0
320 
321  p0 = log(press[0]);
322  p1 = log(press[1]);
323  p2 = log(press[2]);
324 
325  n0 = log(nb[0]);
326  n1 = log(nb[1]);
327  n2 = log(nb[2]);
328 
329  dpdnb = p0*(2*n0-n1-n2)/(n0-n1)/(n0-n2) +
330  p1*(n0-n2)/(n1-n0)/(n1-n2) +
331  p2*(n0-n1)/(n2-n0)/(n2-n1) ;
332 
333  dlpsdlnb->set(0) = dpdnb ;
334 
335  for(int i=1;i<nbp-1;i++) {
336 
337  p0 = log(press[i-1]);
338  p1 = log(press[i]);
339  p2 = log(press[i+1]);
340 
341  n0 = log(nb[i-1]);
342  n1 = log(nb[i]);
343  n2 = log(nb[i+1]);
344 
345  dpdnb = p0*(n1-n2)/(n0-n1)/(n0-n2) +
346  p1*(2*n1-n0-n2)/(n1-n0)/(n1-n2) +
347  p2*(n1-n0)/(n2-n0)/(n2-n1) ;
348 
349  dlpsdlnb->set(i) = dpdnb ;
350 
351  }
352 
353  // special case: i=nbp-1
354 
355  p0 = log(press[nbp-3]);
356  p1 = log(press[nbp-2]);
357  p2 = log(press[nbp-1]);
358 
359  n0 = log(nb[nbp-3]);
360  n1 = log(nb[nbp-2]);
361  n2 = log(nb[nbp-1]);
362 
363  dpdnb = p0*(n2-n1)/(n0-n1)/(n0-n2) +
364  p1*(n2-n0)/(n1-n0)/(n1-n2) +
365  p2*(2*n2-n0-n1)/(n2-n0)/(n2-n1) ;
366 
367  dlpsdlnb->set(nbp-1) = dpdnb ;
368 
369  hmin = pow( double(10), (*logh)(0) ) ;
370  hmax = pow( double(10), (*logh)(nbp-1) ) ;
371 
372 //##
373 // cout << "logh : " << endl ;
374 // cout << *logh << endl ;
375 //
376 //
377 // cout << "logp : " << endl ;
378 // cout << *logp << endl ;
379 //
380 // cout << "dlpsdlh : " << endl ;
381 // cout << *dlpsdlh << endl ;
382 //
383 // cout << "dlpsdlnb : " << endl ;
384 // cout << *dlpsdlnb << endl ;
385 //
386 //
387  cout << "hmin, hmax : " << hmin << " " << hmax << endl ;
388 //##
389 
390  fich.close();
391 
392  delete [] press ;
393  delete [] nb ;
394  delete [] ro ;
395 
396 }
397 
398 
399  //------------------------------//
400  // Computational routines //
401  //------------------------------//
402 
403 // Baryon density from enthalpy
404 //------------------------------
405 
406 double Eos_tabul::nbar_ent_p(double ent, const Param* ) const {
407 
408  static int i_near = logh->get_taille() / 2 ;
409 
410  if ( ent > hmin ) {
411  if (ent > hmax) {
412  cout << "Eos_tabul::nbar_ent_p : ent > hmax !" << endl ;
413  abort() ;
414  }
415  double logh0 = log10( ent ) ;
416  double logp0 ;
417  double dlpsdlh0 ;
418  interpol_herm(*logh, *logp, *dlpsdlh, logh0, i_near, logp0,
419  dlpsdlh0) ;
420 
421  double pp = pow(double(10), logp0) ;
422 
423  double resu = pp / ent * dlpsdlh0 * exp(-ent) ;
424  return resu ;
425  }
426  else{
427  return 0 ;
428  }
429 }
430 
431 // Energy density from enthalpy
432 //------------------------------
433 
434 double Eos_tabul::ener_ent_p(double ent, const Param* ) const {
435 
436  static int i_near = logh->get_taille() / 2 ;
437 
438  if ( ent > hmin ) {
439  if (ent > hmax) {
440  cout << "Eos_tabul::ener_ent_p : ent > hmax !" << endl ;
441  abort() ;
442  }
443  double logh0 = log10( ent ) ;
444  double logp0 ;
445  double dlpsdlh0 ;
446  interpol_herm(*logh, *logp, *dlpsdlh, logh0, i_near, logp0,
447  dlpsdlh0) ;
448 
449  double pp = pow(double(10), logp0) ;
450 
451  double resu = pp / ent * dlpsdlh0 - pp ;
452  return resu ;
453  }
454  else{
455  return 0 ;
456  }
457 }
458 
459 // Pressure from enthalpy
460 //------------------------
461 
462 double Eos_tabul::press_ent_p(double ent, const Param* ) const {
463 
464  static int i_near = logh->get_taille() / 2 ;
465 
466  if ( ent > hmin ) {
467  if (ent > hmax) {
468  cout << "Eos_tabul::press_ent_p : ent > hmax !" << endl ;
469  abort() ;
470  }
471 
472  double logh0 = log10( ent ) ;
473  double logp0 ;
474  double dlpsdlh0 ;
475  interpol_herm(*logh, *logp, *dlpsdlh, logh0, i_near, logp0,
476  dlpsdlh0) ;
477  return pow(double(10), logp0) ;
478  }
479  else{
480  return 0 ;
481  }
482 }
483 
484 // dln(n)/ln(H) from enthalpy
485 //---------------------------
486 
487 double Eos_tabul::der_nbar_ent_p(double ent, const Param* ) const {
488 
489  if ( ent > hmin ) {
490  if (ent > hmax) {
491  cout << "Eos_tabul::der_nbar_ent_p : ent > hmax !" << endl ;
492  abort() ;
493  }
494 
495  double zeta = der_press_ent_p(ent) / der_press_nbar_p(ent) ;
496 
497  return zeta ;
498 
499  }
500  else
501 
502  return 1./(der_press_nbar_p(hmin)-1.) ; // to ensure continuity
503 
504 }
505 
506 
507 // dln(e)/ln(H) from enthalpy
508 //---------------------------
509 
510 double Eos_tabul::der_ener_ent_p(double ent, const Param* ) const {
511 
512  if ( ent > hmin ) {
513 
514  if (ent > hmax) {
515  cout << "Eos_tabul::der_ener_ent_p : ent > hmax !" << endl ;
516  abort() ;
517  }
518 
519  return ( der_nbar_ent_p(ent)
520  *( double(1.) + press_ent_p(ent)/ener_ent_p(ent) )) ;
521 
522  } else return der_nbar_ent_p(hmin) ;
523 
524 }
525 
526 
527 // dln(p)/ln(H) from enthalpy
528 //---------------------------
529 
530 double Eos_tabul::der_press_ent_p(double ent, const Param* ) const {
531 
532  static int i_near = logh->get_taille() / 2 ;
533 
534  if ( ent > hmin ) {
535  if (ent > hmax) {
536  cout << "Eos_tabul::der_press_ent_p : ent > hmax !" << endl ;
537  abort() ;
538  }
539 
540  double logh0 = log10( ent ) ;
541  double logp0 ;
542  double dlpsdlh0 ;
543  interpol_herm(*logh, *logp, *dlpsdlh, logh0, i_near, logp0,
544  dlpsdlh0) ;
545 
546  return dlpsdlh0 ;
547 
548  }
549  else{
550 
551  return 0 ;
552  // return der_press_ent_p(hmin) ; // to ensure continuity
553  }
554 }
555 
556 
557 // dln(p)/dln(n) from enthalpy
558 //---------------------------
559 
560 double Eos_tabul::der_press_nbar_p(double ent, const Param*) const {
561 
562  static int i_near = logh->get_taille() / 2 ;
563 
564  if ( ent > hmin ) {
565  if (ent > hmax) {
566  cout << "Eos_tabul::der_press_nbar_p : ent > hmax !" << endl ;
567  abort() ;
568  }
569 
570  double logh0 = log10( ent ) ;
571  double dlpsdlnb0 ;
572 
573  interpol_linear(*logh, *dlpsdlnb, logh0, i_near, dlpsdlnb0) ;
574 
575 // cout << "gamma: " << dlpsdlnb0 << " ent: " << ent << endl;
576 
577  return dlpsdlnb0 ;
578 
579  }
580  else{
581 
582  return (*dlpsdlnb)(0) ;
583  // return der_press_nbar_p(hmin) ; // to ensure continuity
584 
585  }
586 
587 
588 }
589 }
Cmp log(const Cmp &)
Neperian logarithm.
Definition: cmp_math.C:296
Cmp exp(const Cmp &)
Exponential.
Definition: cmp_math.C:270
Lorene prototypes.
Definition: app_hor.h:64
Standard units of space, time and mass.
Equation of state base class.
Definition: eos.h:190
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:281
Tbl * logp
Table of .
Definition: eos_tabul.h:170
double hmin
Lower boundary of the enthalpy interval.
Definition: eos_tabul.h:161
Tbl * dlpsdlnb
Table of .
Definition: eos_tabul.h:179
Tbl * dlpsdlh
Table of .
Definition: eos_tabul.h:173
virtual void sauve(FILE *) const
Save in a file.
Definition: eos.C:179
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:361
virtual ~Eos_tabul()
Destructor.
Definition: eos_tabul.C:201
virtual double der_ener_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
Definition: eos_tabul.C:510
string tablename
Name of the file containing the tabulated data.
Definition: eos_tabul.h:156
string authors
Authors - reference for the table.
Definition: eos_tabul.h:158
virtual void sauve(FILE *) const
Save in a file.
Definition: eos_tabul.C:213
virtual double der_press_nbar_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
Definition: eos_tabul.C:560
Tbl * lognb
Table of .
Definition: eos_tabul.h:176
double hmax
Upper boundary of the enthalpy interval.
Definition: eos_tabul.h:164
virtual double der_nbar_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
Definition: eos_tabul.C:487
Parameter storage.
Definition: param.h:125
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
void read_table()
Reads the file containing the table and initializes in the arrays logh , logp and dlpsdlh ...
Definition: eos_tabul.C:226
Eos_tabul(const char *name_i, const char *table, const char *path)
Standard constructor.
Definition: eos_tabul.C:133
Cmp log10(const Cmp &)
Basis 10 logarithm.
Definition: cmp_math.C:322
int get_taille() const
Gives the total size (ie dim.taille)
Definition: tbl.h:397
virtual double der_press_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
Definition: eos_tabul.C:530
Basic array class.
Definition: tbl.h:161
virtual double press_ent_p(double ent, const Param *par=0x0) const
Computes the pressure from the log-enthalpy.
Definition: eos_tabul.C:462
Tbl * logh
Table of .
Definition: eos_tabul.h:167
virtual double ener_ent_p(double ent, const Param *par=0x0) const
Computes the total energy density from the log-enthalpy.
Definition: eos_tabul.C:434
virtual double nbar_ent_p(double ent, const Param *par=0x0) const
Computes the baryon density from the log-enthalpy.
Definition: eos_tabul.C:406