80 #include "eos_multi_poly.h" 82 #include "utilitaires.h" 86 double logp(
double,
double,
double,
double,
double,
double) ;
87 double dlpsdlh(
double,
double,
double,
double,
double,
double) ;
88 double dlpsdlnb(
double,
double,
double,
double,
double) ;
98 double logP1_i,
double* logRho_i,
100 :
Eos(
"Multi-polytropic EOS"),
101 npeos(npoly), kappa0(kappa0_i), logP1(logP1_i), m0(double(1)) {
107 for (
int l=0; l<
npeos; l++) {
108 gamma[l] = gamma_i[l] ;
113 for (
int l=0; l<
npeos-1; l++) {
119 for (
int l=0; l<
npeos-1; l++) {
131 npeos(eosmp.npeos), kappa0(eosmp.kappa0), logP1(eosmp.logP1),
136 for (
int l=0; l<
npeos; l++) {
142 for (
int l=0; l<
npeos-1; l++) {
148 for (
int l=0; l<
npeos; l++) {
154 for (
int l=0; l<
npeos-1; l++) {
160 for (
int l=0; l<
npeos-1; l++) {
166 for (
int l=0; l<
npeos-1; l++) {
172 for (
int l=0; l<
npeos; l++) {
185 for (
int l=0; l<
npeos; l++) {
194 for (
int l=0; l<
npeos-1; l++) {
200 for (
int l=0; l<
npeos-1; l++) {
215 fich >>
npeos ; fich.getline(blabla, 80) ;
219 for (
int l=0; l<
npeos; l++) {
220 fich >>
gamma[l] ; fich.getline(blabla, 80) ;
223 fich >>
kappa0 ; fich.getline(blabla, 80) ;
224 fich >>
logP1 ; fich.getline(blabla, 80) ;
228 for (
int l=0; l<
npeos-1; l++) {
229 fich >>
logRho[l] ; fich.getline(blabla, 80) ;
234 for (
int l=0; l<
npeos-1; l++) {
235 fich >>
decInc[l] ; fich.getline(blabla, 80) ;
263 cout <<
"Eos_multi_poly::operator= : not implemented yet !" << endl ;
277 double* kappa_cgs =
new double [
npeos] ;
285 kappa_cgs[2] = kappa_cgs[1]
290 for (
int l=3; l<
npeos; l++) {
292 kappa_cgs[l] = kappa_cgs[l-1]
303 double rhonuc_cgs = rhonuc_si * 1.e-3 ;
305 for (
int l=0; l<
npeos; l++) {
306 kappa[l] = kappa_cgs[l] *
pow( rhonuc_cgs,
gamma[l] -
double(1) ) ;
310 delete [] kappa_cgs ;
319 for (
int l=0; l<
npeos-1; l++) {
326 / (
gamma[l]-double(1))
328 / (
gamma[l+1]-
double(1)) ) ;
333 / (
gamma[l]-
double(1)) /
m0 ) ;
349 cout <<
"The second EOS is not of type Eos_multi_poly !" << endl ;
358 cout <<
"The two Eos_multi_poly have " 359 <<
"different number of polytropes : " 364 for (
int l=0; l<
npeos; l++) {
366 cout <<
"The two Eos_multi_poly have different gamma " 367 <<
gamma[l] <<
" <-> " 373 for (
int l=0; l<
npeos; l++) {
375 cout <<
"The two Eos_multi_poly have different kappa " 376 <<
kappa[l] <<
" <-> " 404 for (
int l=0; l<
npeos; l++) {
411 for (
int l=0; l<
npeos-1; l++) {
415 for (
int l=0; l<
npeos-1; l++) {
426 ost <<
"EOS of class Eos_multi_poly " 427 <<
"(multiple polytropic equation of state) : " << endl ;
429 ost <<
" Number of polytropes : " 430 <<
npeos << endl << endl ;
432 double rhonuc_cgs = rhonuc_si * 1.e-3 ;
435 for (
int l=0; l<
npeos; l++) {
436 ost <<
" EOS in region " << l <<
" : " << endl ;
437 ost <<
" ---------------" << endl ;
438 ost <<
" gamma : " <<
gamma[l] << endl ;
439 ost <<
" kappa : " <<
kappa[l]
440 <<
" [Lorene units: rho_nuc c^2 / n_nuc^gamma]" << endl ;
442 double kappa_cgs =
kappa[l]
443 *
pow( rhonuc_cgs,
double(1) -
gamma[l] ) ;
445 ost <<
" : " << kappa_cgs
446 <<
" [(g/cm^3)^{1-gamma}]" << endl ;
450 ost <<
" Exponent of the pressure at the fiducial density rho_1" 452 ost <<
" ------------------------------------------------------" 454 ost <<
" log P1 : " <<
logP1 << endl ;
457 ost <<
" Exponent of fiducial densities" << endl ;
458 ost <<
" ------------------------------" << endl ;
459 for (
int l=0; l<
npeos-1; l++) {
460 ost <<
" log rho[" << l <<
"] : " <<
logRho[l] << endl ;
464 for (
int l=0; l<
npeos-1; l++) {
465 ost <<
" Critical density and enthalpy between domains " 466 << l <<
" and " << l+1 <<
" : " << endl ;
467 ost <<
" -----------------------------------------------------" 469 ost <<
" num. dens. : " <<
nbCrit[l] <<
" [Lorene units: n_nuc]" 471 ost <<
" density : " <<
nbCrit[l] * rhonuc_cgs <<
" [g/cm^3]" 474 ost <<
" ln(ent) : " <<
entCrit[l] << endl ;
478 for (
int l=0; l<
npeos; l++) {
479 ost <<
" Relat. chem. pot. in region " << l <<
" : " << endl ;
480 ost <<
" -----------------------------" << endl ;
481 ost <<
" mu : " <<
mu0[l] <<
" [m_B c^2]" << endl ;
501 for (
int l=0; l<
npeos-1; l++) {
507 double mgam1skapgam = 1. ;
510 if ( ent >
double(0) ) {
514 return pow( mgam1skapgam*(
exp(ent)-
double(1)),
515 double(1)/(
gamma[0]-
double(1)) ) ;
528 if ( ent < entLarge ) {
530 double log10H =
log10( ent ) ;
531 double log10HSmall =
log10( entSmall ) ;
532 double log10HLarge =
log10( entLarge ) ;
533 double dH = log10HLarge - log10HSmall ;
534 double uu = (log10H - log10HSmall) / dH ;
536 double mgam1skapgamSmall =
m0*(
gamma[i-1]-double(1))
538 double mgam1skapgamLarge =
m0*(
gamma[i]-double(1))
541 double nnSmall =
pow( mgam1skapgamSmall
543 double(1)/(
gamma[i-1]-
double(1)) ) ;
544 double nnLarge =
pow( mgam1skapgamLarge
546 double(1)/(
gamma[i]-
double(1)) ) ;
551 double eeSmall =
mu0[i-1] * nnSmall
552 + ppSmall / (
gamma[i-1] - double(1)) ;
553 double eeLarge =
mu0[i] * nnLarge
554 + ppLarge / (
gamma[i] - double(1)) ;
556 double log10PSmall =
log10( ppSmall ) ;
557 double log10PLarge =
log10( ppLarge ) ;
559 double dlpsdlhSmall = entSmall
560 * (double(1) + eeSmall / ppSmall) ;
561 double dlpsdlhLarge = entLarge
562 * (double(1) + eeLarge / ppLarge) ;
564 double log10PInterpolate = logp(log10PSmall, log10PLarge,
565 dlpsdlhSmall, dlpsdlhLarge,
568 double dlpsdlhInterpolate = dlpsdlh(log10PSmall, log10PLarge,
569 dlpsdlhSmall, dlpsdlhLarge,
572 double pp =
pow(
double(10), log10PInterpolate) ;
574 return pp / ent * dlpsdlhInterpolate *
exp(-ent) /
m0 ;
583 double(1)/(
gamma[i]-
double(1)) ) ;
599 for (
int l=0; l<
npeos-1; l++) {
605 double mgam1skapgam = 1. ;
611 if ( ent >
double(0) ) {
615 nn =
pow( mgam1skapgam*(
exp(ent)-
double(1)),
616 double(1)/(
gamma[0]-
double(1)) ) ;
620 return mu0[0] * nn + pp / (
gamma[0] - double(1)) ;
633 if ( ent < entLarge ) {
635 double log10H =
log10( ent ) ;
636 double log10HSmall =
log10( entSmall ) ;
637 double log10HLarge =
log10( entLarge ) ;
638 double dH = log10HLarge - log10HSmall ;
639 double uu = (log10H - log10HSmall) / dH ;
641 double mgam1skapgamSmall =
m0*(
gamma[i-1]-double(1))
643 double mgam1skapgamLarge =
m0*(
gamma[i]-double(1))
646 double nnSmall =
pow( mgam1skapgamSmall
648 double(1)/(
gamma[i-1]-
double(1)) ) ;
649 double nnLarge =
pow( mgam1skapgamLarge
651 double(1)/(
gamma[i]-
double(1)) ) ;
656 double eeSmall =
mu0[i-1] * nnSmall
657 + ppSmall / (
gamma[i-1] - double(1)) ;
658 double eeLarge =
mu0[i] * nnLarge
659 + ppLarge / (
gamma[i] - double(1)) ;
661 double log10PSmall =
log10( ppSmall ) ;
662 double log10PLarge =
log10( ppLarge ) ;
664 double dlpsdlhSmall = entSmall
665 * (double(1) + eeSmall / ppSmall) ;
666 double dlpsdlhLarge = entLarge
667 * (double(1) + eeLarge / ppLarge) ;
669 double log10PInterpolate = logp(log10PSmall, log10PLarge,
670 dlpsdlhSmall, dlpsdlhLarge,
673 double dlpsdlhInterpolate = dlpsdlh(log10PSmall, log10PLarge,
674 dlpsdlhSmall, dlpsdlhLarge,
677 pp =
pow(
double(10), log10PInterpolate) ;
679 return pp / ent * dlpsdlhInterpolate - pp ;
687 double(1)/(
gamma[i]-
double(1)) ) ;
691 return mu0[i] * nn + pp / (
gamma[i] - double(1)) ;
708 for (
int l=0; l<
npeos-1; l++) {
714 double mgam1skapgam = 1. ;
719 if ( ent >
double(0) ) {
723 nn =
pow( mgam1skapgam*(
exp(ent)-
double(1)),
724 double(1)/(
gamma[0]-
double(1)) ) ;
739 if ( ent < entLarge ) {
741 double log10H =
log10( ent ) ;
742 double log10HSmall =
log10( entSmall ) ;
743 double log10HLarge =
log10( entLarge ) ;
744 double dH = log10HLarge - log10HSmall ;
745 double uu = (log10H - log10HSmall) / dH ;
747 double mgam1skapgamSmall =
m0*(
gamma[i-1]-double(1))
749 double mgam1skapgamLarge =
m0*(
gamma[i]-double(1))
752 double nnSmall =
pow( mgam1skapgamSmall
754 double(1)/(
gamma[i-1]-
double(1)) ) ;
755 double nnLarge =
pow( mgam1skapgamLarge
757 double(1)/(
gamma[i]-
double(1)) ) ;
762 double eeSmall =
mu0[i-1] * nnSmall
763 + ppSmall / (
gamma[i-1] - double(1)) ;
764 double eeLarge =
mu0[i] * nnLarge
765 + ppLarge / (
gamma[i] - double(1)) ;
767 double log10PSmall =
log10( ppSmall ) ;
768 double log10PLarge =
log10( ppLarge ) ;
770 double dlpsdlhSmall = entSmall
771 * (double(1) + eeSmall / ppSmall) ;
772 double dlpsdlhLarge = entLarge
773 * (double(1) + eeLarge / ppLarge) ;
775 double log10PInterpolate = logp(log10PSmall, log10PLarge,
776 dlpsdlhSmall, dlpsdlhLarge,
779 return pow(
double(10), log10PInterpolate) ;
787 double(1)/(
gamma[i]-
double(1)) ) ;
806 for (
int l=0; l<
npeos-1; l++) {
814 if ( ent >
double(0) ) {
816 if ( ent < 1.e-13 ) {
818 return (
double(1) + ent/
double(2) + ent*ent/
double(12) )
819 / (
gamma[0] - double(1)) ;
824 return ent / (double(1) -
exp(-ent))
825 / (
gamma[0] -
double(1)) ;
832 return double(1) / (
gamma[0] - double(1)) ;
848 return ent / (double(1) - (
mu0[i]/
m0) *
exp(-ent))
849 / (
gamma[i] -
double(1)) ;
865 for (
int l=0; l<
npeos-1; l++) {
871 double mgam1skapgam = 1. ;
878 if ( ent >
double(0) ) {
882 nn =
pow( mgam1skapgam*(
exp(ent)-
double(1)),
883 double(1)/(
gamma[0]-
double(1)) ) ;
887 ee =
mu0[0] * nn + pp / (
gamma[0] - double(1)) ;
889 if ( ent < 1.e-13 ) {
891 return (
double(1) + ent/double(2) + ent*ent/double(12) )
892 / (
gamma[0] -
double(1)) * (
double(1) + pp / ee) ;
897 return ent / (double(1) -
exp(-ent))
898 / (
gamma[0] -
double(1)) * (
double(1) + pp / ee) ;
906 return double(1) / (
gamma[0] - double(1)) ;
916 if ( ent < entLarge ) {
918 double log10H =
log10( ent ) ;
919 double log10HSmall =
log10( entSmall ) ;
920 double log10HLarge =
log10( entLarge ) ;
921 double dH = log10HLarge - log10HSmall ;
922 double uu = (log10H - log10HSmall) / dH ;
924 double mgam1skapgamSmall =
m0*(
gamma[i-1]-double(1))
926 double mgam1skapgamLarge =
m0*(
gamma[i]-double(1))
929 double nnSmall =
pow( mgam1skapgamSmall
931 double(1)/(
gamma[i-1]-
double(1)) ) ;
932 double nnLarge =
pow( mgam1skapgamLarge
934 double(1)/(
gamma[i]-
double(1)) ) ;
939 double eeSmall =
mu0[i-1] * nnSmall
940 + ppSmall / (
gamma[i-1] - double(1)) ;
941 double eeLarge =
mu0[i] * nnLarge
942 + ppLarge / (
gamma[i] - double(1)) ;
944 double log10PSmall =
log10( ppSmall ) ;
945 double log10PLarge =
log10( ppLarge ) ;
947 double dlpsdlhSmall = entSmall
948 * (double(1) + eeSmall / ppSmall) ;
949 double dlpsdlhLarge = entLarge
950 * (double(1) + eeLarge / ppLarge) ;
952 double log10PInterpolate = logp(log10PSmall, log10PLarge,
953 dlpsdlhSmall, dlpsdlhLarge,
956 double dlpsdlhInterpolate = dlpsdlh(log10PSmall, log10PLarge,
957 dlpsdlhSmall, dlpsdlhLarge,
960 pp =
pow(
double(10), log10PInterpolate) ;
962 ee = pp / ent * dlpsdlhInterpolate - pp ;
975 double(1)/(
gamma[i]-
double(1)) ) ;
979 ee =
mu0[i] * nn + pp / (
gamma[i] - double(1)) ;
981 return ent / (double(1) - (
mu0[i]/
m0) *
exp(-ent))
982 / (
gamma[i] -
double(1)) * (
double(1) + pp / ee) ;
998 for (
int l=0; l<
npeos-1; l++) {
1006 if ( ent >
double(0) ) {
1008 if ( ent < 1.e-13 ) {
1011 * ( double(1) + ent/double(2) + ent*ent/double(12) )
1012 / (
gamma[0] -
double(1)) ;
1017 return gamma[0] * ent / (double(1) -
exp(-ent))
1018 / (
gamma[0] -
double(1)) ;
1035 if ( ent < entLarge ) {
1037 double log10H =
log10( ent ) ;
1038 double log10HSmall =
log10( entSmall ) ;
1039 double log10HLarge =
log10( entLarge ) ;
1040 double dH = log10HLarge - log10HSmall ;
1041 double uu = (log10H - log10HSmall) / dH ;
1043 double mgam1skapgamSmall =
m0*(
gamma[i-1]-double(1))
1045 double mgam1skapgamLarge =
m0*(
gamma[i]-double(1))
1048 double nnSmall =
pow( mgam1skapgamSmall
1050 double(1)/(
gamma[i-1]-
double(1)) ) ;
1051 double nnLarge =
pow( mgam1skapgamLarge
1053 double(1)/(
gamma[i]-
double(1)) ) ;
1058 double eeSmall =
mu0[i-1] * nnSmall
1059 + ppSmall / (
gamma[i-1] - double(1)) ;
1060 double eeLarge =
mu0[i] * nnLarge
1061 + ppLarge / (
gamma[i] - double(1)) ;
1063 double log10PSmall =
log10( ppSmall ) ;
1064 double log10PLarge =
log10( ppLarge ) ;
1066 double dlpsdlhSmall = entSmall
1067 * (double(1) + eeSmall / ppSmall) ;
1068 double dlpsdlhLarge = entLarge
1069 * (double(1) + eeLarge / ppLarge) ;
1071 double dlpsdlhInterpolate = dlpsdlh(log10PSmall, log10PLarge,
1072 dlpsdlhSmall, dlpsdlhLarge,
1074 return dlpsdlhInterpolate ;
1079 return gamma[i] * ent / (double(1) - (
mu0[i]/
m0) *
exp(-ent))
1080 / (
gamma[i] -
double(1)) ;
1096 for (
int l=0; l<
npeos-1; l++) {
1112 if ( ent < entLarge ) {
1114 double log10H =
log10( ent ) ;
1115 double log10HSmall =
log10( entSmall ) ;
1116 double log10HLarge =
log10( entLarge ) ;
1118 double dlpsdlnbInterpolate = dlpsdlnb(log10HSmall, log10HLarge,
1122 return dlpsdlnbInterpolate ;
1139 double logp(
double log10PressSmall,
double log10PressLarge,
1140 double dlpsdlhSmall,
double dlpsdlhLarge,
1141 double dx,
double u) {
1143 double resu = log10PressSmall * (double(2) * u * u * u
1144 - double(3) * u * u + double(1))
1145 + log10PressLarge * (
double(3) * u * u - double(2) * u * u * u)
1146 + dlpsdlhSmall * dx * (u * u * u -
double(2) * u * u + u)
1147 - dlpsdlhLarge * dx * (u * u - u * u * u) ;
1153 double dlpsdlh(
double log10PressSmall,
double log10PressLarge,
1154 double dlpsdlhSmall,
double dlpsdlhLarge,
1155 double dx,
double u) {
1157 double resu = double(6) * (log10PressSmall - log10PressLarge)
1159 + dlpsdlhSmall * (double(3) * u * u - double(4) * u + double(1))
1160 + dlpsdlhLarge * (
double(3) * u * u - double(2) * u) ;
1166 double dlpsdlnb(
double log10HSmall,
double log10HLarge,
1167 double dlpsdlnbSmall,
double dlpsdlnbLarge,
1170 double resu = log10H * (dlpsdlnbSmall - dlpsdlnbLarge)
1171 / (log10HSmall - log10HLarge)
1172 + (log10HSmall * dlpsdlnbLarge - log10HLarge * dlpsdlnbSmall)
1173 / (log10HSmall - log10HLarge) ;
Cmp log(const Cmp &)
Neperian logarithm.
Cmp exp(const Cmp &)
Exponential.
virtual ~Eos_multi_poly()
Destructor.
virtual double der_press_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
virtual bool operator==(const Eos &) const
Read/write kappa.
double * kappa
Array (size: npeos) of pressure coefficient [unit: ], where and .
Standard units of space, time and mass.
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 double ener_ent_p(double ent, const Param *par=0x0) const
Computes the total energy density from the log-enthalpy.
const int & get_npeos() const
Returns the number of polytropes npeos.
virtual double der_press_nbar_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
void operator=(const Eos_multi_poly &)
Assignment to another Eos_multi_poly.
Eos_multi_poly(int npoly, double *gamma_i, double kappa0_i, double logP1_i, double *logRho_i, double *decInc_i)
Standard constructor (sets m0 to 1).
virtual void sauve(FILE *) const
Save in a file.
double kappa0
Pressure coefficient for the crust [unit: ].
double * gamma
Array (size: npeos) of adiabatic index .
virtual double der_ener_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
double * entCrit
Array (size npeos - 1) of the critical enthalpy at which the polytropic EOS changes its index and con...
virtual double nbar_ent_p(double ent, const Param *par=0x0) const
Computes the baryon density from the log-enthalpy.
double logP1
Exponent of the pressure at the fiducial density .
double * nbCrit
Array (size npeos - 1) of the number density at which the polytropic EOS changes its index and consta...
int fwrite_be(const int *aa, int size, int nb, FILE *fich)
Writes integer(s) into a binary file according to the big endian convention.
const double & get_kappa(int n) const
Returns the pressure coefficient [unit: ], where and .
Cmp pow(const Cmp &, int)
Power .
double * mu0
Array (size: npeos) of the relativistic chemical potential at zero pressure [unit: ...
Base class for a multiple polytropic equation of state.
virtual ostream & operator>>(ostream &) const
Operator >>
int fread_be(int *aa, int size, int nb, FILE *fich)
Reads integer(s) from a binary file according to the big endian convention.
double * decInc
Array (size npeos - 1) of the percentage which detemines the terminating enthalpy for lower density a...
virtual bool operator!=(const Eos &) const
Comparison operator (difference)
int npeos
Number of polytropic equations of state.
void set_auxiliary()
Computes the auxiliary quantities.
Cmp log10(const Cmp &)
Basis 10 logarithm.
double m0
Individual particule mass [unit: ].
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 * logRho
Array (size: npeos - 1) of the exponent of fiducial densities.
const double & get_gamma(int n) const
Returns the adiabatic index .
virtual int identify() const
Returns a number to identify the sub-classe of Eos the object belongs to.
virtual double der_nbar_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.