29char eos_multi_poly_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Eos/eos_multi_poly.C,v 1.10 2014/12/09 14:07:14 j_novak Exp $" ;
77#include "eos_multi_poly.h"
79#include "utilitaires.h"
83double logp(
double,
double,
double,
double,
double,
double) ;
84double dlpsdlh(
double,
double,
double,
double,
double,
double) ;
85double dlpsdlnb(
double,
double,
double,
double,
double) ;
95 double logP1_i,
double* logRho_i,
97 :
Eos(
"Multi-polytropic EOS"),
98 npeos(npoly), kappa0(kappa0_i), logP1(logP1_i), m0(double(1)) {
104 for (
int l=0; l<
npeos; l++) {
105 gamma[l] = gamma_i[l] ;
110 for (
int l=0; l<
npeos-1; l++) {
116 for (
int l=0; l<
npeos-1; l++) {
128 npeos(eosmp.npeos), kappa0(eosmp.kappa0), logP1(eosmp.logP1),
133 for (
int l=0; l<
npeos; l++) {
139 for (
int l=0; l<
npeos-1; l++) {
145 for (
int l=0; l<
npeos; l++) {
151 for (
int l=0; l<
npeos-1; l++) {
157 for (
int l=0; l<
npeos-1; l++) {
163 for (
int l=0; l<
npeos-1; l++) {
169 for (
int l=0; l<
npeos; l++) {
182 for (
int l=0; l<
npeos; l++) {
191 for (
int l=0; l<
npeos-1; l++) {
197 for (
int l=0; l<
npeos-1; l++) {
212 fich >>
npeos ; fich.getline(blabla, 80) ;
216 for (
int l=0; l<
npeos; l++) {
217 fich >>
gamma[l] ; fich.getline(blabla, 80) ;
220 fich >>
kappa0 ; fich.getline(blabla, 80) ;
221 fich >>
logP1 ; fich.getline(blabla, 80) ;
225 for (
int l=0; l<
npeos-1; l++) {
226 fich >>
logRho[l] ; fich.getline(blabla, 80) ;
231 for (
int l=0; l<
npeos-1; l++) {
232 fich >>
decInc[l] ; fich.getline(blabla, 80) ;
260 cout <<
"Eos_multi_poly::operator= : not implemented yet !" << endl ;
274 double* kappa_cgs =
new double [
npeos] ;
282 kappa_cgs[2] = kappa_cgs[1]
287 for (
int l=3; l<
npeos; l++) {
289 kappa_cgs[l] = kappa_cgs[l-1]
300 double rhonuc_cgs = rhonuc_si * 1.e-3 ;
302 for (
int l=0; l<
npeos; l++) {
303 kappa[l] = kappa_cgs[l] *
pow( rhonuc_cgs,
gamma[l] -
double(1) ) ;
307 delete [] kappa_cgs ;
316 for (
int l=0; l<
npeos-1; l++) {
323 / (
gamma[l]-double(1))
325 / (
gamma[l+1]-
double(1)) ) ;
330 / (
gamma[l]-
double(1)) /
m0 ) ;
346 cout <<
"The second EOS is not of type Eos_multi_poly !" << endl ;
355 cout <<
"The two Eos_multi_poly have "
356 <<
"different number of polytropes : "
361 for (
int l=0; l<
npeos; l++) {
363 cout <<
"The two Eos_multi_poly have different gamma "
364 <<
gamma[l] <<
" <-> "
370 for (
int l=0; l<
npeos; l++) {
372 cout <<
"The two Eos_multi_poly have different kappa "
373 <<
kappa[l] <<
" <-> "
401 for (
int l=0; l<
npeos; l++) {
408 for (
int l=0; l<
npeos-1; l++) {
412 for (
int l=0; l<
npeos-1; l++) {
423 ost <<
"EOS of class Eos_multi_poly "
424 <<
"(multiple polytropic equation of state) : " << endl ;
426 ost <<
" Number of polytropes : "
427 <<
npeos << endl << endl ;
429 double rhonuc_cgs = rhonuc_si * 1.e-3 ;
432 for (
int l=0; l<
npeos; l++) {
433 ost <<
" EOS in region " << l <<
" : " << endl ;
434 ost <<
" ---------------" << endl ;
435 ost <<
" gamma : " <<
gamma[l] << endl ;
436 ost <<
" kappa : " <<
kappa[l]
437 <<
" [Lorene units: rho_nuc c^2 / n_nuc^gamma]" << endl ;
439 double kappa_cgs =
kappa[l]
440 *
pow( rhonuc_cgs,
double(1) -
gamma[l] ) ;
442 ost <<
" : " << kappa_cgs
443 <<
" [(g/cm^3)^{1-gamma}]" << endl ;
447 ost <<
" Exponent of the pressure at the fiducial density rho_1"
449 ost <<
" ------------------------------------------------------"
451 ost <<
" log P1 : " <<
logP1 << endl ;
454 ost <<
" Exponent of fiducial densities" << endl ;
455 ost <<
" ------------------------------" << endl ;
456 for (
int l=0; l<
npeos-1; l++) {
457 ost <<
" log rho[" << l <<
"] : " <<
logRho[l] << endl ;
461 for (
int l=0; l<
npeos-1; l++) {
462 ost <<
" Critical density and enthalpy between domains "
463 << l <<
" and " << l+1 <<
" : " << endl ;
464 ost <<
" -----------------------------------------------------"
466 ost <<
" num. dens. : " <<
nbCrit[l] <<
" [Lorene units: n_nuc]"
468 ost <<
" density : " <<
nbCrit[l] * rhonuc_cgs <<
" [g/cm^3]"
471 ost <<
" ln(ent) : " <<
entCrit[l] << endl ;
475 for (
int l=0; l<
npeos; l++) {
476 ost <<
" Relat. chem. pot. in region " << l <<
" : " << endl ;
477 ost <<
" -----------------------------" << endl ;
478 ost <<
" mu : " <<
mu0[l] <<
" [m_B c^2]" << endl ;
498 for (
int l=0; l<
npeos-1; l++) {
504 double mgam1skapgam = 1. ;
507 if ( ent >
double(0) ) {
511 return pow( mgam1skapgam*(
exp(ent)-
double(1)),
512 double(1)/(
gamma[0]-
double(1)) ) ;
525 if ( ent < entLarge ) {
527 double log10H =
log10( ent ) ;
528 double log10HSmall =
log10( entSmall ) ;
529 double log10HLarge =
log10( entLarge ) ;
530 double dH = log10HLarge - log10HSmall ;
531 double uu = (log10H - log10HSmall) / dH ;
533 double mgam1skapgamSmall =
m0*(
gamma[i-1]-double(1))
535 double mgam1skapgamLarge =
m0*(
gamma[i]-double(1))
538 double nnSmall =
pow( mgam1skapgamSmall
540 double(1)/(
gamma[i-1]-
double(1)) ) ;
541 double nnLarge =
pow( mgam1skapgamLarge
543 double(1)/(
gamma[i]-
double(1)) ) ;
548 double eeSmall =
mu0[i-1] * nnSmall
549 + ppSmall / (
gamma[i-1] - double(1)) ;
550 double eeLarge =
mu0[i] * nnLarge
551 + ppLarge / (
gamma[i] - double(1)) ;
553 double log10PSmall =
log10( ppSmall ) ;
554 double log10PLarge =
log10( ppLarge ) ;
556 double dlpsdlhSmall = entSmall
557 * (double(1) + eeSmall / ppSmall) ;
558 double dlpsdlhLarge = entLarge
559 * (double(1) + eeLarge / ppLarge) ;
561 double log10PInterpolate = logp(log10PSmall, log10PLarge,
562 dlpsdlhSmall, dlpsdlhLarge,
565 double dlpsdlhInterpolate = dlpsdlh(log10PSmall, log10PLarge,
566 dlpsdlhSmall, dlpsdlhLarge,
569 double pp =
pow(
double(10), log10PInterpolate) ;
571 return pp / ent * dlpsdlhInterpolate *
exp(-ent) /
m0 ;
580 double(1)/(
gamma[i]-
double(1)) ) ;
596 for (
int l=0; l<
npeos-1; l++) {
602 double mgam1skapgam = 1. ;
608 if ( ent >
double(0) ) {
612 nn =
pow( mgam1skapgam*(
exp(ent)-
double(1)),
613 double(1)/(
gamma[0]-
double(1)) ) ;
617 return mu0[0] * nn + pp / (
gamma[0] - double(1)) ;
630 if ( ent < entLarge ) {
632 double log10H =
log10( ent ) ;
633 double log10HSmall =
log10( entSmall ) ;
634 double log10HLarge =
log10( entLarge ) ;
635 double dH = log10HLarge - log10HSmall ;
636 double uu = (log10H - log10HSmall) / dH ;
638 double mgam1skapgamSmall =
m0*(
gamma[i-1]-double(1))
640 double mgam1skapgamLarge =
m0*(
gamma[i]-double(1))
643 double nnSmall =
pow( mgam1skapgamSmall
645 double(1)/(
gamma[i-1]-
double(1)) ) ;
646 double nnLarge =
pow( mgam1skapgamLarge
648 double(1)/(
gamma[i]-
double(1)) ) ;
653 double eeSmall =
mu0[i-1] * nnSmall
654 + ppSmall / (
gamma[i-1] - double(1)) ;
655 double eeLarge =
mu0[i] * nnLarge
656 + ppLarge / (
gamma[i] - double(1)) ;
658 double log10PSmall =
log10( ppSmall ) ;
659 double log10PLarge =
log10( ppLarge ) ;
661 double dlpsdlhSmall = entSmall
662 * (double(1) + eeSmall / ppSmall) ;
663 double dlpsdlhLarge = entLarge
664 * (double(1) + eeLarge / ppLarge) ;
666 double log10PInterpolate = logp(log10PSmall, log10PLarge,
667 dlpsdlhSmall, dlpsdlhLarge,
670 double dlpsdlhInterpolate = dlpsdlh(log10PSmall, log10PLarge,
671 dlpsdlhSmall, dlpsdlhLarge,
674 pp =
pow(
double(10), log10PInterpolate) ;
676 return pp / ent * dlpsdlhInterpolate - pp ;
684 double(1)/(
gamma[i]-
double(1)) ) ;
688 return mu0[i] * nn + pp / (
gamma[i] - double(1)) ;
705 for (
int l=0; l<
npeos-1; l++) {
711 double mgam1skapgam = 1. ;
716 if ( ent >
double(0) ) {
720 nn =
pow( mgam1skapgam*(
exp(ent)-
double(1)),
721 double(1)/(
gamma[0]-
double(1)) ) ;
736 if ( ent < entLarge ) {
738 double log10H =
log10( ent ) ;
739 double log10HSmall =
log10( entSmall ) ;
740 double log10HLarge =
log10( entLarge ) ;
741 double dH = log10HLarge - log10HSmall ;
742 double uu = (log10H - log10HSmall) / dH ;
744 double mgam1skapgamSmall =
m0*(
gamma[i-1]-double(1))
746 double mgam1skapgamLarge =
m0*(
gamma[i]-double(1))
749 double nnSmall =
pow( mgam1skapgamSmall
751 double(1)/(
gamma[i-1]-
double(1)) ) ;
752 double nnLarge =
pow( mgam1skapgamLarge
754 double(1)/(
gamma[i]-
double(1)) ) ;
759 double eeSmall =
mu0[i-1] * nnSmall
760 + ppSmall / (
gamma[i-1] - double(1)) ;
761 double eeLarge =
mu0[i] * nnLarge
762 + ppLarge / (
gamma[i] - double(1)) ;
764 double log10PSmall =
log10( ppSmall ) ;
765 double log10PLarge =
log10( ppLarge ) ;
767 double dlpsdlhSmall = entSmall
768 * (double(1) + eeSmall / ppSmall) ;
769 double dlpsdlhLarge = entLarge
770 * (double(1) + eeLarge / ppLarge) ;
772 double log10PInterpolate = logp(log10PSmall, log10PLarge,
773 dlpsdlhSmall, dlpsdlhLarge,
776 return pow(
double(10), log10PInterpolate) ;
784 double(1)/(
gamma[i]-
double(1)) ) ;
803 for (
int l=0; l<
npeos-1; l++) {
811 if ( ent >
double(0) ) {
813 if ( ent < 1.e-13 ) {
815 return (
double(1) + ent/
double(2) + ent*ent/
double(12) )
816 / (
gamma[0] - double(1)) ;
821 return ent / (double(1) -
exp(-ent))
822 / (
gamma[0] -
double(1)) ;
829 return double(1) / (
gamma[0] - double(1)) ;
845 return ent / (double(1) - (
mu0[i]/
m0) *
exp(-ent))
846 / (
gamma[i] -
double(1)) ;
862 for (
int l=0; l<
npeos-1; l++) {
868 double mgam1skapgam = 1. ;
875 if ( ent >
double(0) ) {
879 nn =
pow( mgam1skapgam*(
exp(ent)-
double(1)),
880 double(1)/(
gamma[0]-
double(1)) ) ;
884 ee =
mu0[0] * nn + pp / (
gamma[0] - double(1)) ;
886 if ( ent < 1.e-13 ) {
888 return (
double(1) + ent/double(2) + ent*ent/double(12) )
889 / (
gamma[0] -
double(1)) * (
double(1) + pp / ee) ;
894 return ent / (double(1) -
exp(-ent))
895 / (
gamma[0] -
double(1)) * (
double(1) + pp / ee) ;
903 return double(1) / (
gamma[0] - double(1)) ;
913 if ( ent < entLarge ) {
915 double log10H =
log10( ent ) ;
916 double log10HSmall =
log10( entSmall ) ;
917 double log10HLarge =
log10( entLarge ) ;
918 double dH = log10HLarge - log10HSmall ;
919 double uu = (log10H - log10HSmall) / dH ;
921 double mgam1skapgamSmall =
m0*(
gamma[i-1]-double(1))
923 double mgam1skapgamLarge =
m0*(
gamma[i]-double(1))
926 double nnSmall =
pow( mgam1skapgamSmall
928 double(1)/(
gamma[i-1]-
double(1)) ) ;
929 double nnLarge =
pow( mgam1skapgamLarge
931 double(1)/(
gamma[i]-
double(1)) ) ;
936 double eeSmall =
mu0[i-1] * nnSmall
937 + ppSmall / (
gamma[i-1] - double(1)) ;
938 double eeLarge =
mu0[i] * nnLarge
939 + ppLarge / (
gamma[i] - double(1)) ;
941 double log10PSmall =
log10( ppSmall ) ;
942 double log10PLarge =
log10( ppLarge ) ;
944 double dlpsdlhSmall = entSmall
945 * (double(1) + eeSmall / ppSmall) ;
946 double dlpsdlhLarge = entLarge
947 * (double(1) + eeLarge / ppLarge) ;
949 double log10PInterpolate = logp(log10PSmall, log10PLarge,
950 dlpsdlhSmall, dlpsdlhLarge,
953 double dlpsdlhInterpolate = dlpsdlh(log10PSmall, log10PLarge,
954 dlpsdlhSmall, dlpsdlhLarge,
957 pp =
pow(
double(10), log10PInterpolate) ;
959 ee = pp / ent * dlpsdlhInterpolate - pp ;
972 double(1)/(
gamma[i]-
double(1)) ) ;
976 ee =
mu0[i] * nn + pp / (
gamma[i] - double(1)) ;
978 return ent / (double(1) - (
mu0[i]/
m0) *
exp(-ent))
979 / (
gamma[i] -
double(1)) * (
double(1) + pp / ee) ;
995 for (
int l=0; l<
npeos-1; l++) {
1003 if ( ent >
double(0) ) {
1005 if ( ent < 1.e-13 ) {
1008 * ( double(1) + ent/double(2) + ent*ent/double(12) )
1009 / (
gamma[0] -
double(1)) ;
1014 return gamma[0] * ent / (double(1) -
exp(-ent))
1015 / (
gamma[0] -
double(1)) ;
1032 if ( ent < entLarge ) {
1034 double log10H =
log10( ent ) ;
1035 double log10HSmall =
log10( entSmall ) ;
1036 double log10HLarge =
log10( entLarge ) ;
1037 double dH = log10HLarge - log10HSmall ;
1038 double uu = (log10H - log10HSmall) / dH ;
1040 double mgam1skapgamSmall =
m0*(
gamma[i-1]-double(1))
1042 double mgam1skapgamLarge =
m0*(
gamma[i]-double(1))
1045 double nnSmall =
pow( mgam1skapgamSmall
1047 double(1)/(
gamma[i-1]-
double(1)) ) ;
1048 double nnLarge =
pow( mgam1skapgamLarge
1050 double(1)/(
gamma[i]-
double(1)) ) ;
1055 double eeSmall =
mu0[i-1] * nnSmall
1056 + ppSmall / (
gamma[i-1] - double(1)) ;
1057 double eeLarge =
mu0[i] * nnLarge
1058 + ppLarge / (
gamma[i] - double(1)) ;
1060 double log10PSmall =
log10( ppSmall ) ;
1061 double log10PLarge =
log10( ppLarge ) ;
1063 double dlpsdlhSmall = entSmall
1064 * (double(1) + eeSmall / ppSmall) ;
1065 double dlpsdlhLarge = entLarge
1066 * (double(1) + eeLarge / ppLarge) ;
1068 double dlpsdlhInterpolate = dlpsdlh(log10PSmall, log10PLarge,
1069 dlpsdlhSmall, dlpsdlhLarge,
1071 return dlpsdlhInterpolate ;
1076 return gamma[i] * ent / (double(1) - (
mu0[i]/
m0) *
exp(-ent))
1077 / (
gamma[i] -
double(1)) ;
1093 for (
int l=0; l<
npeos-1; l++) {
1109 if ( ent < entLarge ) {
1111 double log10H =
log10( ent ) ;
1112 double log10HSmall =
log10( entSmall ) ;
1113 double log10HLarge =
log10( entLarge ) ;
1115 double dlpsdlnbInterpolate = dlpsdlnb(log10HSmall, log10HLarge,
1119 return dlpsdlnbInterpolate ;
1136double logp(
double log10PressSmall,
double log10PressLarge,
1137 double dlpsdlhSmall,
double dlpsdlhLarge,
1138 double dx,
double u) {
1140 double resu = log10PressSmall * (double(2) * u * u * u
1141 - double(3) * u * u + double(1))
1142 + log10PressLarge * (
double(3) * u * u - double(2) * u * u * u)
1143 + dlpsdlhSmall * dx * (u * u * u -
double(2) * u * u + u)
1144 - dlpsdlhLarge * dx * (u * u - u * u * u) ;
1150double dlpsdlh(
double log10PressSmall,
double log10PressLarge,
1151 double dlpsdlhSmall,
double dlpsdlhLarge,
1152 double dx,
double u) {
1154 double resu = double(6) * (log10PressSmall - log10PressLarge)
1156 + dlpsdlhSmall * (double(3) * u * u - double(4) * u + double(1))
1157 + dlpsdlhLarge * (
double(3) * u * u - double(2) * u) ;
1163double dlpsdlnb(
double log10HSmall,
double log10HLarge,
1164 double dlpsdlnbSmall,
double dlpsdlnbLarge,
1167 double resu = log10H * (dlpsdlnbSmall - dlpsdlnbLarge)
1168 / (log10HSmall - log10HLarge)
1169 + (log10HSmall * dlpsdlnbLarge - log10HLarge * dlpsdlnbSmall)
1170 / (log10HSmall - log10HLarge) ;
Base class for a multiple polytropic equation of state.
virtual void sauve(FILE *) const
Save in a file.
const double & get_gamma(int n) const
Returns the adiabatic index .
int npeos
Number of polytropic equations of state.
double * decInc
Array (size npeos - 1) of the percentage which detemines the terminating enthalpy for lower density a...
double m0
Individual particule mass [unit: ].
void operator=(const Eos_multi_poly &)
Assignment to another Eos_multi_poly.
double * nbCrit
Array (size npeos - 1) of the number density at which the polytropic EOS changes its index and consta...
const double & get_kappa(int n) const
Returns the pressure coefficient [unit: ], where and .
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.
double * mu0
Array (size: npeos) of the relativistic chemical potential at zero pressure [unit: ,...
virtual double nbar_ent_p(double ent, const Param *par=0x0) const
Computes the baryon density from the log-enthalpy.
double * gamma
Array (size: npeos) of adiabatic index .
const int & get_npeos() const
Returns the number of polytropes npeos.
double * logRho
Array (size: npeos - 1) of the exponent of fiducial densities.
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.
double logP1
Exponent of the pressure at the fiducial density .
virtual bool operator!=(const Eos &) const
Comparison operator (difference)
void set_auxiliary()
Computes the auxiliary quantities.
virtual double der_press_nbar_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
virtual int identify() const
Returns a number to identify the sub-classe of Eos the object belongs to.
virtual ostream & operator>>(ostream &) const
Operator >>
double * kappa
Array (size: npeos) of pressure coefficient [unit: ], where and .
double kappa0
Pressure coefficient for the crust [unit: ].
virtual ~Eos_multi_poly()
Destructor.
double * entCrit
Array (size npeos - 1) of the critical enthalpy at which the polytropic EOS changes its index and con...
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
Read/write kappa.
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).
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.
Cmp log10(const Cmp &)
Basis 10 logarithm.
Cmp exp(const Cmp &)
Exponential.
Cmp pow(const Cmp &, int)
Power .
Cmp log(const Cmp &)
Neperian logarithm.
int fread_be(int *aa, int size, int nb, FILE *fich)
Reads integer(s) from a binary file according to the big endian convention.
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.
Standard units of space, time and mass.