hoteos_tabul.C
/*
* Methods of class Hoteos_tabul
*
* (see file hoteos.h for documentation).
*
*/
/*
* Copyright (c) 2015 Jerome Novak
*
* This file is part of LORENE.
*
* LORENE is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* LORENE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with LORENE; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*/
char hoteos_tabul_C[] = "$Header: /cvsroot/Lorene/C++/Source/Eos/hoteos_tabul.C,v 1.2 2015/12/08 15:42:17 j_novak Exp $" ;
/*
* $Id: hoteos_tabul.C,v 1.2 2015/12/08 15:42:17 j_novak Exp $
* $Log: hoteos_tabul.C,v $
* Revision 1.2 2015/12/08 15:42:17 j_novak
* Low/zero entropy is set to the lowest value in the table in computational functions.
*
* Revision 1.1 2015/12/08 10:52:18 j_novak
* New class Hoteos_tabul for tabulated temperature-dependent EoSs.
*
*
*
* $Header: /cvsroot/Lorene/C++/Source/Eos/hoteos_tabul.C,v 1.2 2015/12/08 15:42:17 j_novak Exp $
*
*/
// headers C
#include <cstdlib>
#include <cstring>
#include <cmath>
// Headers Lorene
#include "hoteos.h"
#include "eos.h"
#include "utilitaires.h"
#include "unites.h"
namespace Lorene {
void interpol_herm_2d(const Tbl&, const Tbl&, const Tbl&, const Tbl&, const Tbl&,
const Tbl&, double, double, double&, double&, double&) ;
//----------------------------//
// Constructors //
//----------------------------//
// Standard constructor
// --------------------
Hoteos_tabul::Hoteos_tabul(const string& filename):
Hot_eos("Tabulated hot EoS"), tablename(filename), authors("Unknown"),
hmin(0.), hmax(1.), sbmin(0.), sbmax(1.)
{
set_arrays_0x0() ;
read_table() ;
}
// Constructor from binary file
// ----------------------------
Hoteos_tabul::Hoteos_tabul(FILE* fich) : Hot_eos(fich) {
char tmp_string[160] ;
fread(tmp_string, sizeof(char), 160, fich) ;
tablename = tmp_string ;
set_arrays_0x0() ;
read_table() ;
}
// Constructor from a formatted file
// ---------------------------------
Hoteos_tabul::Hoteos_tabul(ifstream& fich) :
Hot_eos(fich) {
fich >> tablename ;
set_arrays_0x0() ;
read_table() ;
}
//Sets the arrays to the null pointer
void Hoteos_tabul::set_arrays_0x0() {
hhh = 0x0 ;
s_B = 0x0 ;
ppp = 0x0 ;
dpdh = 0x0 ;
dpds = 0x0 ;
d2p = 0x0 ;
}
//--------------//
// Destructor //
//--------------//
Hoteos_tabul::~Hoteos_tabul(){
if (hhh != 0x0) delete hhh ;
if (s_B != 0x0) delete s_B ;
if (ppp != 0x0) delete ppp ;
if (dpdh != 0x0) delete dpdh ;
if (dpds != 0x0) delete dpds ;
if (d2p != 0x0) delete d2p ;
}
//------------//
// Outputs //
//------------//
void Hoteos_tabul::sauve(FILE* fich) const {
Hot_eos::sauve(fich) ;
char tmp_string[160] ;
strcpy(tmp_string, tablename.c_str()) ;
fwrite(tmp_string, sizeof(char), 160, fich) ;
}
ostream& Hoteos_tabul::operator>>(ostream & ost) const {
ost << "Hot EOS of class Hoteos_tabul (tabulated hot beta-equilibrium EoS) : "
<< endl ;
ost << "Built from file " << tablename << endl ;
ost << "Authors : " << authors << endl ;
ost << "Number of points in file : " << hhh->get_dim(0)
<< " in enthalpy, and " << hhh->get_dim(1)
<< " in entropy." << endl ;
return ost ;
}
//------------------------//
// Reading of the table //
//------------------------//
void Hoteos_tabul::read_table() {
using namespace Unites ;
ifstream fich(tablename.data()) ;
if (!fich) {
cerr << "Hoteos_tabul::read_table(): " << endl ;
cerr << "Problem in opening the EOS file!" << endl ;
cerr << "While trying to open " << tablename << endl ;
cerr << "Aborting..." << endl ;
abort() ;
}
fich.ignore(1000, '\n') ; // Jump over the first header
fich.ignore(1) ;
getline(fich, authors, '\n') ;
for (int i=0; i<3; i++) { // jump over the file
fich.ignore(1000, '\n') ; // header
} //
int nbp1, nbp2 ;
fich >> nbp1 >> nbp2 ; fich.ignore(1000, '\n') ; // number of data
if ( (nbp1<=0) || (nbp2<=0) ) {
cerr << "Hoteos_tabul::read_table(): " << endl ;
cerr << "Wrong value for the number of lines!" << endl ;
cerr << "nbp1 = " << nbp1 << ", nbp2 = " << nbp2 << endl ;
cerr << "Aborting..." << endl ;
abort() ;
}
for (int i=0; i<3; i++) { // jump over the table
fich.ignore(1000, '\n') ; // description
}
ppp = new Tbl(nbp2, nbp1) ;
hhh = new Tbl(nbp2, nbp1) ;
s_B = new Tbl(nbp2, nbp1) ;
dpdh = new Tbl(nbp2, nbp1) ;
dpds = new Tbl(nbp2, nbp1) ;
d2p = new Tbl(nbp2, nbp1) ;
ppp->set_etat_qcq() ;
hhh->set_etat_qcq() ;
s_B->set_etat_qcq() ;
dpdh->set_etat_qcq() ;
dpds->set_etat_qcq() ;
d2p->set_etat_qcq() ;
double c2 = c_si * c_si ;
double dummy, nb_fm3, rho_cgs, p_cgs, mu_MeV, entr, temp, der2 ;
double ww = 0. ;
for (int j=0; j<nbp2; j++) {
for (int i=0; i<nbp1; i++) {
fich >> mu_MeV >> entr >> nb_fm3 >> temp >> p_cgs >> der2
>> dummy >> rho_cgs ;
fich.ignore(1000,'\n') ;
if ( (nb_fm3<0) || (rho_cgs<0) || (p_cgs < 0) ){
cerr << "Eos_mag::read_table(): " << endl ;
cerr << "Negative value in table!" << endl ;
cerr << "nb = " << nb_fm3 << ", rho = " << rho_cgs <<
", p = " << p_cgs << endl ;
cerr << "Aborting..." << endl ;
abort() ;
}
double psc2 = 0.1*p_cgs/c2 ; // in kg m^-3
double rho_si = rho_cgs*1000. ; // in kg m^-3
double h_read = log(mu_MeV) ;
if ( (i==0) && (j==0) ) ww = h_read ;
double h_new = h_read - ww ;
ppp->set(j, i) = psc2/rhonuc_si ;
hhh->set(j, i) = h_new ;
s_B->set(j, i) = entr ; // in Lorene units (k_B)
dpdh->set(j, i) = (rho_si + psc2)/rhonuc_si ;
dpds->set(j, i) = -temp*nb_fm3*mevpfm3 ;
d2p->set(j, i) = 0.1*der2*mu_MeV/(c2*rhonuc_si) ;
}
}
hmin = (*hhh)(0, 0) ;
hmax = (*hhh)(0, nbp1-1) ;
sbmin = (*s_B)(0, 0) ;
sbmax = (*s_B)(nbp2-1, 0) ;
cout << "hmin: " << hmin << ", hmax: " << hmax << endl ;
cout << "sbmin: " << sbmin << ", sbmax: " << sbmax << endl ;
fich.close();
}
//-------------------------------//
// The corresponding cold Eos //
//-------------------------------//
const Eos& Hoteos_tabul::new_cold_Eos() const {
if (p_cold_eos == 0x0) {
cerr << "Warning: Hoteos_tabul::new_cold_Eos " <<
"The corresponding cold EoS is likely not to function." << endl ;
p_cold_eos = new Eos_CompOSE(tablename.c_str()) ;
}
return *p_cold_eos ;
}
//------------------------//
// Comparison operators //
//------------------------//
bool Hoteos_tabul::operator==(const Hot_eos& eos_i) const {
bool resu = true ;
if ( eos_i.identify() != identify() ) {
cout << "The second EOS is not of type Hoteos_tabul !" << endl ;
resu = false ;
}
return resu ;
}
bool Hoteos_tabul::operator!=(const Hot_eos& eos_i) const {
return !(operator==(eos_i)) ;
}
//------------------------------//
// Computational routines //
//------------------------------//
// Baryon density from enthalpy and entropy
//------------------------------------------
double Hoteos_tabul::nbar_Hs_p(double ent, double sb) const {
if ((ent > hmin - 1.e-12) && (ent < hmin))
ent = hmin ;
if (sb < sbmin) sb = sbmin ;
if ( ent >= hmin ) {
if (ent > hmax) {
cout << "Hoteos_tabul::nbar_Hs_p : ent > hmax !" << endl ;
abort() ;
}
if (sb > sbmax) {
cerr << "Hoteos_tabul::nbar_Hs_p : s_B not in the tabulated interval !"
<< endl ;
cerr << "s_B = " << sb << ", sbmin = " << sbmin << ", sbmax = " << sbmax
<< endl ;
abort() ;
}
double p_int, dpds_int, dpdh_int ;
interpol_herm_2d(*s_B, *hhh, *ppp, *dpds, *dpdh, *d2p, sb, ent, p_int,
dpds_int, dpdh_int) ;
double nbar_int = dpdh_int * exp(-ent) ;
return nbar_int ;
}
else{
return 0 ;
}
}
// Energy density from enthalpy and entropy
//-----------------------------------------
double Hoteos_tabul::ener_Hs_p(double ent, double sb) const {
if ((ent > hmin - 1.e-12) && (ent < hmin))
ent = hmin ;
if (sb < sbmin) sb = sbmin ;
if ( ent >= hmin ) {
if (ent > hmax) {
cout << "Hoteos_tabul::ener_Hs_p : ent > hmax !" << endl ;
abort() ;
}
if (sb > sbmax) {
cerr << "Hoteos_tabul::ener_Hs_p : s_B not in the tabulated interval !"
<< endl ;
cerr << "s_B = " << sb << ", sbmin = " << sbmin << ", sbmax = " << sbmax
<< endl ;
abort() ;
}
double p_int, dpds_int, dpdh_int ;
interpol_herm_2d(*s_B, *hhh, *ppp, *dpds, *dpdh, *d2p, sb, ent, p_int,
dpds_int, dpdh_int) ;
double nbar_int = dpdh_int * exp(-ent) ;
double f_int = - p_int + exp(ent) * nbar_int;
return f_int ;
}
else{
return 0 ;
}
}
// Pressure from enthalpy and entropy
//-----------------------------------
double Hoteos_tabul::press_Hs_p(double ent, double sb) const {
if ((ent > hmin - 1.e-12) && (ent < hmin))
ent = hmin ;
if (sb < sbmin) sb = sbmin ;
if ( ent >= hmin ) {
if (ent > hmax) {
cout << "Hoteos_tabul::press_Hs_p : ent > hmax !" << endl ;
abort() ;
}
if (sb > sbmax) {
cerr << "Hoteos_tabul::press_Hs_p : s_B not in the tabulated interval !"
<< endl ;
cerr << "s_B = " << sb << ", sbmin = " << sbmin << ", sbmax = " << sbmax
<< endl ;
abort() ;
}
double p_int, dpds_int, dpdh_int ;
interpol_herm_2d(*s_B, *hhh, *ppp, *dpds, *dpdh, *d2p, sb, ent, p_int,
dpds_int, dpdh_int) ;
return p_int ;
}
else{
return 0 ;
}
}
// Temperature from enthalpy and entropy
//--------------------------------------
double Hoteos_tabul::temp_Hs_p(double ent, double sb) const {
if ((ent > hmin - 1.e-12) && (ent < hmin))
ent = hmin ;
if (sb < sbmin) sb = sbmin ;
if ( ent >= hmin ) {
if (ent > hmax) {
cout << "Hoteos_tabul::temp_Hs_p : ent > hmax !" << endl ;
abort() ;
}
if (sb > sbmax) {
cerr << "Hoteos_tabul::temp_Hs_p : s_B not in the tabulated interval !"
<< endl ;
cerr << "s_B = " << sb << ", sbmin = " << sbmin << ", sbmax = " << sbmax
<< endl ;
abort() ;
}
double p_int, dpds_int, dpdh_int ;
interpol_herm_2d(*s_B, *hhh, *ppp, *dpds, *dpdh, *d2p, sb, ent, p_int,
dpds_int, dpdh_int) ;
double nbar_int = dpdh_int * exp(-ent) ;
double temp_int = -dpds_int / nbar_int ;
return temp_int ;
}
else {
return 0 ;
}
}
} //End of namespace Lorene
Lorene::Eos_CompOSE
Equation of state for the CompOSE database.
Definition eos_compose.h:77
Lorene::Eos
Equation of state base class.
Definition eos.h:190
Lorene::Hot_eos::identify
virtual int identify() const =0
Returns a number to identify the sub-classe of Hot_eos the object belongs to.
Lorene::Hot_eos::p_cold_eos
Eos * p_cold_eos
Corresponding cold Eos.
Definition hoteos.h:108
Lorene::Hot_eos::Hot_eos
Hot_eos()
Standard constructor.
Definition hoteos.C:56
Lorene::Hot_eos::sauve
virtual void sauve(FILE *) const
Save in a file.
Definition hoteos.C:129
Lorene::Hoteos_tabul::hmin
double hmin
Lower boundary of the enthalpy interval.
Definition hoteos.h:561
Lorene::Hoteos_tabul::sauve
virtual void sauve(FILE *) const
Save in a file.
Definition hoteos_tabul.C:127
Lorene::Hoteos_tabul::~Hoteos_tabul
virtual ~Hoteos_tabul()
Destructor.
Definition hoteos_tabul.C:114
Lorene::Hoteos_tabul::sbmax
double sbmax
Upper boundary of the entropy interval.
Definition hoteos.h:570
Lorene::Hoteos_tabul::nbar_Hs_p
virtual double nbar_Hs_p(double ent, double sb) const
Computes the baryon density from the log-enthalpy and entropy per baryon (virtual function implemente...
Definition hoteos_tabul.C:294
Lorene::Hoteos_tabul::Hoteos_tabul
Hoteos_tabul(const string &filename)
Standard constructor from a filename.
Definition hoteos_tabul.C:70
Lorene::Hoteos_tabul::operator!=
virtual bool operator!=(const Hot_eos &) const
Comparison operator (difference).
Definition hoteos_tabul.C:282
Lorene::Hoteos_tabul::hhh
Tbl * hhh
Table of .
Definition hoteos.h:573
Lorene::Hoteos_tabul::d2p
Tbl * d2p
Table of .
Definition hoteos.h:588
Lorene::Hoteos_tabul::set_arrays_0x0
void set_arrays_0x0()
Sets all the arrays to the null pointer.
Definition hoteos_tabul.C:101
Lorene::Hoteos_tabul::operator>>
virtual ostream & operator>>(ostream &) const
Operator >>.
Definition hoteos_tabul.C:136
Lorene::Hoteos_tabul::operator==
virtual bool operator==(const Hot_eos &) const
Comparison operator (egality).
Definition hoteos_tabul.C:269
Lorene::Hoteos_tabul::press_Hs_p
virtual double press_Hs_p(double ent, double sb) const
Computes the pressure from the log-enthalpy and entropy per baryon (virtual function implemented in t...
Definition hoteos_tabul.C:369
Lorene::Hoteos_tabul::hmax
double hmax
Upper boundary of the enthalpy interval.
Definition hoteos.h:564
Lorene::Hoteos_tabul::new_cold_Eos
virtual const Eos & new_cold_Eos() const
Returns the corresponding cold Eos.
Definition hoteos_tabul.C:252
Lorene::Hoteos_tabul::tablename
string tablename
Name of the file containing the tabulated data.
Definition hoteos.h:556
Lorene::Hoteos_tabul::ener_Hs_p
virtual double ener_Hs_p(double ent, double sb) const
Computes the total energy density from the log-enthalpy and entropy per baryon (virtual function impl...
Definition hoteos_tabul.C:330
Lorene::Hoteos_tabul::dpdh
Tbl * dpdh
Table of .
Definition hoteos.h:582
Lorene::Hoteos_tabul::temp_Hs_p
virtual double temp_Hs_p(double ent, double sb) const
Computes the temperature from the log-enthalpy and entropy per baryon (virtual function implemented i...
Definition hoteos_tabul.C:403
Lorene::Hoteos_tabul::identify
virtual int identify() const
Returns a number to identify the sub-classe of Hot_eos the object belongs to.
Definition hoteos_from_file.C:60
Lorene::Hoteos_tabul::dpds
Tbl * dpds
Table of .
Definition hoteos.h:585
Lorene::Hoteos_tabul::authors
string authors
Authors - reference for the table.
Definition hoteos.h:558
Lorene::Hoteos_tabul::s_B
Tbl * s_B
Table of , entropy per baryon (in units of Boltzmann constant).
Definition hoteos.h:576
Lorene::Hoteos_tabul::read_table
void read_table()
Reads the file containing the table and initializes in the arrays hhh , s_B, ppp, ....
Definition hoteos_tabul.C:153
Lorene::Hoteos_tabul::ppp
Tbl * ppp
Table of pressure $P$.
Definition hoteos.h:579
Lorene::Hoteos_tabul::sbmin
double sbmin
Lower boundary of the entropy interval.
Definition hoteos.h:567
Lorene::Tbl
Basic array class.
Definition tbl.h:161
Lorene::exp
Cmp exp(const Cmp &)
Exponential.
Definition cmp_math.C:270
Lorene::log
Cmp log(const Cmp &)
Neperian logarithm.
Definition cmp_math.C:296
Lorene
Lorene prototypes.
Definition app_hor.h:64
Unites
Standard units of space, time and mass.