31 #include <boost/numeric/ublas/vector.hpp> 32 #include <o2scl/eos_had_apr.h> 33 #include <o2scl/table.h> 34 #include <o2scl/constants.h> 35 #include <o2scl/fermion_nonrel.h> 36 #include <o2scl/tov_solve.h> 37 #include <o2scl/deriv_cern.h> 38 #include <o2scl/mroot_cern.h> 39 #include <o2scl/test_mgr.h> 40 #include <o2scl/hdf_io.h> 43 using namespace o2scl;
127 static const int low_phase=1;
128 static const int mixed_phase=2;
129 static const int high_phase=3;
148 if (x[0]<0.0 || x[1]<0.0 || x[2]<0.0 || x[3]<0.0)
return 1;
151 ap.
pion=eos_had_apr::ldp;
163 ap.
pion=eos_had_apr::hdp;
181 y[3]=p2.
n-e2.
n-mu2.
n;
194 ap.
pion=eos_had_apr::ldp;
205 ap.
pion=eos_had_apr::hdp;
216 y[0]=nb-chi*(n.
n+p.
n)-(1.0-chi)*(n2.
n+p2.
n);
220 y[4]=p2.
n-e2.
n-mu2.
n;
235 ap.
pion=eos_had_apr::ldp;
238 ap.
pion=eos_had_apr::hdp;
252 if (n.
n<0.0 || p.
n<0.0) {
256 ap.
pion=eos_had_apr::ldp;
281 if (n2.
n<0.0 || p2.
n<0.0)
return 1;
283 ap.
pion=eos_had_apr::hdp;
312 if (phase==low_phase) chi=1.0;
313 else if (phase==high_phase) chi=0.0;
316 if (n.
n<0.0 || n2.
n<0.0)
return 1;
317 if (p.
n<0.0 || p2.
n<0.0)
return 1;
319 ap.
pion=eos_had_apr::ldp;
322 ap.
pion=eos_had_apr::hdp;
331 y[0]=nb-chi*(n.
n+p.
n)-(1.0-chi)*(n2.
n+p2.
n);
332 y[1]=chi*p.
n+(1.0-chi)*p2.
n-e.
n-mu.
n;
337 if (phase==mixed_phase) y[5]=hb.
pr-hb2.
pr;
339 if (phase==low_phase) {
342 }
else if (phase==mixed_phase) {
344 tot.
ed=hb.
ed*chi+hb2.
ed*(1.0-chi)+l.
ed;
356 std::vector<double> line;
357 line.push_back(tot.
ed);
358 line.push_back(tot.
pr);
360 line.push_back((chi*n.
n+(1.0-chi)*n2.
n));
361 line.push_back((chi*p.
n+(1.0-chi)*p2.
n));
364 line.push_back(n2.
n);
365 line.push_back(p2.
n);
368 line.push_back(mu.
n);
369 line.push_back(n.
mu);
370 line.push_back(p.
mu);
371 line.push_back(e.
mu);
372 line.push_back(mu.
mu);
373 line.push_back(n.
ms);
374 line.push_back(p.
ms);
375 line.push_back(n2.
ms);
376 line.push_back(p2.
ms);
377 line.push_back(n.
kf);
378 line.push_back(p.
kf);
379 line.push_back(n2.
kf);
380 line.push_back(p2.
kf);
381 line.push_back(e.
kf);
382 line.push_back(mu.
kf);
397 if (phase==low_phase) chi=1.0;
398 else if (phase==high_phase) chi=0.0;
401 if (n.
n<0.0 || n2.
n<0.0)
return 1;
406 ap.
pion=eos_had_apr::ldp;
409 ap.
pion=eos_had_apr::hdp;
413 y[1]=nb-chi*(n.
n+p.
n)-(1.0-chi)*(n2.
n+p2.
n);
415 if (phase==mixed_phase) y[2]=hb.
pr-hb2.
pr;
417 if (phase==low_phase) {
420 }
else if (phase==mixed_phase) {
422 tot.
ed=hb.
ed*chi+hb2.
ed*(1.0-chi);
438 if (phase==low_phase) chi=1.0;
439 else if (phase==high_phase) chi=0.0;
442 if (n.
n<0.0 || n2.
n<0.0)
return 1;
447 ap.
pion=eos_had_apr::ldp;
450 ap.
pion=eos_had_apr::hdp;
454 y[1]=nb-chi*(n.
n+p.
n)-(1.0-chi)*(n2.
n+p2.
n);
456 if (phase==mixed_phase) y[2]=hb.
pr-hb2.
pr;
458 if (phase==low_phase) {
461 }
else if (phase==mixed_phase) {
463 tot.
ed=hb.
ed*chi+hb2.
ed*(1.0-chi);
475 double nn1,np1,mun1,mup1;
476 double nn2,np2,mun2,mup2;
482 ap.
pion=eos_had_apr::ldp;
491 ap.
pion=eos_had_apr::ldp;
506 ey[0]=(mun1-mun2)/1.0;
507 ey[1]=(th1.
pr-th2.
pr)/1.0;
509 if (nv>3)ey[3]=
barn-u*(np1+nn1)-(1.0-u)*nn2;
512 hb.
ed=u*th1.
ed+(1.0-u)*th2.
ed;
513 tot.
ed=u*th1.
ed+(1.0-u)*th2.
ed+e.
ed;
517 ap.
pion=eos_had_apr::ldp;
526 double nn1,np1,mun1,mup1;
527 double nn2,np2,mun2,mup2;
531 if (ex[0]<0.0 || ex[1]<0.0 || ex[2]<0.0 || ex[3]<0.0) {
537 ap.
pion=eos_had_apr::ldp;
546 ap.
pion=eos_had_apr::ldp;
561 ey[0]=(mun1-mun2)/1.0;
562 ey[1]=(th1.
pr-th2.
pr)/1.0;
563 ey[2]=e.
n-u*np1-(1.0-u)*np2;
564 ey[3]=(mup1-mup2)/1.0;
565 if (nv>4) ey[4]=
barn-u*(np1+nn1)-(1.0-u)*(nn2+np2);
568 hb.
ed=u*th1.
ed+(1.0-u)*th2.
ed;
569 tot.
ed=u*th1.
ed+(1.0-u)*th2.
ed+e.
ed;
573 ap.
pion=eos_had_apr::ldp;
614 std::string prefix=
"ex_eos_had_apr_";
615 cout <<
"Set output prefix to '" << prefix <<
"' ." << endl;
619 double nbstart, nb_end, dnb;
640 std::placeholders::_1,std::placeholders::_2,std::placeholders::_3);
644 std::placeholders::_1,std::placeholders::_2,std::placeholders::_3);
648 std::placeholders::_1,std::placeholders::_2,std::placeholders::_3);
652 std::placeholders::_1,std::placeholders::_2,std::placeholders::_3);
656 std::placeholders::_1,std::placeholders::_2,std::placeholders::_3);
660 std::placeholders::_1,std::placeholders::_2,std::placeholders::_3);
664 std::placeholders::_1,std::placeholders::_2,std::placeholders::_3);
668 std::placeholders::_1,std::placeholders::_2,std::placeholders::_3);
682 solver_trans_density.
tol_abs=1.0e-10;
683 solver_trans_density.
tol_rel=1.0e-12;
687 "nn np nn1 np1 nn2 np2 chi ne nmu "+
689 "msn msp msn2 msp2 kfn kfp kfn2 kfp2 kfe kfmu");
724 cout <<
"--------- Saturation properties --------------------\n" << endl;
746 cout <<
"--------- Nuclear matter ---------------------------\n" << endl;
752 ap.
pion=eos_had_apr::ldp;
755 for(nb=nbstart;nb<nb_switch;nb+=dnb) {
776 for(nb=nb_switch;nb<=nb_end;nb+=dnb) {
778 if (phase!=mixed_phase) ret=solver.
msolve(2,x,f_nucmixed);
779 else ret=solver.
msolve(3,x,f_nucmixed);
782 O2SCL_ERR(
"Solving nuclear matter failed.",
786 if (hb.
pr<hb2.
pr && phase==low_phase) {
787 cout <<
"Mixed phase begins near nb=" << nb <<
" fm^{-3}." << endl;
792 }
else if (phase==mixed_phase && x[2]<0.0) {
793 cout <<
"Mixed phase ends near nb=" << nb <<
" fm^{-3}." << endl;
801 string fn1=((string)prefix)+
"nuc.o2";
805 cout <<
"Generated file '" << fn1 <<
"'." << endl;
812 cout <<
"--------- Neutron matter ---------------------------\n" << endl;
814 ap.
pion=eos_had_apr::ldp;
822 for(nb=nbstart;nb<nb_switch;nb+=dnb) {
844 for(nb=nb_switch;nb<=nb_end;nb+=dnb) {
846 if (phase!=mixed_phase) {
847 ret=solver.
msolve(2,x,f_neutmixed);
849 ret=solver.
msolve(3,x,f_neutmixed);
853 O2SCL_ERR(
"Solving neutron matter failed.",
857 if (hb.
pr<hb2.
pr && phase==low_phase) {
858 cout <<
"Mixed phase begins near nb=" << nb <<
" fm^{-3}." << endl;
863 }
else if (phase==mixed_phase && x[2]<0.0) {
864 cout <<
"Mixed phase ends near nb=" << nb <<
" fm^{-3}." << endl;
872 string fn2=((string)prefix)+
"neut.o2";
876 cout <<
"Generated file '" << fn2 <<
"'." << endl;
883 cout <<
"--------- Neutron star matter ----------------------\n" << endl;
888 cout <<
"Maxwell construction." << endl;
893 filenm+=
"fig7.1.txt";
894 fout.open(filenm.c_str());
895 fout.setf(ios::scientific);
897 for(f7x=0.5;f7x>=-0.001;f7x-=0.025) {
898 ret=solver.
msolve(1,x,f_fig7fun);
902 fout << f7x <<
" " << x[0] << endl;
905 cout <<
"Generated file '" << filenm <<
"'." << endl;
914 solver.
msolve(4,x,f_maxwell_fig7);
916 double nb_low=x[0]+x[1], nb_high=x[2]+x[3];
917 cout <<
"Mixed phase begins at nb=" << nb_low <<
" fm^{-3}." << endl;
918 cout <<
"Mixed phase ends at nb=" << nb_high <<
" fm^{-3}." << endl;
922 filenm+=
"fig7.2.txt";
923 fout.open(filenm.c_str());
927 for(nb=0.02;nb<nb_low*1.00001;nb+=(nb_low-nbstart)/20.0) {
928 ret=solver.
msolve(1,x,f_nstar_low);
932 cout << x[0] <<
" " << y[0] << endl;
933 O2SCL_ERR(
"Solving Maxwell construction failed.",
936 fout << p.
n/nb <<
" " << nb << endl;
941 x[0]=(1.0-0.07)*nb_low;
943 x[2]=(1.0-0.06)*nb_high;
946 dnb=(nb_high-nb_low)/20.0;
947 for(nb=nb_low+dnb;nb<=nb_high*1.00001;nb+=dnb) {
948 ret=solver.
msolve(5,x,f_mixedmaxwell);
951 O2SCL_ERR(
"Solving Maxwell construction (part 2) failed.",
954 fout << (chi*p.
n+(1.0-chi)*p2.
n)/nb <<
" " << nb << endl;
960 for(nb=nb_high;nb<nb_end;nb+=(nb_end-nb_high)/40.0) {
961 ret=solver.
msolve(1,x,f_nstar_high);
964 O2SCL_ERR(
"Solving Maxwell construction (part 3) failed.",
967 fout << p2.
n/nb <<
" " << nb << endl;
970 cout <<
"Generated file '" << filenm <<
"'." << endl;
976 cout <<
"\nGibbs construction." << endl;
979 ap.
pion=eos_had_apr::ldp;
984 for(nb=nbstart;nb<=0.1701;nb+=dnb) {
985 ret=solver.
msolve(1,x,f_nstar_low);
988 O2SCL_ERR(
"Solving Gibbs construction failed.",
1002 for(nb=0.17;nb<=nb_end;nb+=dnb) {
1004 if (phase!=mixed_phase) {
1005 ret=solver.
msolve(5,x,f_nstar_mixed);
1008 ret=solver.
msolve(6,x,f_nstar_mixed);
1014 if (hb.
pr<hb2.
pr && phase==low_phase) {
1015 cout <<
"Mixed phase begins at nb=" << nb <<
" fm^{-3}." << endl;
1019 }
else if (phase==mixed_phase && x[5]<0.0) {
1020 cout <<
"Mixed phase ends at nb=" << nb <<
" fm^{-3}." << endl;
1031 cout <<
"\nEstimate of transition density." << endl;
1037 std::placeholders::_1,std::placeholders::_2,std::placeholders::_3);
1041 std::placeholders::_1,std::placeholders::_2,std::placeholders::_3);
1043 solver_trans_density.
tol_abs/=1.0e4;
1044 solver_trans_density.
tol_rel/=1.0e4;
1049 ret=solver_trans_density.
msolve(3,newx,nuclei_f);
1051 nucleimat(3,newx,newy);
1052 if (fabs(newy[0])>1.0e-8 || fabs(newy[1])>1.0e-8 ||
1053 fabs(newy[2])>1.0e-8) {
1055 cout <<
"Problem in transition density (1)." << endl;
1057 cout << newx[0] <<
" " << newy[0] << endl;
1058 cout << newx[1] <<
" " << newy[1] << endl;
1059 cout << newx[2] <<
" " << newy[2] << endl;
1063 cout <<
"Without proton drip: density: " << newx[0]+newx[1]
1064 <<
" fm^{-3},\n pressure: " 1065 << at.
interp(
"nb",newx[0]+newx[1],
"pr") <<
" fm^{-4}." << endl;
1067 solver_trans_density.
tol_abs=5.0e-8;
1068 solver_trans_density.
tol_rel=5.0e-8;
1071 ret=solver_trans_density.
msolve(4,newx,nucleip_f);
1073 nucleimat_pdrip(4,newx,newy);
1074 if (fabs(newy[0])>1.0e-8 || fabs(newy[1])>1.0e-8 ||
1075 fabs(newy[2])>1.0e-8 || fabs(newy[3])>1.0e-8) {
1077 cout <<
"Problem in transition density (2)." << endl;
1079 cout << newx[0] <<
" " << newy[0] << endl;
1080 cout << newx[1] <<
" " << newy[1] << endl;
1081 cout << newx[2] <<
" " << newy[2] << endl;
1082 cout << newx[3] <<
" " << newy[3] << endl;
1086 nucleimat_pdrip(4,newx,newy);
1087 cout << newx[0] <<
" " << newy[0] << endl;
1088 cout << newx[1] <<
" " << newy[1] << endl;
1089 cout << newx[2] <<
" " << newy[2] << endl;
1090 cout << newx[3] <<
" " << newy[3] << endl;
1091 cout <<
"With proton drip: density: " << newx[0]+newx[1]
1092 <<
" fm^{-3},\n pressure: " 1093 << at.
interp(
"nb",newx[0]+newx[1],
"pr") <<
" fm^{-4}." << endl;
1096 solver_trans_density.
tol_abs=1.0e-16;
1097 solver_trans_density.
tol_rel=1.0e-16;
1101 string fn3=((string)prefix)+
"nstar.o2";
1105 cout <<
"Generated file '" << fn3 <<
"'." << endl;
1112 cout <<
"--------- TOV solver, M vs. R. ---------------------\n" << endl;
1118 atov.
set_units(
"1/fm^4",
"1/fm^4",
"1/fm^3");
1127 std::shared_ptr<table_units<> > tov_tmp=atov.
get_results();
1128 cout <<
"Maximum mass is " << tov_tmp->max(
"gm")
1129 <<
" solar masses." << endl;
1130 t.
test_rel(tov_tmp->max(
"gm"),2.191338,5.0e-4,
"max mass.");
1134 string fn4=((string)prefix)+
"mvsr.o2";
1138 cout <<
"Generated file '" << fn4 <<
"'." << endl;
1141 cout <<
"--------- TOV solver, 1.4 Msun ---------------------\n" << endl;
1144 std::shared_ptr<table_units<> > tov_tmp2=atov.
get_results();
1145 cout <<
"Aprpoximate radius of a 1.4 solar mass neutron star is " 1146 << tov_tmp2->get(
"r",tov_tmp2->lookup(
"gm",1.4)) <<
" km." << endl;
1147 t.
test_rel(tov_tmp2->get(
"r",tov_tmp2->lookup(
"gm",1.4)),11.46630,
1150 string fn5=((string)prefix)+
"m14.o2";
1154 cout <<
"Generated file '" << fn5 <<
"'." << endl;
1166 cout.setf(ios::scientific);
void set_eos(eos_tov &ter)
Set the EOS to use.
virtual int mvsr()
Calculate the mass vs. radius curve.
thermo hb
Baryon thermodynamics for low-density phase.
double chi
Volume fraction of low-density phase.
int nstar_mixed(size_t nv, const ubvector &x, ubvector &y)
Solve for neutron star matter (mixed phase)
int fig7fun(size_t nv, const ubvector &x, ubvector &y)
Function to construct Fig. 7.
double comp
Compression modulus in .
bool calc_gpot
calculate the gravitational potential (default false)
virtual int fixed(double target_mass, double pmax=1.0e20)
Calculate the profile of a star with fixed mass.
void default_low_dens_eos()
Standard crust EOS from Negele73 and Baym71tg.
fermion e2
Electron for high-density phase.
void run()
Main driver, computing the APR EOS and the associated M vs. R curve.
int choice
Choice of model from APR.
std::function< int(size_t, const boost::numeric::ublas::vector< double > &, boost::numeric::ublas::vector< double > &) > mm_funct
fermion p2
Proton for high-density phase.
std::shared_ptr< table_units<> > get_results()
Return the results data table.
thermo l2
Leptonic thermodynamics for high-density phase.
virtual void calc_mu_zerot(fermion &f)
virtual double convert(std::string from, std::string to, double val)
fermion mu
Muon for low-density phase.
fermion p
Proton for low-density phase.
int verbose
control for output (default 1)
lib_settings_class o2scl_settings
int pion
Choice of phase (default best)
thermo tot
Total thermodynamics.
void store_data()
Write a line of data to the table.
double muq
Charge chemical potential.
int nstar_low(size_t nv, const ubvector &x, ubvector &y)
Solve for neutron star matter (low-density phase)
fermion mu2
Muon for high-density phase.
void set_units(double s_efactor=1.0, double s_pfactor=1.0, double s_nbfactor=1.0)
Set output units for the table.
eos_had_apr ap
Base APR EOS.
bool test_rel(data_t result, data_t expected, data_t rel_error, std::string description)
int nucmixed(size_t nv, const ubvector &x, ubvector &y)
Solve for nuclear matter (mixed phase)
double esym
Symmetry energy in .
void read_table(table_units<> &eosat, std::string s_cole, std::string s_colp, std::string s_colnb="")
Specify the EOS through a table.
fermion n2
Neutron for high-density phase.
table_units at
Table for output.
fermion e
Electron for low-density phase.
double n0
Saturation density in .
virtual int saturation()
Calculates some of the EOS properties at the saturation density.
fermion n
Neutron for low-density phase.
double f7x
Proton fraction for Fig. 7.
virtual int calc_e(fermion &n, fermion &p, thermo &th)
Equation of state as a function of density.
mroot_hybrids solver
General solver.
int nucleimat_pdrip(size_t nv, const ubvector &ex, ubvector &ey)
Solve for phase transition to nuclei with a proton drip.
void set_unit(std::string scol, std::string unit)
int nucleimat(size_t nv, const ubvector &ex, ubvector &ey)
Solve for phase transition to nuclei.
double interp(std::string sx, double x0, std::string sy)
void open_or_create(std::string fname)
virtual const char * get_str()=0
int mixedmaxwell(size_t nv, const ubvector &x, ubvector &y)
Maxwell construction of the nuclear matter mixed phase.
void line_of_names(std::string newheads)
An EOS for the TOV solver using simple linear interpolation and an optional crust EOS...
void select(int model_index=1)
Select model.
Additional functions to read and write EOS data to HDF5 files.
convert_units & get_convert_units()
Compute the APR EOS with a Gibbs construction and the mass versus radius curve [Example class]...
int maxwell_fig7(size_t nv, const ubvector &x, ubvector &y)
Function for the Maxwell construction in Fig. 7.
void line_of_data(size_t nv, const vec2_t &v)
thermo hb2
Baryon thermodynamics for high-density phase.
virtual void init(double mass, double dof)
Class to solve the Tolman-Oppenheimer-Volkov equations.
size_t get_ncolumns() const
thermo l
Leptonic thermodynamics for low-density phase.
int nstar_high(size_t nv, const ubvector &x, ubvector &y)
Solve for neutron star matter (high-density phase)
void hdf_output(hdf_file &hf, o2scl::uniform_grid< double > &h, std::string name)
double mub
Baryon chemical potential.
int neutmixed(size_t nv, const ubvector &x, ubvector &y)
Solve for neutron matter (mixed phase)
mroot_hybrids solver_trans_density
Solver for transition densities (lower tolerances)
hdf_file hf
HDF file for output.
virtual int msolve(size_t nn, vec_t &xx, func_t &ufunc)
EOS from Akmal, Pandharipande, and Ravenhall.
void set_output_level(int l)
virtual void set_thermo(thermo &th)
Set class thermo object.
int verbose
Control for output (default 1)
double msom
Effective mass (neutron)
const double mass_electron
fermion_zerot fzt
Compute zero-temperature thermodynamics.
double eoa
Binding energy (without the rest mass) in .
virtual void set_mroot(mroot<> &mr)
Set class mroot object for use in calculating chemical potentials from densities. ...
virtual void set_n_and_p(fermion &n, fermion &p)
Set neutron and proton.
void inc_maxlines(size_t llines)
deriv_cern cd
Derivative object.