eos_mag.C
/* Methods of class Eos_mag
*
* (see file eos_mag.h for documentation).
*
*/
/*
* Copyright (c) 2011 Thomas Elghozi & Jerome Novak
* Copyright (c) 2013 Debarati Chatterjee
*
* 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 eos_mag_C[] = "$Header: /cvsroot/Lorene/C++/Source/Eos/eos_mag.C,v 1.13 2014/10/13 08:52:53 j_novak Exp $" ;
/*
* $Id: eos_mag.C,v 1.13 2014/10/13 08:52:53 j_novak Exp $
* $Log: eos_mag.C,v $
* Revision 1.13 2014/10/13 08:52:53 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.12 2014/09/22 16:13:10 j_novak
* Minor modif.
*
* Revision 1.11 2014/05/27 12:32:28 j_novak
* Added possibility to converge to a given magnetic moment.
*
* Revision 1.10 2014/05/13 15:37:12 j_novak
* Updated to new magnetic units.
*
* Revision 1.9 2014/04/28 12:48:13 j_novak
* Minor modifications.
*
* Revision 1.8 2014/03/11 14:27:26 j_novak
* Corrected a missing 4pi term.
*
* Revision 1.7 2014/03/03 16:23:08 j_novak
* Updated error message
*
* Revision 1.6 2013/12/12 16:07:30 j_novak
* interpol_herm_2d outputs df/dx, used to get the magnetization.
*
* Revision 1.5 2013/11/25 15:00:52 j_novak
* Correction of memory error.
*
* Revision 1.4 2013/11/14 16:12:55 j_novak
* Corrected a mistake in the units.
*
* Revision 1.2 2011/10/04 16:05:19 j_novak
* Update of Eos_mag class. Suppression of loge, re-definition of the derivatives
* and use of interpol_herm_2d.
*
* Revision 1.1 2011/06/16 10:49:18 j_novak
* New class Eos_mag for EOSs depending on density and magnetic field.
*
*
* $Header: /cvsroot/Lorene/C++/Source/Eos/eos_mag.C,v 1.13 2014/10/13 08:52:53 j_novak Exp $
*
*/
// headers C
#include <cmath>
// Headers Lorene
#include "eos.h"
#include "cmp.h"
#include "param.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
// --------------------
Eos_mag::Eos_mag(const char* name_i, const char* table,
const char* path) : Eos(name_i), tablename(path) {
tablename += '/' ;
tablename += table ;
read_table() ;
}
// Standard constructor with full filename
// ---------------------------------------
Eos_mag::Eos_mag(const char* name_i, const char* file_name)
: Eos(name_i), tablename(file_name) {
read_table() ;
}
// Constructor from binary file
// ----------------------------
Eos_mag::Eos_mag(FILE* fich) : Eos(fich) {
char tmp_name[160] ;
fread(tmp_name, sizeof(char), 160, fich) ;
tablename += tmp_name ;
read_table() ;
}
// Constructor from a formatted file
// ---------------------------------
Eos_mag::Eos_mag(ifstream& fich, const char* table) : Eos(fich) {
fich >> tablename ;
tablename += '/' ;
tablename += table ;
read_table() ;
}
Eos_mag::Eos_mag(ifstream& fich) : Eos(fich) {
fich >> tablename ;
read_table() ;
}
//--------------//
// Destructor //
//--------------//
Eos_mag::~Eos_mag(){
delete d2lp ;
delete dlpsdB ;
delete dlpsdlh ;
delete Bfield ;
delete logh ;
delete logp ;
}
//------------//
// Outputs //
//------------//
void Eos_mag::sauve(FILE* fich) const {
Eos::sauve(fich) ;
fwrite(tablename.c_str(), sizeof(char), tablename.size(), fich) ;
}
//------------------------//
// Comparison operators //
//------------------------//
bool Eos_mag::operator==(const Eos& eos_i) const {
bool resu = true ;
if ( eos_i.identify() != identify() ) {
cerr << "The second EOS is not of type Eos_mag !" << endl ;
resu = false ;
}
return resu ;
}
bool Eos_mag::operator!=(const Eos& eos_i) const {
return !(operator==(eos_i)) ;
}
//------------//
// Outputs //
//------------//
ostream& Eos_mag::operator>>(ostream & ost) const {
ost <<
"EOS of class Eos_mag : tabulated EOS depending on two variables: enthalpy and magnetic field."
<< '\n' ;
ost << "Table read from " << tablename << endl ;
return ost ;
}
//------------------------//
// Reading of the table //
//------------------------//
void Eos_mag::read_table() {
using namespace Unites_mag ;
ifstream fich(tablename.data()) ;
if (!fich) {
cerr << "Eos_mag::read_table(): " << endl ;
cerr << "Problem in opening the EOS file!" << endl ;
cerr << "Aborting..." << endl ;
abort() ;
}
for (int i=0; i<5; 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 << "Eos_mag::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') ;
}
logp = new Tbl(nbp2, nbp1) ;
logh = new Tbl(nbp2, nbp1) ;
Bfield = new Tbl(nbp2, nbp1) ;
dlpsdlh = new Tbl(nbp2, nbp1) ;
dlpsdB = new Tbl(nbp2, nbp1) ;
d2lp = new Tbl(nbp2, nbp1) ;
logp->set_etat_qcq() ;
logh->set_etat_qcq() ;
Bfield->set_etat_qcq() ;
dlpsdlh->set_etat_qcq() ;
dlpsdB->set_etat_qcq() ;
d2lp->set_etat_qcq() ;
double c2 = c_si * c_si ;
double B_unit = mag_unit / 1.e11 ;
double M_unit = B_unit*mu0/(4*M_PI) ;
int no1 ;
double nb_fm3, rho_cgs, p_cgs, mu_MeV, magB_PG, magM_PG, chi_PGpMeV ;
double ww = 0. ;
for (int j=0; j<nbp2; j++) {
for (int i=0; i<nbp1; i++) {
fich >> no1 >> nb_fm3 >> rho_cgs >> p_cgs >> mu_MeV
>> magB_PG >> magM_PG >> chi_PGpMeV ;
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 ;
logp->set(j, i) = psc2/rhonuc_si ;
logh->set(j, i) = h_new ;
Bfield->set(j, i) = magB_PG / B_unit ; // in Lorene units
dlpsdlh->set(j, i) = (rho_si + psc2)/rhonuc_si ;
dlpsdB->set(j, i) = magM_PG / M_unit ;
d2lp->set(j, i) = mu_MeV*chi_PGpMeV / (M_unit) ;
}
}
hmin = (*logh)(0, 0) ;
hmax = (*logh)(0, nbp1-1) ;
Bmax = (*Bfield)(nbp2-1, 0) ;
// cout << "hmin: " << hmin << ", hmax: " << hmax << endl ;
// cout << "Bmax: " << Bmax << endl ;
fich.close();
}
//------------------------------//
// Computational routines //
//------------------------------//
// Baryon density from enthalpy
//------------------------------
double Eos_mag::nbar_ent_p(double ent, const Param* par ) const {
using namespace Unites_mag ;
if ((ent > hmin - 1.e-12) && (ent < hmin))
ent = hmin ;
if ( ent >= hmin) {
if (ent > hmax) {
cerr << "Eos_tabul::nbar_ent_p : ent > hmax !" << endl ;
abort() ;
}
// recuperer magB0 (input)
double magB0 = 0. ;
if (par != 0x0)
if (par->get_n_double_mod() > 0) {
magB0 = par->get_double_mod() ;
}
if ( magB0 > Bmax) {
cerr << "Eos_tabul::nbar_ent_p : B > Bmax !" << endl ;
cerr << "B = " << magB0*mag_unit << ", Bmax = " << Bmax*mag_unit << endl ;
abort() ;
}
double p_int, dpdb_int, dpdh_int ;
interpol_herm_2d(*Bfield, *logh, *logp, *dlpsdB, *dlpsdlh, *d2lp, magB0, ent,
p_int, dpdb_int, dpdh_int) ;
double nbar_int = dpdh_int * exp(-ent) ;
return nbar_int ;
}
else{
return 0 ;
}
}
// Energy density from enthalpy
//------------------------------
double Eos_mag::ener_ent_p(double ent, const Param* par ) const {
using namespace Unites_mag ;
if ((ent > hmin - 1.e-12) && (ent < hmin))
ent = hmin ;
if ( ent >= hmin ) {
if (ent > hmax) {
cerr << "Eos_tabul::ener_ent_p : ent > hmax !" << endl ;
abort() ;
}
// recuperer magB0 (input)
double magB0 = 0. ;
if (par != 0x0)
if (par->get_n_double_mod() > 0) {
magB0 = par->get_double_mod() ;
}
if ( magB0 > Bmax) {
cerr << "Eos_tabul::ener_ent_p : B > Bmax !" << endl ;
cerr << "B = " << magB0*mag_unit << ", Bmax = " << Bmax*mag_unit << endl ;
abort() ;
}
double p_int, dpdb_int, dpdh_int ;
interpol_herm_2d(*Bfield, *logh, *logp, *dlpsdB, *dlpsdlh, *d2lp, magB0, ent,
p_int, dpdb_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
//------------------------
double Eos_mag::press_ent_p(double ent, const Param* par ) const {
using namespace Unites_mag ;
if ((ent > hmin - 1.e-12) && (ent < hmin))
ent = hmin ;
if ( ent >= hmin ) {
if (ent > hmax) {
cout << "Eos_mag::press_ent_p : ent > hmax !" << endl ;
abort() ;
}
// recuperer magB0 (input)
double magB0 = 0. ;
if (par != 0x0)
if (par->get_n_double_mod() > 0) {
magB0 = par->get_double_mod() ;
}
if ( magB0 > Bmax) {
cerr << "Eos_tabul::press_ent_p : B > Bmax !" << endl ;
cerr << "B = " << magB0*mag_unit << ", Bmax = " << Bmax*mag_unit << endl ;
abort() ;
}
double p_int, dpdb_int, dpdh_int ;
interpol_herm_2d(*Bfield, *logh, *logp, *dlpsdB, *dlpsdlh, *d2lp, magB0, ent,
p_int, dpdb_int, dpdh_int) ;
return p_int;
}
else{
return 0 ;
}
}
// Magnetization from enthalpy
//----------------------------
double Eos_mag::mag_ent_p(double ent, const Param* par) const {
using namespace Unites_mag ;
if ((ent > hmin - 1.e-12) && (ent < hmin))
ent = hmin ;
if ( ent >= hmin ) {
if (ent > hmax) {
cout << "Eos_mag::mag_ent_p : ent > hmax !" << endl ;
abort() ;
}
// recuperer magB0 (input)
double magB0 = 0. ;
if (par != 0x0)
if (par->get_n_double_mod() > 0) {
magB0 = par->get_double_mod() ;
}
if ( magB0 > Bmax) {
cerr << "Eos_tabul::mag_ent_p : B > Bmax !" << endl ;
cerr << "B = " << magB0*mag_unit << ", Bmax = " << Bmax*mag_unit << endl ;
abort() ;
}
double p_int, dpdb_int, dpdh_int ;
interpol_herm_2d(*Bfield, *logh, *logp, *dlpsdB, *dlpsdlh, *d2lp, magB0, ent,
p_int, dpdb_int, dpdh_int) ;
double magnetization = dpdb_int ;
if (magB0 == 0.)
return 0 ;
else
return mu0*magnetization / magB0 ;
}
else
return 0. ;
}
// dln(n)/ln(H) from enthalpy
//---------------------------
double Eos_mag::der_nbar_ent_p(double ent, const Param* ) const {
c_est_pas_fait("Eos_mag::der_nbar_ent_p" ) ;
return ent ;
}
// dln(e)/ln(H) from enthalpy
//---------------------------
double Eos_mag::der_ener_ent_p(double ent, const Param* ) const {
c_est_pas_fait("Eos_mag::der_ener_enr_p" ) ;
return ent ;
}
// dln(p)/ln(H) from enthalpy
//---------------------------
double Eos_mag::der_press_ent_p(double ent, const Param* ) const {
c_est_pas_fait("Eos_mag::" ) ;
return ent ;
}
// dln(p)/dln(n) from enthalpy
//---------------------------
double Eos_mag::der_press_nbar_p(double ent, const Param*) const {
c_est_pas_fait("Eos_mag::der_press_nbar_p" ) ;
return ent ;
}
}
Lorene::Eos_mag::dlpsdB
Tbl * dlpsdB
Table of .
Definition eos_mag.h:109
Lorene::Eos_mag::mag_ent_p
double mag_ent_p(double ent, const Param *par=0x0) const
Computes the magnetisation.
Definition eos_mag.C:460
Lorene::Eos_mag::d2lp
Tbl * d2lp
Table of .
Definition eos_mag.h:112
Lorene::Eos_mag::operator!=
virtual bool operator!=(const Eos &) const
Comparison operator (difference).
Definition eos_mag.C:196
Lorene::Eos_mag::identify
virtual int identify() const
Returns a number to identify the sub-classe of Eos the object belongs to.
Definition eos_from_file.C:157
Lorene::Eos_mag::logp
Tbl * logp
Table of .
Definition eos_mag.h:97
Lorene::Eos_mag::logh
Tbl * logh
Table of .
Definition eos_mag.h:100
Lorene::Eos_mag::der_press_ent_p
virtual double der_press_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
Definition eos_mag.C:533
Lorene::Eos_mag::press_ent_p
virtual double press_ent_p(double ent, const Param *par=0x0) const
Computes the pressure from the log-enthalpy.
Definition eos_mag.C:419
Lorene::Eos_mag::sauve
virtual void sauve(FILE *) const
Save in a file.
Definition eos_mag.C:171
Lorene::Eos_mag::Bmax
double Bmax
Upper boundary of the magnetic field interval.
Definition eos_mag.h:94
Lorene::Eos_mag::~Eos_mag
virtual ~Eos_mag()
Destructor.
Definition eos_mag.C:158
Lorene::Eos_mag::nbar_ent_p
virtual double nbar_ent_p(double ent, const Param *par=0x0) const
Computes the baryon density from the log-enthalpy.
Definition eos_mag.C:330
Lorene::Eos_mag::Bfield
Tbl * Bfield
Table of .
Definition eos_mag.h:103
Lorene::Eos_mag::der_nbar_ent_p
virtual double der_nbar_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
Definition eos_mag.C:508
Lorene::Eos_mag::operator==
virtual bool operator==(const Eos &) const
Comparison operator (egality).
Definition eos_mag.C:183
Lorene::Eos_mag::hmin
double hmin
Lower boundary of the log-enthalpy interval.
Definition eos_mag.h:88
Lorene::Eos_mag::der_ener_ent_p
virtual double der_ener_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
Definition eos_mag.C:520
Lorene::Eos_mag::ener_ent_p
virtual double ener_ent_p(double ent, const Param *par=0x0) const
Computes the total energy density from the log-enthalpy.
Definition eos_mag.C:373
Lorene::Eos_mag::tablename
string tablename
Name of the file containing the tabulated data.
Definition eos_mag.h:85
Lorene::Eos_mag::read_table
void read_table()
Reads the file containing the table and initializes in the arrays logh , logp and dlpsdlh .
Definition eos_mag.C:224
Lorene::Eos_mag::dlpsdlh
Tbl * dlpsdlh
Table of .
Definition eos_mag.h:106
Lorene::Eos_mag::der_press_nbar_p
virtual double der_press_nbar_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
Definition eos_mag.C:545
Lorene::Eos_mag::Eos_mag
Eos_mag(const char *name_i, const char *table, const char *path)
Standard constructor.
Definition eos_mag.C:98
Lorene::Eos_mag::hmax
double hmax
Upper boundary of the log-enthalpy interval.
Definition eos_mag.h:91
Lorene::Eos_mag::operator>>
virtual ostream & operator>>(ostream &) const
Operator >>.
Definition eos_mag.C:207
Lorene::Eos::identify
virtual int identify() const =0
Returns a number to identify the sub-classe of Eos the object belongs to.
Lorene::Eos::sauve
virtual void sauve(FILE *) const
Save in a file.
Definition eos.C:179
Lorene::Eos::Eos
Eos()
Standard constructor.
Definition eos.C:108
Lorene::Param
Parameter storage.
Definition param.h:125
Lorene::Param::get_n_double_mod
int get_n_double_mod() const
Returns the number of stored modifiable double 's addresses.
Definition param.C:446
Lorene::Param::get_double_mod
double & get_double_mod(int position=0) const
Returns the reference of a stored modifiable double .
Definition param.C:498
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::c_est_pas_fait
void c_est_pas_fait(const char *)
Helpful function to say something is not implemented yet.
Definition c_est_pas_fait.C:73
Lorene
Lorene prototypes.
Definition app_hor.h:64
Unites_mag
Standard electro-magnetic units.