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