LORENE
eos_compose.C
1/*
2 * Methods of class Eos_CompOSE
3 *
4 * (see file eos_compose.h for documentation).
5 *
6 */
7
8/*
9 * Copyright (c) 2010, 2014 Jerome Novak
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_compose_C[] = "$Header: /cvsroot/Lorene/C++/Source/Eos/eos_compose.C,v 1.7 2015/12/04 16:27:05 j_novak Exp $ " ;
31
32/*
33 * $Id: eos_compose.C,v 1.7 2015/12/04 16:27:05 j_novak Exp $
34 * $Log: eos_compose.C,v $
35 * Revision 1.7 2015/12/04 16:27:05 j_novak
36 * Correction of constructor calling.
37 *
38 * Revision 1.6 2015/08/04 14:41:29 j_novak
39 * Back to previous version for Eos_CompOSE. Enthalpy-consistent EoS can be accessed using Eos_consistent class (derived from Eos_CompOSE).
40 *
41 * Revision 1.5 2015/01/27 14:22:38 j_novak
42 * New methods in Eos_tabul to correct for EoS themro consistency (optional).
43 *
44 * Revision 1.4 2014/10/13 08:52:52 j_novak
45 * Lorene classes and functions now belong to the namespace Lorene.
46 *
47 * Revision 1.3 2014/07/01 09:26:21 j_novak
48 * Improvement of comments
49 *
50 * Revision 1.2 2014/06/30 16:13:18 j_novak
51 * New methods for reading directly from CompOSE files.
52 *
53 * Revision 1.1 2014/03/06 15:53:34 j_novak
54 * Eos_compstar is now Eos_compOSE. Eos_tabul uses strings and contains informations about authors.
55 *
56 * Revision 1.1 2010/02/03 14:56:45 j_novak
57 * *** empty log message ***
58 *
59 *
60 * $Header: /cvsroot/Lorene/C++/Source/Eos/eos_compose.C,v 1.7 2015/12/04 16:27:05 j_novak Exp $
61 *
62 */
63
64#include <string>
65
66// Headers Lorene
67#include "headcpp.h"
68#include "eos.h"
69#include "tbl.h"
70#include "utilitaires.h"
71#include "unites.h"
72
73 //----------------------------//
74 // Constructors //
75 //----------------------------//
76
77// Standard constructor
78// --------------------
79namespace Lorene {
80 Eos_CompOSE::Eos_CompOSE(const char* file_name)
81 : Eos_tabul("Tabulated EoS", file_name)
82{}
83
84
85// Constructor from binary file
86// ----------------------------
88{}
89
90
91// Constructor from a formatted file
92// ---------------------------------
93 Eos_CompOSE::Eos_CompOSE(ifstream& fich) : Eos_tabul(fich)
94{}
95
96
97// Constructor from CompOSE data files
98// ------------------------------------
99 Eos_CompOSE::Eos_CompOSE(const string& files) : Eos_tabul("CompOSE Eos")
100{
101
102 using namespace Unites ;
103
104 // Files containing data and a test
105 //---------------------------------
106 tablename = files ;
107 string file_nb = files + ".nb" ;
108 string file_thermo = files + ".thermo" ;
109
110 ifstream in_nb(file_nb.data()) ;
111 if (!in_nb) {
112 cerr << "Eos_CompOSE::Eos_CompOSE(string): " << endl ;
113 cerr << "Problem in opening the EOS file!" << endl ;
114 cerr << "While trying to open " << file_nb << endl ;
115 cerr << "Aborting..." << endl ;
116 abort() ;
117 }
118
119 // obtaining the size of the tables for memory allocation
120 //-------------------------------------------------------
121 int index1, index2 ;
122 in_nb >> index1 >> index2 ;
123 int nbp = index2 - index1 + 1 ;
124 assert(nbp > 0) ;
125
126 press = new double[nbp] ;
127 nb = new double[nbp] ;
128 ro = new double[nbp] ;
129
130 logh = new Tbl(nbp) ;
131 logp = new Tbl(nbp) ;
132 dlpsdlh = new Tbl(nbp) ;
133 lognb = new Tbl(nbp) ;
134 dlpsdlnb = new Tbl(nbp) ;
135
136 logh->set_etat_qcq() ;
137 logp->set_etat_qcq() ;
141
142 // Variables and conversion
143 //-------------------------
144 double nb_fm3, rho_cgs, p_cgs, p_over_nb_comp, eps_comp ;
145 double dummy_x ;
146 int dummy_n ;
147
148 double rhonuc_cgs = rhonuc_si * 1e-3 ;
149 double c2_cgs = c_si * c_si * 1e4 ;
150 double m_neutron_MeV, m_proton_MeV ;
151
152 ifstream in_p_rho (file_thermo.data()) ;
153 if (!in_p_rho) {
154 cerr << "Eos_CompOSE::Eos_CompOSE(string): " << endl ;
155 cerr << "Problem in opening the EOS file!" << endl ;
156 cerr << "While trying to open " << file_thermo << endl ;
157 cerr << "Aborting..." << endl ;
158 abort() ;
159 }
160 in_p_rho >> m_neutron_MeV >> m_proton_MeV ; //Neutron and proton masses
161 in_p_rho.ignore(1000, '\n') ;
162
163 double p_convert = mev_si * 1.e45 * 10. ; // Conversion from MeV/fm^3 to cgs
164 double eps_convert = mev_si * 1.e42 / (c_si*c_si) ; //From meV/fm^3 to g/cm^3
165
166 // Main loop reading the table
167 //----------------------------
168 for (int i=0; i<nbp; i++) {
169 in_nb >> nb_fm3 ;
170 in_p_rho >> dummy_n >> dummy_n >> dummy_n >> p_over_nb_comp ;
171 in_p_rho >> dummy_x >> dummy_x >> dummy_x >> dummy_x >> dummy_x >> eps_comp ;
172 in_p_rho.ignore(1000, '\n') ;
173 p_cgs = p_over_nb_comp * nb_fm3 * p_convert ;
174 rho_cgs = ( eps_comp + 1. ) * m_neutron_MeV * nb_fm3 * eps_convert ;
175
176 if ( (nb_fm3<0) || (rho_cgs<0) || (p_cgs < 0) ){
177 cout << "Eos_CompOSE::Eos_CompOSE(string): " << endl ;
178 cout << "Negative value in table!" << endl ;
179 cout << "nb = " << nb_fm3 << ", rho = " << rho_cgs <<
180 ", p = " << p_cgs << endl ;
181 cout << "Aborting..." << endl ;
182 abort() ;
183 }
184
185 press[i] = p_cgs / c2_cgs ;
186 nb[i] = nb_fm3 ;
187 ro[i] = rho_cgs ;
188 }
189
190 double ww = 0. ;
191 for (int i=0; i<nbp; i++) {
192 double h = log( (ro[i] + press[i]) /
193 (10 * nb[i] * rhonuc_cgs) ) ;
194
195 if (i==0) { ww = h ; }
196 h = h - ww + 1.e-14 ;
197
198 logh->set(i) = log10( h ) ;
199 logp->set(i) = log10( press[i] / rhonuc_cgs ) ;
200 dlpsdlh->set(i) = h * (ro[i] + press[i]) / press[i] ;
201 lognb->set(i) = log10(nb[i]) ;
202 }
203
204 // Computation of dpdnb
205 //---------------------
206 double p0, p1, p2, n0, n1, n2, dpdnb;
207
208 // special case: i=0
209 p0 = log(press[0]);
210 p1 = log(press[1]);
211 p2 = log(press[2]);
212
213 n0 = log(nb[0]);
214 n1 = log(nb[1]);
215 n2 = log(nb[2]);
216
217 dpdnb = p0*(2*n0-n1-n2)/(n0-n1)/(n0-n2) +
218 p1*(n0-n2)/(n1-n0)/(n1-n2) +
219 p2*(n0-n1)/(n2-n0)/(n2-n1) ;
220
221 dlpsdlnb->set(0) = dpdnb ;
222
223 for(int i=1;i<nbp-1;i++) {
224
225 p0 = log(press[i-1]);
226 p1 = log(press[i]);
227 p2 = log(press[i+1]);
228
229 n0 = log(nb[i-1]);
230 n1 = log(nb[i]);
231 n2 = log(nb[i+1]);
232
233 dpdnb = p0*(n1-n2)/(n0-n1)/(n0-n2) +
234 p1*(2*n1-n0-n2)/(n1-n0)/(n1-n2) +
235 p2*(n1-n0)/(n2-n0)/(n2-n1) ;
236
237 dlpsdlnb->set(i) = dpdnb ;
238
239 }
240
241 // special case: i=nbp-1
242 p0 = log(press[nbp-3]);
243 p1 = log(press[nbp-2]);
244 p2 = log(press[nbp-1]);
245
246 n0 = log(nb[nbp-3]);
247 n1 = log(nb[nbp-2]);
248 n2 = log(nb[nbp-1]);
249
250 dpdnb = p0*(n2-n1)/(n0-n1)/(n0-n2) +
251 p1*(n2-n0)/(n1-n0)/(n1-n2) +
252 p2*(2*n2-n0-n1)/(n2-n0)/(n2-n1) ;
253
254 dlpsdlnb->set(nbp-1) = dpdnb ;
255
256 hmin = pow( double(10), (*logh)(0) ) ;
257 hmax = pow( double(10), (*logh)(nbp-1) ) ;
258
259 // Cleaning
260 //---------
261 delete [] press ;
262 delete [] nb ;
263 delete [] ro ;
264
265
266}
267
268
269
270 //--------------//
271 // Destructor //
272 //--------------//
273
275
276 // does nothing
277
278}
279
280
281 //------------------------//
282 // Comparison operators //
283 //------------------------//
284
285
286bool Eos_CompOSE::operator==(const Eos& eos_i) const {
287
288 bool resu = true ;
289
290 if ( eos_i.identify() != identify() ) {
291 cout << "The second EOS is not of type Eos_CompOSE !" << endl ;
292 resu = false ;
293 }
294
295 return resu ;
296
297}
298
299bool Eos_CompOSE::operator!=(const Eos& eos_i) const {
300
301 return !(operator==(eos_i)) ;
302
303}
304
305 //------------//
306 // Outputs //
307 //------------//
308
309
310ostream& Eos_CompOSE::operator>>(ostream & ost) const {
311
312 ost << "EOS of class Eos_CompOSE." << endl ;
313 ost << "Built from file " << tablename << endl ;
314 ost << "Authors : " << authors << endl ;
315 ost << "Number of points in file : " << logh->get_taille() << endl ;
316 return ost ;
317
318}
319
320
321}
virtual ~Eos_CompOSE()
Destructor.
virtual bool operator==(const Eos &) const
Comparison operator (egality)
virtual bool operator!=(const Eos &) const
Comparison operator (difference)
Eos_CompOSE(const string &files_path)
Constructor from CompOSE data.
Definition eos_compose.C:99
virtual ostream & operator>>(ostream &) const
Operator >>
virtual int identify() const
Returns a number to identify the sub-classe of Eos the object belongs to.
Base class for tabulated equations of state.
Definition eos_tabul.h:149
Tbl * lognb
Table of .
Definition eos_tabul.h:176
Tbl * dlpsdlh
Table of .
Definition eos_tabul.h:173
Tbl * dlpsdlnb
Table of .
Definition eos_tabul.h:179
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
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 int identify() const =0
Returns a number to identify the sub-classe of Eos the object belongs to.
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 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.