30char eos_mag_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Eos/eos_mag.C,v 1.13 2014/10/13 08:52:53 j_novak Exp $" ;
84#include "utilitaires.h"
89void interpol_herm_2d(
const Tbl&,
const Tbl&,
const Tbl&,
const Tbl&,
const Tbl&,
const Tbl&,
double,
double,
double&,
double&,
double&) ;
99 const char* path) :
Eos(name_i), tablename(path) {
111 :
Eos(name_i), tablename(file_name) {
124 fread(tmp_name,
sizeof(
char), 160, fich) ;
188 cerr <<
"The second EOS is not of type Eos_mag !" << endl ;
210 "EOS of class Eos_mag : tabulated EOS depending on two variables: enthalpy and magnetic field."
212 ost <<
"Table read from " <<
tablename << endl ;
231 cerr <<
"Eos_mag::read_table(): " << endl ;
232 cerr <<
"Problem in opening the EOS file!" << endl ;
233 cerr <<
"Aborting..." << endl ;
237 for (
int i=0; i<5; i++) {
238 fich.ignore(1000,
'\n') ;
242 fich >> nbp1 >> nbp2 ; fich.ignore(1000,
'\n') ;
243 if ((nbp1<=0) || (nbp2<=0)) {
244 cerr <<
"Eos_mag::read_table(): " << endl ;
245 cerr <<
"Wrong value for the number of lines!" << endl ;
246 cerr <<
"nbp1 = " << nbp1 <<
", nbp2 = " << nbp2 << endl ;
247 cerr <<
"Aborting..." << endl ;
251 for (
int i=0; i<3; i++) {
252 fich.ignore(1000,
'\n') ;
269 double c2 = c_si * c_si ;
270 double B_unit = mag_unit / 1.e11 ;
271 double M_unit = B_unit*mu0/(4*M_PI) ;
274 double nb_fm3, rho_cgs, p_cgs, mu_MeV, magB_PG, magM_PG, chi_PGpMeV ;
278 for (
int j=0; j<nbp2; j++) {
280 for (
int i=0; i<nbp1; i++) {
281 fich >> no1 >> nb_fm3 >> rho_cgs >> p_cgs >> mu_MeV
282 >> magB_PG >> magM_PG >> chi_PGpMeV ;
283 fich.ignore(1000,
'\n') ;
284 if ( (nb_fm3<0) || (rho_cgs<0)) {
285 cerr <<
"Eos_mag::read_table(): " << endl ;
286 cerr <<
"Negative value in table!" << endl ;
287 cerr <<
"nb = " << nb_fm3 <<
", rho = " << rho_cgs <<
288 ", p = " << p_cgs << endl ;
289 cerr <<
"Aborting..." << endl ;
293 double psc2 = 0.1*p_cgs/c2 ;
294 double rho_si = rho_cgs*1000. ;
296 double h_read =
log(mu_MeV) ;
297 if ( (i==0) && (j==0) ) ww = h_read ;
298 double h_new = h_read - ww ;
300 logp->
set(j, i) = psc2/rhonuc_si ;
303 dlpsdlh->
set(j, i) = (rho_si + psc2)/rhonuc_si ;
305 d2lp->
set(j, i) = mu_MeV*chi_PGpMeV / (M_unit) ;
310 hmin = (*logh)(0, 0) ;
311 hmax = (*logh)(0, nbp1-1) ;
313 Bmax = (*Bfield)(nbp2-1, 0) ;
334 if ((ent >
hmin - 1.e-12) && (ent <
hmin))
339 cerr <<
"Eos_tabul::nbar_ent_p : ent > hmax !" << endl ;
350 cerr <<
"Eos_tabul::nbar_ent_p : B > Bmax !" << endl ;
351 cerr <<
"B = " << magB0*mag_unit <<
", Bmax = " <<
Bmax*mag_unit << endl ;
355 double p_int, dpdb_int, dpdh_int ;
358 p_int, dpdb_int, dpdh_int) ;
360 double nbar_int = dpdh_int *
exp(-ent) ;
377 if ((ent >
hmin - 1.e-12) && (ent <
hmin))
382 cerr <<
"Eos_tabul::ener_ent_p : ent > hmax !" << endl ;
394 cerr <<
"Eos_tabul::ener_ent_p : B > Bmax !" << endl ;
395 cerr <<
"B = " << magB0*mag_unit <<
", Bmax = " <<
Bmax*mag_unit << endl ;
400 double p_int, dpdb_int, dpdh_int ;
403 p_int, dpdb_int, dpdh_int) ;
405 double nbar_int = dpdh_int *
exp(-ent) ;
407 double f_int = - p_int +
exp(ent) * nbar_int;
423 if ((ent >
hmin - 1.e-12) && (ent <
hmin))
428 cout <<
"Eos_mag::press_ent_p : ent > hmax !" << endl ;
440 cerr <<
"Eos_tabul::press_ent_p : B > Bmax !" << endl ;
441 cerr <<
"B = " << magB0*mag_unit <<
", Bmax = " <<
Bmax*mag_unit << endl ;
445 double p_int, dpdb_int, dpdh_int ;
448 p_int, dpdb_int, dpdh_int) ;
464 if ((ent >
hmin - 1.e-12) && (ent <
hmin))
469 cout <<
"Eos_mag::mag_ent_p : ent > hmax !" << endl ;
481 cerr <<
"Eos_tabul::mag_ent_p : B > Bmax !" << endl ;
482 cerr <<
"B = " << magB0*mag_unit <<
", Bmax = " <<
Bmax*mag_unit << endl ;
486 double p_int, dpdb_int, dpdh_int ;
489 p_int, dpdb_int, dpdh_int) ;
491 double magnetization = dpdb_int ;
496 return mu0*magnetization / magB0 ;
double mag_ent_p(double ent, const Param *par=0x0) const
Computes the magnetisation.
virtual bool operator!=(const Eos &) const
Comparison operator (difference)
virtual int identify() const
Returns a number to identify the sub-classe of Eos the object belongs to.
virtual double der_press_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
virtual double press_ent_p(double ent, const Param *par=0x0) const
Computes the pressure from the log-enthalpy.
virtual void sauve(FILE *) const
Save in a file.
double Bmax
Upper boundary of the magnetic field interval.
virtual ~Eos_mag()
Destructor.
virtual double nbar_ent_p(double ent, const Param *par=0x0) const
Computes the baryon density from the log-enthalpy.
virtual double der_nbar_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
virtual bool operator==(const Eos &) const
Comparison operator (egality)
double hmin
Lower boundary of the log-enthalpy interval.
virtual double der_ener_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
virtual double ener_ent_p(double ent, const Param *par=0x0) const
Computes the total energy density from the log-enthalpy.
string tablename
Name of the file containing the tabulated data.
void read_table()
Reads the file containing the table and initializes in the arrays logh , logp and dlpsdlh .
virtual double der_press_nbar_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
Eos_mag(const char *name_i, const char *table, const char *path)
Standard constructor.
double hmax
Upper boundary of the log-enthalpy interval.
virtual ostream & operator>>(ostream &) const
Operator >>
Equation of state base class.
virtual int identify() const =0
Returns a number to identify the sub-classe of Eos the object belongs to.
virtual void sauve(FILE *) const
Save in a file.
int get_n_double_mod() const
Returns the number of stored modifiable double 's addresses.
double & get_double_mod(int position=0) const
Returns the reference of a stored modifiable double .
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
double & set(int i)
Read/write of a particular element (index i) (1D case)
Cmp exp(const Cmp &)
Exponential.
Cmp log(const Cmp &)
Neperian logarithm.
void c_est_pas_fait(const char *)
Helpful function to say something is not implemented yet.
Standard electro-magnetic units.