30 char eos_bf_tabul_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Eos/eos_bf_tabul.C,v 1.3 2015/06/11 14:41:59 a_sourie Exp $" ;
57 #include "eos_bifluid.h" 61 #include "utilitaires.h" 68 void interpol_linear(
const Tbl& xtab,
const Tbl& ytab,
69 double x,
int& i,
double& y) ;
71 void interpol_herm(
const Tbl& xtab,
const Tbl& ytab,
const Tbl& dytab,
72 double x,
int& i,
double& y,
double& dy) ;
74 void interpol_mixed_3d(
const Tbl&,
const Tbl&,
const Tbl&,
const Tbl&,
75 const Tbl&,
const Tbl&,
const Tbl&,
76 double,
double,
double,
double&,
double&,
double&) ;
79 void interpol_mixed_3d_new(
double m1,
double m2,
const Tbl& xtab,
const Tbl& ytab,
const Tbl& ztab,
80 const Tbl& ftab,
const Tbl& dfdytab,
const Tbl& dfdztab,
const Tbl& d2fdydztab,
81 const Tbl& dlpsddelta_car,
const Tbl& d2lpsdlent1ddelta_car,
const Tbl& d2lpsdlent2ddelta_car,
82 const Tbl& mu2_P,
const Tbl& n_p_P,
const Tbl& press_P,
83 const Tbl& mu1_N,
const Tbl& n_n_N,
const Tbl& press_N,
84 const Tbl& delta_car_n0,
const Tbl& mu1_n0,
const Tbl& mu2_n0,
85 const Tbl& delta_car_p0,
const Tbl& mu1_p0,
const Tbl& mu2_p0,
86 double x,
double y,
double z,
double& f,
double& dfdy,
double& dfdz,
double& alpha) ;
88 void interpolation_brutale(
double x,
double y,
double z,
90 double delta111,
double delta211,
91 double mu1_111,
double mu1_121,
double mu2_111,
double mu2_112,
double mu1_211,
double mu1_221,
double mu2_211,
double mu2_212,
92 double p_111,
double p_121,
double p_112,
double p_122,
double p_211,
double p_221,
double p_212,
double p_222,
93 double n1_111,
double n1_121,
double n1_112,
double n1_122,
double n1_211,
double n1_221,
double n1_212,
double n1_222,
94 double n2_111,
double n2_121,
double n2_112,
double n2_122,
double n2_211,
double n2_221,
double n2_212,
double n2_222,
95 double cross_111,
double cross_121,
double cross_112,
double cross_122,
double cross_211,
double cross_221,
double cross_212,
double cross_222,
96 double& f,
double& dfdy,
double& dfdz) ;
98 void interpolation_brutale_mod(
double x,
double y,
double z,
100 double delta111,
double delta211,
101 double mu1_111,
double mu1_121,
double mu2_111,
double mu2_112,
double mu1_211,
double mu1_221,
double mu2_211,
double mu2_212,
102 double p_111,
double p_121,
double p_112,
double p_122,
double p_211,
double p_221,
double p_212,
double p_222,
103 double n1_111,
double n1_121,
double n1_112,
double n1_122,
double n1_211,
double n1_221,
double n1_212,
double n1_222,
104 double n2_111,
double n2_121,
double n2_112,
double n2_122,
double n2_211,
double n2_221,
double n2_212,
double n2_222,
105 double& f,
double& dfdy,
double& dfdz) ;
116 void interpol_mixed_3d_mod(
const Tbl&,
const Tbl&,
const Tbl&,
const Tbl&,
117 const Tbl&,
const Tbl&,
118 double,
double,
double,
double&,
double&,
double&) ;
129 const char* path,
double mass1,
double mass2) :
142 double mass1,
double mass2) :
155 char tmp_string[160] ;
156 fread(tmp_string,
sizeof(
char), 160, fich) ;
200 delete delta_car_n0 ;
203 delete delta_car_p0 ;
242 delta_car_n0= eosi.delta_car_n0 ;
243 mu1_n0= eosi.mu1_n0 ;
244 mu2_n0= eosi.mu2_n0 ;
245 delta_car_p0= eosi.delta_car_p0 ;
246 mu1_p0 = eosi.mu1_p0;
247 mu2_p0 = eosi.mu2_p0 ;
250 press_N = eosi.press_N;
253 press_P = eosi.press_P;
266 char tmp_string[160] ;
268 fwrite(tmp_string,
sizeof(
char), 160, fich) ;
276 "EOS of class Eos_bf_tabul : tabulated EOS depending on three variables (relative velocity and enthalpies of neutrons and protons)." 278 ost <<
"Table read from " <<
tablename << endl ;
293 cout <<
"The second EOS is not of type Eos_bf_tabul !" << endl ;
302 "The two Eos_bf_tabul have different names and not refer to the same tables! " << endl ;
326 double m_b_si = 1.66E-27;
335 cerr <<
"Eos_bf_tabul::read_table(): " << endl ;
336 cerr <<
"Problem in opening the EOS file!" << endl ;
337 cerr <<
"While trying to open " <<
tablename << endl ;
338 cerr <<
"Aborting..." << endl ;
342 fich.ignore(1000,
'\n') ;
346 for (
int i=0; i<5; i++) {
347 fich.ignore(1000,
'\n') ;
351 fich >> nbp ; fich.ignore(1000,
'\n') ;
354 int n_delta, n_mu1, n_mu2;
355 fich >> n_delta; fich.ignore(1000,
'\n') ;
356 fich >> n_mu1; fich.ignore(1000,
'\n') ;
357 fich >> n_mu2; fich.ignore(1000,
'\n') ;
363 cerr <<
"Eos_bf_tabul::read_table(): " << endl ;
364 cerr <<
"Wrong value for the number of lines!" << endl ;
365 cerr <<
"nbp = " << nbp << endl ;
366 cerr <<
"Aborting..." << endl ;
370 cerr <<
"Eos_bf_tabul::read_table(): " << endl ;
371 cerr <<
"Wrong value for the number of values of delta_car!" << endl ;
372 cerr <<
"n_delta = " << n_delta << endl ;
373 cerr <<
"Aborting..." << endl ;
377 cerr <<
"Eos_bf_tabul::read_table(): " << endl ;
378 cerr <<
"Wrong value for the number of values of mu_n!" << endl ;
379 cerr <<
"n_mu1 = " << n_mu1 << endl ;
380 cerr <<
"Aborting..." << endl ;
384 cerr <<
"Eos_bf_tabul::read_table(): " << endl ;
385 cerr <<
"Wrong value for the number of values of mu_p!" << endl ;
386 cerr <<
"n_mu2 = " << n_mu2 << endl ;
387 cerr <<
"Aborting..." << endl ;
390 if (n_mu2*n_mu1*n_delta != nbp ) {
391 cerr <<
"Eos_bf_tabul::read_table(): " << endl ;
392 cerr <<
"Wrong value for the number of lines!" << endl ;
393 cerr <<
"Aborting..." << endl ;
397 for (
int i=0; i<3; i++) {
398 fich.ignore(1000,
'\n') ;
404 logp =
new Tbl(n_delta, n_mu1, n_mu2) ;
433 double mu1_MeV, mu2_MeV, n1_fm3, n2_fm3;
434 double Knp_Mev_2, press_MeV_fm3;
435 double d2press_dmu1dmu2_MeV_fm3, dn1_ddelta_car_fm3, dn2_ddelta_car_fm3;
438 double delta_car_adim, mu1_adim, mu2_adim, n1_01fm3, n2_01fm3, Knp_adim ;
439 double press_adim, dpress_ddelta_car_adim, dn1_ddelta_car_adim, dn2_ddelta_car_adim ;
440 double d2press_dmu1dmu2_adim;
441 double hbarc = 197.33 ;
442 double hbarc3 = hbarc * hbarc * hbarc ;
445 for (
int j=0 ; j < n_delta ; j++) {
446 for (
int k=0 ; k < n_mu1 ; k++) {
447 for (
int l=0 ; l < n_mu2 ; l++) {
449 fich >> delta_car_adim ;
451 mu1_adim = mu1_MeV * mev_si / (m_b_si * v_unit * v_unit ) ;
453 mu2_adim = mu2_MeV * mev_si / (m_b_si * v_unit * v_unit ) ;
455 n1_01fm3 = 10. * n1_fm3 ;
457 n2_01fm3 = 10. * n2_fm3 ;
459 Knp_adim = Knp_Mev_2 / ( m_b_si * v_unit * v_unit *10. )
460 * (mev_si *hbarc3 ) ;
461 dpress_ddelta_car_adim = - Knp_adim * n1_01fm3 * n2_01fm3
462 *
pow( 1.-delta_car_adim, -1.5) / 2. ;
464 fich >> press_MeV_fm3 ;
465 press_adim = press_MeV_fm3 * mevpfm3 ;
466 fich >> d2press_dmu1dmu2_MeV_fm3 ;
467 d2press_dmu1dmu2_adim = d2press_dmu1dmu2_MeV_fm3
468 * (10. * m_b_si * v_unit * v_unit ) / mev_si ;
469 fich >> dn1_ddelta_car_fm3 ;
470 dn1_ddelta_car_adim = dn1_ddelta_car_fm3 * 10. ;
471 fich >> dn2_ddelta_car_fm3 ;
472 dn2_ddelta_car_adim = dn2_ddelta_car_fm3 * 10. ;
477 fich.ignore(1000,
'\n') ;
505 if (press_adim<=0) {press_adim = 0. ;}
506 logp->
set(j,k,l) = press_adim ;
512 if ((n1_01fm3 < 1e-16) && (n2_01fm3 <1e-16 )) {
513 d2press_dmu1dmu2_adim = 0. ;
514 dpress_ddelta_car_adim = 0. ;
515 dn1_ddelta_car_adim =0. ;
516 dn2_ddelta_car_adim = 0. ;
526 fich.ignore(1000,
'\n') ;
530 fich.ignore(1000,
'\n') ;
552 ent1_max = (*logent1)(0, n_mu1-1, 0) ;
554 ent2_max = (*logent2)(0, 0, n_mu2-1);
577 fichN.open(
"eos_bf_test_1_fluide_N.d") ;
580 cerr <<
"Eos_bf_tabul::read_table(): " << endl ;
581 cerr <<
"Problem in opening the EOS file!" << endl ;
582 cerr <<
"While trying to open " <<
tablename << endl ;
583 cerr <<
"Aborting..." << endl ;
587 fichN.ignore(1000,
'\n') ;
589 fichN.ignore(1000,
'\n') ;
591 for (
int i=0; i<5; i++) {
592 fichN.ignore(1000,
'\n') ;
596 fichN >> nbp_N ; fichN.ignore(1000,
'\n') ;
597 cout<<
"nbp_N = " << nbp_N ;
599 fichN >> n_mu1_N;fichN.ignore(1000,
'\n') ;
604 cerr <<
"Eos_bf_tabul::read_table(): " << endl ;
605 cerr <<
"Wrong value for the number of lines!" << endl ;
606 cerr <<
"nbp = " << nbp << endl ;
607 cerr <<
"Aborting..." << endl ;
611 cerr <<
"Eos_bf_tabul::read_table(): " << endl ;
612 cerr <<
"Wrong value for the number of values of mu_n!" << endl ;
613 cerr <<
"n_mu1 = " << n_mu1 << endl ;
614 cerr <<
"Aborting..." << endl ;
617 for (
int i=0; i<3; i++) {
618 fichN.ignore(1000,
'\n') ;
621 mu1_N =
new Tbl(n_mu1_N) ;
622 n_n_N =
new Tbl(n_mu1_N) ;
623 press_N =
new Tbl(n_mu1_N) ;
625 mu1_N ->set_etat_qcq() ;
626 n_n_N->set_etat_qcq() ;
627 press_N ->set_etat_qcq() ;
630 double mu1_MeV_N, n1_fm3_N,press_MeV_fm3_N;
633 double mu1_adimN, n1_01fm3N,press_adimN;
635 for (
int k=0 ; k < n_mu1_N ; k++) {
638 mu1_adimN = mu1_MeV_N * mev_si / (m_b_si * v_unit * v_unit ) ;
640 n1_01fm3N = 10. * n1_fm3_N ;
641 fichN >> press_MeV_fm3_N;
642 press_adimN = press_MeV_fm3_N * mevpfm3 ;
643 fichN.ignore(1000,
'\n') ;
647 if ( (n1_01fm3N<0) || (press_adimN < 0)){
648 cout <<
"Eos_tabul::read_table(): " << endl ;
649 cout <<
"Negative value in table!" << endl ;
650 cout <<
"n_neutrons = " << n1_01fm3N <<
651 ", p = " << press_adimN <<
", "<< endl ;
652 cout <<
"Aborting..." << endl ;
656 mu1_N ->set(k) = mu1_adimN ;
657 n_n_N->set(k) = n1_01fm3N ;
658 press_N ->set(k) = press_adimN ;
680 fichP.open(
"eos_bf_test_1_fluide_P.d") ;
683 cerr <<
"Eos_bf_tabul::read_table(): " << endl ;
684 cerr <<
"Problem in opening the EOS file!" << endl ;
685 cerr <<
"While trying to open " <<
tablename << endl ;
686 cerr <<
"Aborting..." << endl ;
690 fichP.ignore(1000,
'\n') ;
692 fichP.ignore(1000,
'\n') ;
694 for (
int i=0; i<5; i++) {
695 fichP.ignore(1000,
'\n') ;
699 fichP >> nbp_P ; fichP.ignore(1000,
'\n') ;
702 fichP >> n_mu2_P;fichP.ignore(1000,
'\n') ;
705 cerr <<
"Eos_bf_tabul::read_table(): " << endl ;
706 cerr <<
"Wrong value for the number of lines!" << endl ;
707 cerr <<
"nbp = " << nbp << endl ;
708 cerr <<
"Aborting..." << endl ;
712 cerr <<
"Eos_bf_tabul::read_table(): " << endl ;
713 cerr <<
"Wrong value for the number of values of mu_p!" << endl ;
714 cerr <<
"n_mu2 = " << n_mu2 << endl ;
715 cerr <<
"Aborting..." << endl ;
719 for (
int i=0; i<3; i++) {
720 fichP.ignore(1000,
'\n') ;
723 mu2_P =
new Tbl(n_mu2_P) ;
724 n_p_P =
new Tbl(n_mu2_P) ;
725 press_P =
new Tbl(n_mu2_P) ;
727 mu2_P ->set_etat_qcq() ;
728 n_p_P->set_etat_qcq() ;
729 press_P ->set_etat_qcq() ;
732 double mu2_MeV_P, n2_fm3_P,press_MeV_fm3_P;
735 double mu2_adimP, n2_01fm3P, press_adimP;
737 for (
int l=0 ; l < n_mu2_P ; l++) {
740 mu2_adimP = mu2_MeV_P * mev_si / (m_b_si * v_unit * v_unit ) ;
742 n2_01fm3P = 10. * n2_fm3_P ;
743 fichP >> press_MeV_fm3_P;
744 press_adimP = press_MeV_fm3_P * mevpfm3 ;
748 fichP.ignore(1000,
'\n') ;
749 if ( (n2_01fm3P<0) || (press_adimP < 0)){
750 cout <<
"Eos_tabul::read_table(): " << endl ;
751 cout <<
"Pegative value in table!" << endl ;
752 cout <<
", n_protons " << n2_01fm3P <<
753 ", p = " << press_adimP <<
", "<< endl ;
754 cout <<
"Aborting..." << endl ;
758 mu2_P ->set(l) = mu2_adimP ;
759 n_p_P->set(l) = n2_01fm3P ;
760 press_P ->set(l) = press_adimP ;
772 fich1.open(
"np=0.dat") ;
775 cerr <<
"Eos_bf_tabul::read_table(): " << endl ;
776 cerr <<
"Problem in opening the EOS file!" << endl ;
777 cerr <<
"While trying to open " <<
tablename << endl ;
778 cerr <<
"Aborting..." << endl ;
781 int n_delta_n0, n_mu1_n0;
782 fich1 >> n_delta_n0;fich1.ignore(1000,
'\n') ;
783 fich1 >> n_mu1_n0;fich1.ignore(1000,
'\n') ;
784 fich1.ignore(1000,
'\n') ;
786 delta_car_n0 =
new Tbl(n_delta_n0, n_mu1_n0) ;
787 mu1_n0 =
new Tbl(n_delta_n0, n_mu1_n0) ;
788 mu2_n0 =
new Tbl(n_delta_n0, n_mu1_n0) ;
790 delta_car_n0 ->set_etat_qcq() ;
791 mu1_n0->set_etat_qcq() ;
792 mu2_n0->set_etat_qcq() ;
794 double delta_car_nn0, mu1_MeV_nn0, mu2_MeV_nn0;
796 for (
int o = 0; o < n_delta_n0 ; o++ ) {
797 for (
int p = 0 ; p < n_mu1_n0 ; p++) {
799 fich1 >> delta_car_nn0 ;
800 fich1 >> mu1_MeV_nn0 ;
801 fich1 >> mu2_MeV_nn0 ;
803 fich1.ignore(1000,
'\n') ;
805 delta_car_n0 ->set(o,p) = delta_car_nn0;
806 mu1_n0 ->set(o,p) = mu1_MeV_nn0 * mev_si / (m_b_si * v_unit * v_unit ) ;
807 mu2_n0 ->set(o,p) = mu2_MeV_nn0 * mev_si / (m_b_si * v_unit * v_unit ) ;
810 fich1.ignore(1000,
'\n') ;
820 fich2.open(
"nn=0.dat") ;
824 cerr <<
"Eos_bf_tabul::read_table(): " << endl ;
825 cerr <<
"Problem in opening the EOS file!" << endl ;
826 cerr <<
"While trying to open " <<
tablename << endl ;
827 cerr <<
"Aborting..." << endl ;
830 int n_delta_p0, n_mu2_p0;
831 fich2 >> n_delta_p0;fich2.ignore(1000,
'\n') ;
832 fich2 >> n_mu2_p0;fich2.ignore(1000,
'\n') ;
833 fich2.ignore(1000,
'\n') ;
835 delta_car_p0 =
new Tbl(n_delta_p0, n_mu2_p0) ;
836 mu1_p0 =
new Tbl(n_delta_p0, n_mu2_p0) ;
837 mu2_p0 =
new Tbl(n_delta_p0, n_mu2_p0) ;
839 delta_car_p0 ->set_etat_qcq() ;
840 mu1_p0->set_etat_qcq() ;
841 mu2_p0 ->set_etat_qcq() ;
843 double delta_car_np0, mu1_MeV_np0, mu2_MeV_np0;
845 for (
int o = 0; o < n_delta_p0 ; o++ ) {
846 for (
int p = 0 ; p < n_mu2_p0 ; p++) {
848 fich2 >> delta_car_np0 ;
849 fich2 >> mu1_MeV_np0 ;
850 fich2 >> mu2_MeV_np0 ;
852 fich2.ignore(1000,
'\n') ;
854 delta_car_p0 ->set(o,p) = delta_car_np0;
855 mu1_p0 ->set(o,p) = mu1_MeV_np0 * mev_si / (m_b_si * v_unit * v_unit ) ;
856 mu2_p0 ->set(o,p) = mu2_MeV_np0 * mev_si / (m_b_si * v_unit * v_unit ) ;
860 fich2.ignore(1000,
'\n') ;
877 const Cmp& delta2,
Cmp& nbar1,
Cmp& nbar2,
879 int nzet,
int l_min)
const {
881 assert(ent1.
get_etat() != ETATNONDEF) ;
882 assert(ent2.
get_etat() != ETATNONDEF) ;
883 assert(delta2.
get_etat() != ETATNONDEF) ;
886 assert(mp == ent2.
get_mp()) ;
887 assert(mp == delta2.
get_mp()) ;
888 assert(mp == ener.
get_mp()) ;
907 const Mg3d* mg = mp->get_mg() ;
914 nbar1.
annule(0, l_min-1) ;
915 nbar2.
annule(0, l_min-1) ;
917 press.
annule(0, l_min-1) ;
923 if (l_min + nzet < nz) {
924 nbar1.
annule(l_min + nzet, nz - 1) ;
925 nbar2.
annule(l_min + nzet, nz - 1) ;
926 ener.
annule(l_min + nzet, nz - 1) ;
927 press.
annule(l_min + nzet, nz - 1) ;
928 K_nn.
annule(l_min + nzet, nz - 1) ;
929 K_np.
annule(l_min + nzet, nz - 1) ;
930 K_pp.
annule(l_min + nzet, nz - 1) ;
933 int np0 = mp->get_mg()->get_np(0) ;
934 int nt0 = mp->get_mg()->get_nt(0) ;
935 for (
int l=l_min; l<l_min+nzet; l++) {
936 assert(mp->get_mg()->get_np(l) == np0) ;
937 assert(mp->get_mg()->get_nt(l) == nt0) ;
941 for (
int k=0; k<np0; k++) {
942 for (
int j=0; j<nt0; j++) {
944 for (
int l=l_min; l< l_min + nzet; l++) {
945 for (
int i=0; i<mp->get_mg()->get_nr(l); i++) {
948 xx1 = ent1(l,k,j,i) ;
949 xx2 = ent2(l,k,j,i) ;
950 xx = delta2(l,k,j,i) ;
1030 cout <<
"Eos_bf_tabul::calcule_tout : delta2 < delta_car_min !" 1035 cout <<
"Eos_bf_tabul::calcule_tout : delta2 > delta_car_max !" 1040 cout <<
"Eos_bf_tabul::calcule_tout : ent1 > ent1_max !" << endl ;
1044 cout <<
"Eos_bf_tabul::calcule_tout : ent2 > ent2_max !" << endl ;
1050 double pressure = 0 ;
1056 double mu1 =
m_1 *
exp(xx1);
1057 double mu2 =
m_2 *
exp(xx2);
1059 if ( (
exp(xx1) < 1.) && (
exp(xx2) < 1.) ) {
1073 double dpsdent1_interpol ;
1074 double dpsdent2_interpol ;
1075 double alpha_interpol ;
1082 *mu2_P, *n_p_P, *press_P,
1083 *mu1_N, *n_n_N, *press_N,
1084 *delta_car_n0, *mu1_n0, *mu2_n0,
1085 *delta_car_p0, *mu1_p0, *mu2_p0,
1087 p_interpol, dpsdent1_interpol, dpsdent2_interpol, alpha_interpol) ;
1093 n1 = dpsdent1_interpol ;
1094 n2 = dpsdent2_interpol ;
1095 pressure = p_interpol;
1096 alpha = alpha_interpol ;
1101 cout <<
"Eos_bf_tabul::nbar_ent_p : nbar1<0 !" << endl ;
1106 cout <<
"Eos_bf_tabul::nbar_ent_p : nbar2<0 !" << endl ;
1113 if (pressure < 0 ) {
1114 cout <<
"Eos_bf_tabul::press_ent_p : pressure < 0 !" << endl ;
1122 K11 = mu1 / n1 - double(2) * alpha * ( 1. - xx) / ( n1 * n1 ) ;
1129 K22 = mu2 / n2 - double(2) * alpha * ( 1. - xx) / ( n2 * n2 ) ;
1134 if ((n1 <= 0.) || (n2 <= 0.) ) {
1141 K12 = double(2) * alpha *
pow(1.-xx, 1.5)/ ( n1 * n2 );
1146 nbar1.
set(l,k,j,i) = n1 ;
1147 nbar2.
set(l,k,j,i) = n2 ;
1148 press.
set(l,k,j,i) = pressure ;
1149 ener.
set(l,k,j,i) = mu1 * n1 + mu2 * n2 - pressure ;
1150 K_nn.
set(l,k,j,i) = K11 ;
1151 K_np.
set(l,k,j,i) = K12;
1152 K_pp.
set(l,k,j,i) = K22 ;
1182 const double delta2,
double& nbar1,
1183 double& nbar2)
const {
1185 bool one_fluid =
false;
1188 cout <<
"Eos_bf_tabul::nbar_ent_p : delta2 < delta_car_min !" << endl ;
1192 cout <<
"Eos_bf_tabul::nbar_ent_p : delta2 > delta_car_max !" << endl ;
1196 cout <<
"Eos_bf_tabul::nbar_ent_p : ent1 > ent1_max !" << endl ;
1200 cout <<
"Eos_bf_tabul::nbar_ent_p : ent2 > ent2_max !" << endl ;
1204 if ( (
exp(ent1) < 1.) && (
exp(ent2) < 1.) ) {
1210 double mu1 =
m_1 *
exp(ent1);
1211 double mu2 =
m_2 *
exp(ent2);
1214 double dpsdent1_interpol ;
1215 double dpsdent2_interpol ;
1221 *mu2_P, *n_p_P, *press_P,
1222 *mu1_N, *n_n_N, *press_N,
1223 *delta_car_n0, *mu1_n0, *mu2_n0,
1224 *delta_car_p0, *mu1_p0, *mu2_p0,
1226 p_interpol, dpsdent1_interpol, dpsdent2_interpol, alpha) ;
1228 nbar1 = dpsdent1_interpol ;
1229 nbar2 = dpsdent2_interpol ;
1234 cout <<
"Eos_bf_tabul::nbar_ent_p : nbar1<0 !" << endl ;
1239 cout <<
"Eos_bf_tabul::nbar_ent_p : nbar2<0 !" << endl ;
1245 one_fluid = ((nbar1 <= 0.)||(nbar2 <= 0.)) ;
1257 double pressN_interpol ;
1258 double n_n_N_interpol ;
1259 double mu1 =
m_1 *
exp(ent1);
1262 if (
exp(ent1) < 1. ) {
1263 n_n_N_interpol = 0. ;
1266 interpol_herm( *mu1_N, *press_N, *n_n_N,
1267 mu1, i, pressN_interpol , n_n_N_interpol) ;
1269 return n_n_N_interpol ;
1276 double pressP_interpol ;
1277 double n_p_P_interpol ;
1278 double mu2 =
m_2 *
exp(ent2);
1280 if (
exp(ent2) < 1. ) {
1281 n_p_P_interpol = 0. ;
1284 interpol_herm( *mu2_P, *press_P, *n_p_P,
1285 mu2, i, pressP_interpol , n_p_P_interpol) ;
1287 return n_p_P_interpol ;
1294 const double delta2)
const{
1298 return nbar1 + nbar2 + delta2;
1305 const double delta2)
const {
1309 return nbar1 + nbar2 + delta2;
1319 cout <<
"Eos_bf_tabul::press_ent_p : delta2 < delta_car_min !" << endl ;
1323 cout <<
"Eos_bf_tabul::press_ent_p : delta2 > delta_car_max !" << endl ;
1327 cout <<
"Eos_bf_tabul::press_ent_p : ent1 > ent1_max !" << endl ;
1331 cout <<
"Eos_bf_tabul::press_ent_p : ent2 > ent2_max !" << endl ;
1337 if ( (
exp(ent1) < 1.) && (
exp(ent2) < 1.)) {
1344 double mu1 =
m_1 *
exp(ent1);
1345 double mu2 =
m_2 *
exp(ent2);
1348 double dpsdent1_interpol ;
1349 double dpsdent2_interpol ;
1355 *mu2_P, *n_p_P, *press_P,
1356 *mu1_N, *n_n_N, *press_N,
1357 *delta_car_n0, *mu1_n0, *mu2_n0,
1358 *delta_car_p0, *mu1_p0, *mu2_p0,
1360 p_interpol, dpsdent1_interpol, dpsdent2_interpol, alpha) ;
1362 pressure = p_interpol;
1366 if (pressure < 0 ) {
1367 cout <<
"Eos_bf_tabul::press_ent_p : pressure < 0 !" << endl ;
1379 cout <<
"Eos_bf_tabul::press_ent_p : ent1 > ent1_max !" << endl ;
1385 double pressN_interpol ;
1386 double n_n_N_interpol ;
1388 double mu1 =
m_1 *
exp(ent1);
1390 interpol_herm( *mu1_N, *press_N, *n_n_N,
1391 mu1, i, pressN_interpol , n_n_N_interpol) ;
1393 pressure = pressN_interpol ;
1405 cout <<
"Eos_bf_tabul::press_ent_p : ent1 > ent1_max !" << endl ;
1410 double pressP_interpol ;
1411 double n_p_P_interpol ;
1412 double mu2 =
m_2*
exp(ent2);
1414 interpol_herm( *mu2_P, *press_P, *n_p_P,
1415 mu2, i, pressP_interpol , n_p_P_interpol) ;
1417 pressure = pressP_interpol ;
1428 const double delta2)
const {
1431 if ( (
exp(ent1) < 1.) && (
exp(ent2) < 1.)) {
1470 double mu1 =
m_1 *
exp(ent1);
1471 double mu2 =
m_2 *
exp(ent2);
1474 double dpsdent1_interpol ;
1475 double dpsdent2_interpol ;
1481 *mu2_P, *n_p_P, *press_P,
1482 *mu1_N, *n_n_N, *press_N,
1483 *delta_car_n0, *mu1_n0, *mu2_n0,
1484 *delta_car_p0, *mu1_p0, *mu2_p0,
1485 delta2, mu1, mu2, p_interpol, dpsdent1_interpol, dpsdent2_interpol, alpha) ;
1487 energy = mu1 * dpsdent1_interpol + mu2 * dpsdent2_interpol - p_interpol ;
1492 cout <<
"Eos_bf_tabul::ener_ent_p : energy < 0 !" << endl ;
1504 const double delta2)
const {
1507 cout <<
"Eos_bf_tabul::alpha_ent_p : delta2 < delta_car_min !" 1512 cout <<
"Eos_bf_tabul::alpha_ent_p : delta2 > delta_car_max !" 1517 cout <<
"Eos_bf_tabul::alpha_ent_p : ent1 > ent1_max !" << endl ;
1521 cout <<
"Eos_bf_tabul::alpha_ent_p : ent2 > ent2_max !" << endl ;
1527 if ((
exp(ent1) <= 1.) && (
exp(ent2) <= 1.) ) {
1533 double mu1 =
m_1 *
exp(ent1);
1534 double mu2 =
m_2 *
exp(ent2);
1537 double dpsdent1_interpol ;
1538 double dpsdent2_interpol ;
1543 *mu2_P, *n_p_P, *press_P,
1544 *mu1_N, *n_n_N, *press_N,
1545 *delta_car_n0, *mu1_n0, *mu2_n0,
1546 *delta_car_p0, *mu1_p0, *mu2_p0,
1547 delta2, mu1, mu2, p_interpol, dpsdent1_interpol, dpsdent2_interpol, alpha) ;
1566 double mu_1 =
m_1 *
exp(ent1);
1570 if ((
exp(ent1) <= 1.) && (
exp(ent2) <= 1.) ){
1596 xx = mu_1 / nbar1 - double(2) * alpha * ( 1. - delta2) / ( nbar1 * nbar1 ) ;
1613 double mu_2 =
m_2 *
exp (ent2) ;
1617 if ((
exp(ent1) <= 1.) && (
exp(ent2) <= 1.) ){
1643 xx = mu_2 / nbar2 - double(2) * alpha * ( 1. - delta2) / ( nbar2 * nbar2 ) ;
1661 if ((
exp(ent1) <= 1.) && (
exp(ent2) <= 1.) ){
1669 if ((nbar1 <= 0.) || (nbar2 <= 0.) ) {
1675 xx = double(2) * alpha *
pow(1.-delta2, 1.5)/ ( nbar1 * nbar2 );
const Map * get_mp() const
Returns the mapping.
virtual void sauve(FILE *) const
Save in a file.
double m_1
Individual particle mass [unit: ].
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Cmp exp(const Cmp &)
Exponential.
Tbl * d3lpsdlent1dlent2ddelta_car
if necessary for the interpolation to find alpha (derivee seconde croisee) ie, if it's possible to ca...
Tbl * logent2
Table of where .
void annule(int l)
Sets the Cmp to zero in a given domain.
Standard units of space, time and mass.
Tbl * d2lpsdlent2ddelta_car
Table of .
Equation of state base class.
virtual int identify() const =0
Returns a number to identify the sub-classe of Eos_bifluid the object belongs to. ...
virtual ostream & operator>>(ostream &) const
Operator >>
int get_etat() const
Returns the logical state.
double & set(int i)
Read/write of a particular element (index i) (1D case)
Base class for coordinate mappings.
double ent1_min
Lower boundary of the log-enthalpy interval (fluid 1 = n)
virtual double nbar_ent_p1(const double ent1) const
Computes baryon density out of the log-enthalpy asuming that only fluid 1 is present.
virtual double get_K11(const double delta2, const double ent1, const double ent2) const
Computes the derivative of the energy with respect to (baryonic density 1) .
virtual int identify() const
Returns a number to identify the sub-classe of Eos the object belongs to.
virtual double press_ent_p2(const double ent2) const
Computes the pressure from the baryonic log-enthalpies assuming that only fluid 2 is present...
Tbl * dlpsdlent1
Table of .
void operator=(const Eos_bifluid &)
Assignment to another Eos_bifluid.
virtual double ener_ent_p(const double ent1, const double ent2, const double delta_car) const
Computes the total energy density from the baryonic log-enthalpies and the relative velocity...
virtual double alpha_ent_p(const double ent1, const double ent2, const double delta_car) const
Computes alpha, the derivative of the total energy density with respect to from the baryonic log-ent...
Tbl * logent1
Table of where .
Tbl * d2lpsdlent1ddelta_car
Table of .
2-fluids equation of state base class.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
double m_2
Individual particle mass [unit: ].
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Tbl * dlpsddelta_car
Table of .
virtual double nbar_ent_p2(const double ent2) const
Computes baryon density out of the log-enthalpy assuming that only fluid 2 is present.
void calcule_interpol(const Cmp &ent1, const Cmp &ent2, const Cmp &delta2, Cmp &nbar1, Cmp &nbar2, Cmp &ener, Cmp &press, Cmp &K_nn, Cmp &K_np, Cmp &K_pp, int nzet, int l_min=0) const
General computational method for Cmp 's, it computes both baryon densities, energy and pressure profi...
double ent2_max
Upper boundary of the log-enthalpy interval (fluid 2 = p)
double delta_car_min
Lower boundary of the relative velocity interval –> 0 ?
Tbl * d2lpsdlent1dlent1
Table of .
virtual bool operator==(const Eos_bifluid &) const
Comparison operator (egality)
Polytropic equation of state (relativistic case).
virtual double press_ent_p1(const double ent1) const
Computes the pressure from the baryonic log-enthalpies asuming that only fluid 1 is present...
int get_nzone() const
Returns the number of domains.
virtual double get_K22(const double delta2, const double ent1, const double ent2) const
Computes the derivative of the energy/(baryonic density 2) .
double ent1_max
Upper boundary of the log-enthalpy interval (fluid 1 = n)
virtual void sauve(FILE *) const
Save in a file.
virtual ~Eos_bf_tabul()
Destructor.
Cmp pow(const Cmp &, int)
Power .
void read_table()
Reads the file containing the table and initializes the arrays logent1, logent2, delta_car, logp, dlpsdlent1, dlpsdlent2, d2lpsdlent1dlent2, dlpsddelta_car, d2lpsdlent1ddelta_car, d2lpsdlent2ddelta_car, d3lpsdlent1dlent2ddelta_car.
virtual bool operator!=(const Eos_bifluid &) const
Comparison operator (difference)
virtual double press_nbar_p(const double nbar1, const double nbar2, const double delta2) const
Computes the pressure from the baryonic densities and the relative velocity.
void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
Tbl & set(int l)
Read/write of the value in a given domain.
virtual double ener_nbar_p(const double nbar1, const double nbar2, const double delta2) const
Computes the total energy density from the baryonic densities and the relative velocity.
virtual double press_ent_p(const double ent1, const double ent2, const double delta_car) const
Computes the pressure from the baryonic log-enthalpies and the relative velocity. ...
void c_est_pas_fait(const char *)
Helpful function to say something is not implemented yet.
Tbl * dlpsdlent2
Table of .
Tbl * delta_car
Table of .
Tbl * d2lpsdlent2dlent2
Table of .
Class for a two-fluid (tabulated) equation of state.
virtual double get_K12(const double delta2, const double ent1, const double ent2) const
Computes the derivative of the energy with respect to .
Tbl * d2lpsdlent1dlent2
Table of .
double delta_car_max
Upper boundary of the relative velocity interval –> 1 ?
void operator=(const Eos_bf_tabul &)
Assignment to another Eos_bf_tabul.
double ent2_min
Lower boundary of the log-enthalpy interval (fluid 2 = p)
virtual Eos * trans2Eos() const
Makes a translation from Eos_bifluid to Eos .
string tablename
Name of the file containing the tabulated data (be careful, Eos_bifluid uses char*) ...
Eos_bf_tabul(const char *name_i, const char *table, const char *path, double mass1, double mass2)
Standard constructor.
virtual bool nbar_ent_p(const double ent1, const double ent2, const double delta2, double &nbar1, double &nbar2) const
Computes both baryon densities from the log-enthalpies.