107 string file_nb = files +
".nb" ;
108 string file_thermo = files +
".thermo" ;
110 ifstream in_nb(file_nb.data()) ;
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 ;
122 in_nb >> index1 >> index2 ;
123 int nbp = index2 - index1 + 1 ;
126 press =
new double[nbp] ;
127 nb =
new double[nbp] ;
128 ro =
new double[nbp] ;
144 double nb_fm3, rho_cgs, p_cgs, p_over_nb_comp, eps_comp ;
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 ;
152 ifstream in_p_rho (file_thermo.data()) ;
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 ;
160 in_p_rho >> m_neutron_MeV >> m_proton_MeV ;
161 in_p_rho.ignore(1000,
'\n') ;
163 double p_convert = mev_si * 1.e45 * 10. ;
164 double eps_convert = mev_si * 1.e42 / (c_si*c_si) ;
168 for (
int i=0; i<nbp; i++) {
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 ;
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 ;
185 press[i] = p_cgs / c2_cgs ;
191 for (
int i=0; i<nbp; i++) {
192 double h =
log( (ro[i] + press[i]) /
193 (10 * nb[i] * rhonuc_cgs) ) ;
195 if (i==0) { ww = h ; }
196 h = h - ww + 1.e-14 ;
200 dlpsdlh->
set(i) = h * (ro[i] + press[i]) / press[i] ;
206 double p0, p1, p2, n0, n1, n2, dpdnb;
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) ;
223 for(
int i=1;i<nbp-1;i++) {
225 p0 =
log(press[i-1]);
227 p2 =
log(press[i+1]);
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) ;
242 p0 =
log(press[nbp-3]);
243 p1 =
log(press[nbp-2]);
244 p2 =
log(press[nbp-1]);
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) ;