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
30char 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
121namespace Lorene {
122void interpol_herm(const Tbl& , const Tbl&, const Tbl&, double, int&,
123 double&, double& ) ;
124
125void interpol_linear(const Tbl&, const Tbl&, double, int&, double&) ;
126
127 //----------------------------//
128 // Constructors //
129 //----------------------------//
130
131// Standard constructor
132// --------------------
133Eos_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
213void 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() ;
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
406double 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
434double 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
462double 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
487double 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
510double 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
530double 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
560double 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}
virtual double press_ent_p(double ent, const Param *par=0x0) const
Computes the pressure from the log-enthalpy.
Definition eos_tabul.C:462
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 der_press_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
Definition eos_tabul.C:530
void read_table()
Reads the file containing the table and initializes in the arrays logh , logp and dlpsdlh .
Definition eos_tabul.C:226
virtual ~Eos_tabul()
Destructor.
Definition eos_tabul.C:201
Tbl * lognb
Table of .
Definition eos_tabul.h:176
Eos_tabul(const char *name_i, const char *table, const char *path)
Standard constructor.
Definition eos_tabul.C:133
Tbl * dlpsdlh
Table of .
Definition eos_tabul.h:173
Tbl * dlpsdlnb
Table of .
Definition eos_tabul.h:179
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
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
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
virtual void sauve(FILE *) const
Save in a file.
Definition eos_tabul.C:213
Tbl * logh
Table of .
Definition eos_tabul.h:167
string tablename
Name of the file containing the tabulated data.
Definition eos_tabul.h:156
double hmin
Lower boundary of the enthalpy interval.
Definition eos_tabul.h:161
Tbl * logp
Table of .
Definition eos_tabul.h:170
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
string authors
Authors - reference for the table.
Definition eos_tabul.h:158
double hmax
Upper boundary of the enthalpy interval.
Definition eos_tabul.h:164
Equation of state base class.
Definition eos.h:190
virtual void sauve(FILE *) const
Save in a file.
Definition eos.C:179
Parameter storage.
Definition param.h:125
Basic array class.
Definition tbl.h:161
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:361
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition tbl.h:281
int get_taille() const
Gives the total size (ie dim.taille)
Definition tbl.h:397
Cmp log10(const Cmp &)
Basis 10 logarithm.
Definition cmp_math.C:322
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
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.