28 char eos_fitting_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Eos/eos_fitting.C,v 1.5 2014/10/13 08:52:53 j_novak Exp $" ;
60 #include "eos_fitting.h" 62 #include "utilitaires.h" 79 const char* path) :
Eos(name_i) {
93 fread(
dataname,
sizeof(
char), 160, fich) ;
105 fich.getline(path, 160) ;
131 fwrite(
dataname,
sizeof(
char), 160, fich) ;
144 for (
int i=0; i<3; i++) {
145 fich.getline(blabla, 120) ;
149 fich >> nb_coef ; fich.getline(blabla, 120) ;
151 for (
int i=0; i<3; i++) {
152 fich.getline(blabla, 120) ;
155 pp =
new double[nb_coef] ;
157 for (
int i=0; i<nb_coef; i++) {
158 fich >>
pp[i] ; fich.getline(blabla, 120) ;
177 if ( ent >
double(0) ) {
186 double trans_dens = msol_si /
pow(g_si*msol_si/c_si/c_si,3.)
189 while (xx > 1.e-15) {
195 while (ent_value > 1.e-15) {
205 double eee = -
pp[7]*yy+
pp[9] ;
206 double fff = -pp[8]*yy+pp[10] ;
207 double ggg = pp[4]*pp[5]*pp[6]*
pow(yy,pp[5]) ;
209 ent_value =
exp(ent) - 1.0
210 -( aaa*bbb - 1. ) * fc(eee)
211 -pp[11]*
pow(yy,pp[12])*fc(-eee)*fc(fff)
212 -pp[13]*
pow(yy,pp[14])*fc(-fff)
213 -( ccc*bbb + aaa*ddd*ggg )*fc(eee)
214 -yy*( aaa*bbb - 1. )*pp[7]*gc(eee)
215 -pp[11]*pp[12]*
pow(yy,pp[12])*fc(-eee)*fc(fff)
216 +pp[11]*
pow(yy,pp[12]+1.)*(pp[7]*gc(-eee)*fc(fff)
217 -pp[8]*fc(-eee)*gc(fff))
218 -pp[13]*
pow(yy,pp[14])*(pp[14]*fc(-fff)
219 -pp[8]*yy*gc(-fff)) ;
227 nb = rhob * trans_dens ;
244 if ( ent >
double(0) ) {
252 double trans_dens = msol_si /
pow(g_si*msol_si/c_si/c_si,3.)
256 double yy = nb / trans_dens ;
260 double eee = -
pp[7]*yy+
pp[9] ;
261 double fff = -pp[8]*yy+pp[10] ;
263 double epsil = ( aaa*bbb - 1. ) * fc(eee)
264 +pp[11]*
pow(yy,pp[12])*fc(-eee)*fc(fff)
265 +pp[13]*
pow(yy,pp[14])*fc(-fff) ;
271 double ee = nb * (1. + epsil) ;
288 if ( ent >
double(0) ) {
296 double trans_dens = msol_si /
pow(g_si*msol_si/c_si/c_si,3.)
300 double yy = nb / trans_dens ;
306 double eee = -
pp[7]*yy+
pp[9] ;
307 double fff = -pp[8]*yy+pp[10] ;
308 double ggg = pp[4]*pp[5]*pp[6]*
pow(yy,pp[5]) ;
310 double ppp = yy*( ccc*bbb + aaa*ddd*ggg )*fc(eee)
311 +yy*yy*( aaa*bbb - 1. )*pp[7]*gc(eee)
312 +pp[11]*pp[12]*
pow(yy,pp[12]+1.)*fc(-eee)*fc(fff)
313 -pp[11]*
pow(yy,pp[12]+2.)*(pp[7]*gc(-eee)*fc(fff)
314 -pp[8]*fc(-eee)*gc(fff))
315 +pp[13]*
pow(yy,pp[14]+1.)*(pp[14]*fc(-fff)-pp[8]*yy*gc(-fff)) ;
321 double pres = ppp * trans_dens ;
338 if ( ent >
double(0) ) {
346 double trans_dens = msol_si /
pow(g_si*msol_si/c_si/c_si,3.)
350 double yy = nb / trans_dens ;
356 double eee = -
pp[7]*yy+
pp[9] ;
357 double fff = -pp[8]*yy+pp[10] ;
358 double ggg = pp[4]*pp[5]*pp[6]*
pow(yy,pp[5]) ;
359 double jjj = pp[0]*pp[1]*pp[1]*
pow(yy,pp[1])
360 +pp[2]*pp[3]*pp[3]*
pow(yy,pp[3]) ;
362 double dlnsdlh =
exp(ent) * ent /
363 ( ( ccc*bbb + aaa*ddd*ggg )*( fc(eee) + 2.*yy*pp[7]*gc(eee) )
364 +yy*pp[7]*( aaa*bbb - 1. )*( 2.*gc(eee) - yy*pp[7]*hc(eee) )
365 +pp[11]*pp[12]*(pp[12]+1.)*
pow(yy,pp[12])*fc(-eee)*fc(fff)
366 +2.*pp[11]*(pp[12]+1.)*
pow(yy,pp[12]+1.)
367 *( -pp[7]*gc(-eee)*fc(fff) + pp[8]*fc(-eee)*gc(fff) )
368 -pp[11]*
pow(yy,pp[12]+2.)*( pp[7]*pp[7]*hc(-eee)*fc(fff)
369 +2.*pp[7]*pp[8]*gc(-eee)*gc(fff)
370 +pp[8]*pp[8]*fc(-eee)*hc(fff) )
371 +pp[13]*pp[14]*(pp[14]+1.)*
pow(yy,pp[14])*fc(-fff)
372 -2.*pp[8]*pp[13]*(pp[14]+1.)*
pow(yy,pp[14]+1.)*gc(-fff)
373 -pp[8]*pp[8]*pp[13]*
pow(yy,pp[14]+2.)*hc(-fff)
374 +( jjj*bbb + 2.*ccc*ddd*ggg
375 + aaa*
pow(1.+pp[4]*
pow(yy,pp[5]),pp[6]-2.)
376 *ggg*(ggg+pp[5]) )*fc(eee)
384 return double(1) /
pp[12] ;
397 if ( ent >
double(0) ) {
405 double trans_dens = msol_si /
pow(g_si*msol_si/c_si/c_si,3.)
409 double yy = nb / trans_dens ;
415 double eee = -
pp[7]*yy+
pp[9] ;
416 double fff = -pp[8]*yy+pp[10] ;
417 double ggg = pp[4]*pp[5]*pp[6]*
pow(yy,pp[5]) ;
421 double dlesdlh = dlnsdlh
422 * (1. + ( ( ccc*bbb + aaa*ddd*ggg )*fc(eee)
423 +yy*pp[7]*( aaa*bbb - 1. )*gc(eee)
424 +pp[11]*pp[12]*
pow(yy,pp[12])*fc(-eee)*fc(fff)
425 +pp[11]*
pow(yy,pp[12]+1.)*( -pp[7]*gc(-eee)*fc(fff)
426 +pp[8]*fc(-eee)*gc(fff) )
427 +pp[13]*pp[14]*
pow(yy,pp[14])*fc(-fff)
428 -pp[8]*pp[13]*
pow(yy,pp[14]+1.)*gc(-fff) )
429 / ( 1. + ( aaa*bbb - 1. )*fc(eee)
430 + pp[11]*
pow(yy,pp[12])*fc(-eee)*fc(fff)
431 + pp[13]*
pow(yy,pp[14])*fc(-fff) ) ) ;
438 return double(1) /
pp[12] ;
451 if ( ent >
double(0) ) {
459 double trans_dens = msol_si /
pow(g_si*msol_si/c_si/c_si,3.)
463 double yy = nb / trans_dens ;
469 double eee = -
pp[7]*yy+
pp[9] ;
470 double fff = -pp[8]*yy+pp[10] ;
471 double ggg = pp[4]*pp[5]*pp[6]*
pow(yy,pp[5]) ;
472 double jjj = pp[0]*pp[1]*pp[1]*
pow(yy,pp[1])
473 +pp[2]*pp[3]*pp[3]*
pow(yy,pp[3]) ;
477 double dlpsdlh = dlnsdlh
478 * ( ( ccc*bbb + aaa*ddd*ggg )*fc(eee)
479 +( jjj*bbb + 2.*ccc*ddd*ggg
480 + aaa*
pow(1.+pp[4]*
pow(yy,pp[5]),pp[6]-2.)*ggg*(ggg+pp[5])
482 +2.*yy*pp[7]*( ccc*bbb + aaa*ddd*ggg )*gc(eee)
483 +yy*pp[7]*( aaa*bbb - 1. )*(2.*gc(eee) - yy*pp[7]*hc(eee))
484 +pp[11]*pp[12]*(pp[12]+1.)*
pow(yy,pp[12])*fc(-eee)*fc(fff)
485 +2.*pp[11]*(pp[12]+1.)*
pow(yy,pp[12]+1.)
486 *( -pp[7]*gc(-eee)*fc(fff) + pp[8]*fc(-eee)*gc(fff) )
487 -pp[11]*
pow(yy,pp[12]+2.)*( pp[7]*pp[7]*hc(-eee)*fc(fff)
488 +2.*pp[7]*pp[8]*gc(-eee)*gc(fff)
489 +pp[8]*pp[8]*fc(-eee)*hc(fff) )
490 +pp[13]*(pp[14]+1.)*
pow(yy,pp[14])*( pp[14]*fc(-fff)
491 -2.*pp[8]*yy*gc(-fff) )
492 -pp[8]*pp[8]*pp[13]*
pow(yy,pp[14]+2.)*hc(-fff) )
493 / ( ( ccc*bbb + aaa*ddd*ggg )*fc(eee)
494 +yy*pp[7]*( aaa*bbb - 1. )*gc(eee)
495 +pp[11]*
pow(yy,pp[12])*( pp[12]*fc(-eee)*fc(fff)
496 -yy*pp[7]*gc(-eee)*fc(fff)
497 +yy*pp[8]*fc(-eee)*gc(fff) )
498 +pp[13]*
pow(yy,pp[14])*( pp[14]*fc(-fff)
499 -yy*pp[8]*gc(-fff) ) ) ;
506 return (
pp[12] + 1.) /
pp[12] ;
516 double fc(
double xx) {
518 double resu = double(1) / (double(1) +
exp(xx)) ;
524 double gc(
double xx) {
529 resu =
exp(-xx) /
pow(
exp(-xx)+
double(1),
double(2)) ;
532 resu =
exp(xx) /
pow(
double(1)+
exp(xx),
double(2)) ;
539 double hc(
double xx) {
544 resu =
exp(-xx) * (
exp(-xx)-double(1)) /
545 pow(
exp(-xx)+double(1),double(3)) ;
548 resu =
exp(xx) * (double(1)-
exp(xx)) /
549 pow(
double(1)+
exp(xx),double(3)) ;
Cmp exp(const Cmp &)
Exponential.
virtual double ener_ent_p(double ent, const Param *par=0x0) const
Computes the total energy density from the log-enthalpy.
Eos_fitting(const char *name_i, const char *data, const char *path)
Standard constructor.
Standard units of space, time and mass.
Equation of state base class.
virtual void sauve(FILE *) const
Save in a file.
virtual void sauve(FILE *) const
Save in a file.
virtual double press_ent_p(double ent, const Param *par=0x0) const
Computes the pressure from the log-enthalpy.
double * pp
Array of the coefficients of the fitting data.
virtual double der_ener_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
virtual double nbar_ent_p(double ent, const Param *par=0x0) const
Computes the baryon density from the log-enthalpy.
void read_coef()
Reading coefficients of the fitting equation for the energy density, the pressure, and the enthalpy.
virtual double der_nbar_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
Cmp pow(const Cmp &, int)
Power .
char dataname[160]
Name of the file containing the fitting data.
virtual ~Eos_fitting()
Destructor.
virtual double der_press_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.