29 #include <o2scl/test_mgr.h>
30 #include <o2scl/eos_had_skyrme.h>
31 #include <o2scl/fermion_nonrel.h>
32 #include <o2scl/nstar_cold.h>
33 #include <o2scl/format_float.h>
34 #include <o2scl/hdf_file.h>
35 #include <o2scl/hdf_io.h>
36 #include <o2scl/hdf_eos_io.h>
37 #include <o2scl/cli.h>
38 #include <o2scl/lib_settings.h>
41 using namespace o2scl;
185 n.non_interacting=
false;
186 p.non_interacting=
false;
197 file_prefix=
"skyrme_";
210 tneut.new_column(
"nb");
211 tneut.new_column(
"pr_neut");
212 tneut.new_column(
"pr_nstar");
214 tneut2.new_column(
"ed");
215 tneut2.new_column(
"pr_neut");
216 tneut2.new_column(
"pr_nstar");
224 tneut2.
set_unit(
"pr_nstar",
"fm^-4");
226 for(
double nbt=0.06;nbt<0.9001;nbt+=0.02) {
229 int ret=sk.
calc_e(n,p,th);
230 double line[3]={nbt,th.
pr,0.0};
231 tneut.line_of_data(3,line);
232 double line2[3]={th.
ed,th.
pr,0.0};
233 tneut2.line_of_data(3,line2);
241 bool failed_once=
false;
247 for(
double nbt=0.005;nbt<0.1601 && failed_once==
false;nbt+=0.01) {
250 int ret=sk.
calc_e(n,p,th);
257 if (th.
ed-n.n*n.m<0.0) {
263 if (failed_once==
false) {
264 cout <<
"Pressure increases and energy of neutrons is positive.\n"
267 cout <<
"Pressure decreased or energy of neutrons is negative.\n"
277 double g_kf[4]={0.06,0.14,0.26,0.54};
278 double g_rat[4]={0.80,0.67,0.57,0.52};
281 ln.line_of_names(
"nb epb kf fermi rat");
282 if (verbose>0) cout <<
"Low density neutron matter: " << endl;
283 for(
double nbt=0.000003;nbt<0.0051;nbt*=1.2) {
286 int ret=sk.
calc_e(n,p,th);
295 cout << nbt <<
" " << epb <<
" "
296 << n.kf <<
" " << fermi <<
" " << epb/fermi << endl;
298 double line[5]={nbt,epb,n.kf,fermi,epb/fermi};
299 ln.line_of_data(5,line);
301 if (verbose>1) cout << endl;
312 for(
size_t i=0;i<4;i++) {
314 cout << g_kf[i] <<
" " << g_rat[i] <<
" "
315 << ln.interp(
"kf",g_kf[i],
"rat") <<
" "
316 << ln.interp(
"kf",g_kf[i],
"epb") << endl;
318 res.
neut_qual+=g_kf[i]*fabs(g_rat[i]-ln.interp(
"kf",g_kf[i],
"rat"));
320 ln.add_constant(
"neut_qual",res.
neut_qual);
322 cout <<
"Quality: " << res.
neut_qual << endl;
345 cout <<
"{\"n0\":" << sk.
n0 <<
"," << endl;
348 cout <<
"\"msom\":" << sk.
msom <<
"," << endl;
350 cout <<
"\"S2\":" << res.
alt_S <<
"," << endl;
351 cout <<
"\"L\":" << res.
L <<
"}" << endl;
352 }
else if (verbose>0) {
353 cout <<
"Saturation: " << endl;
354 cout <<
"n_0=" << sk.
n0 <<
" fm^-3" << endl;
357 cout <<
"M^*/M=" << sk.
msom << endl;
359 cout <<
"S2=" << res.
alt_S <<
" MeV" << endl;
360 cout <<
"L=" << res.
L <<
" MeV" << endl;
368 if (fabs(sk.
n0-0.16)>0.013) {
381 cout <<
"Bad saturation." << endl;
393 cout <<
"EOS:" << endl;
396 for(
double nb=0.16;nb<=2.0001;nb+=0.001) {
402 if (n.mu-p.mu-me<0.0) {
408 cout <<
"Pure neutron matter after nb=" << nb_last << endl;
418 cout <<
"Neutron stars:" << endl;
425 string fn=file_prefix+res.
name+
"_eos.o2";
429 fn=file_prefix+res.
name+
"_mvsr.o2";
436 cout <<
"M_{max} = " << tr->max(
"gm") <<
" R_{max} = "
437 << tr->get(
"r",tr->lookup(
"gm",tr->max(
"gm")))
438 <<
" cent. density = "
439 << tr->get(
"nb",tr->lookup(
"gm",tr->max(
"gm"))) << endl;
441 << tr->get(
"r",tr->lookup(
"gm",1.4)) <<
" cent. density = "
442 << tr->get(
"nb",tr->lookup(
"gm",1.4)) << endl;
445 res.
m_max=tr->max(
"gm");
446 res.
r_max=tr->get(
"r",tr->lookup(
"gm",tr->max(
"gm")));
447 res.
nb_max=tr->get(
"nb",tr->lookup(
"gm",tr->max(
"gm")));
448 res.
r_14=tr->get(
"r",tr->lookup(
"gm",1.4));
449 res.
nb_14=tr->get(
"nb",tr->lookup(
"gm",1.4));
453 if (nbtop>2.0) nbtop=2.0;
454 if (nbtop<0.4) nbtop=0.7;
461 for(
double nb=0.1;nb<=nbtop;nb+=0.01) {
462 if (te->get(
"np",te->lookup(
"nb",nb))<1.0e-5) {
475 if (tr->max(
"gm")<1.6) {
484 nst.
pressure_dec<tr->get(
"nb",tr->lookup(
"gm",tr->max(
"gm")))) {
485 cout <<
"Pressure decreases in maximum mass star" << endl;
490 nst.
acausal<tr->get(
"nb",tr->lookup(
"gm",tr->max(
"gm")))) {
491 cout <<
"Acausal before central density of maximum mass star." << endl;
492 cout <<
"acausal: " << nst.
acausal << endl;
504 int json(vector<string> &sv,
bool itive_com) {
506 cout <<
"No model to summarize." << endl;
510 cout <<
"Model: " << sv[1] << endl;
511 cout <<
"-----------------------------------------------------"
512 <<
"-----------------------" << endl;
536 int summary(vector<string> &sv,
bool itive_com) {
538 cout <<
"No model to summarize." << endl;
542 cout <<
"Model: " << sv[1] << endl;
543 cout <<
"-----------------------------------------------------"
544 <<
"-----------------------" << endl;
564 compare_neut_nstar();
568 cout <<
"-----------------------------------------------------"
569 <<
"-----------------------" << endl;
575 int test(vector<string> &sv,
bool itive_com) {
579 bool of_old=output_files;
581 file_prefix=
"ex_eos_had_skyrme_";
585 args.push_back(
"summary");
586 args.push_back(
"SLy4");
602 int store(vector<string> &sv,
bool itive_com) {
605 cout <<
"No filename specified in 'store'." << endl;
614 cout <<
"Wrote model '" << name <<
"' to file named '"
615 << sv[1] <<
"'." << endl;
622 int unedf(vector<string> &sv,
bool itive_com) {
631 double coups[13][3]={
632 {-706.382928878428856,-779.3730087208653,-650.796319465688839},
633 {240.049520427681131,287.722131583286796,291.664014339185485},
634 {868.871771539645351,891.47789044234969,768.32770588203357},
635 {-69.0518957481631617,-200.587774317884879,-283.187292227492492},
636 {-12.9172408208016112,-0.989915057807676746,9.78520558824309639},
637 {-45.1894169426759476,-33.6320970701835549,-23.3573612977035339},
638 {0.321955989588264435,0.2700180115027076,0.351455132555483607},
639 {-55.2606,-45.135131022237303,-46.8314091470605973},
640 {-55.6226,-145.382167908057,-113.163790795259004},
641 {-79.5308,-74.0263331764599,-64.3088624157838069},
642 {45.6302,-35.6582611147917,-38.6501946851355029},
643 {0.0,0.0,-54.4333635973721002},
644 {0.0,0.0,-65.9030310445938028}
650 for(
size_t i=0;i<3;i++) {
651 cout <<
"UNEDF" << i <<
":" << endl;
654 for(
size_t j=0;j<13;j++) {
660 sk.
alt_params_set(coups[0][i],coups[1][i],coups[2][i],coups[3][i],
661 coups[4][i],coups[5][i],coups[7][i],coups[8][i],
662 coups[9][i],coups[10][i],coups[6][i]);
670 t.
test_rel(sk.x0,0.00974375,1.0e-4,
"");
671 t.
test_rel(sk.x1,-1.77784395,1.0e-4,
"");
672 t.
test_rel(sk.x2,-1.67699035,1.0e-4,
"");
673 t.
test_rel(sk.x3,-0.38079041,1.0e-4,
"");
681 t.
test_rel(sk.x0,0.05375692,1.0e-4,
"");
682 t.
test_rel(sk.x1,-5.07723238,1.0e-4,
"");
683 t.
test_rel(sk.x2,-1.36650561,1.0e-4,
"");
684 t.
test_rel(sk.x3,-0.16249117,1.0e-4,
"");
693 coups[8][i],coups[9][i],coups[10][i]);
696 cout <<
"n0: " << sk.
n0 << endl;
705 cout <<
"alpha: " << sk.alpha << endl;
734 sk.
reference=((string)
"M. Kortelainen, T. Lesinski, ")+
735 "J. More, W. Nazarewicz, J. Sarich, N. Schunck, "+
736 "M. V. Stoitsov, and S. Wild, Phys. Rev. C 82, "+
739 sk.
reference=((string)
"M. Kortelainen, J. McDonnell, ")+
740 "W. Nazarewicz, P.-G. Reinhard, J. Sarich, N. Schunck, "+
741 "M. V. Stoitsov, and S. M. Wild, Phys. Rev. C 85, "+
744 sk.
reference=((string)
"M. Kortelainen, J. McDonnell, ")+
745 "W. Nazarewicz, E. Olsen, P.-G. Reinhard, J. Sarich, "+
746 "N. Schunck, S. M. Wild, D. Davesne, J. Erler, "+
747 "and A. Pastore, Phys. Rev. C 89, 054314 (2014)";
762 int load(vector<string> &sv,
bool itive_com) {
765 cout <<
"No model specified in 'load'." << endl;
773 cout <<
"Loaded model '" << name <<
"'." << endl;
780 int run_all(vector<string> &sv,
bool itive_com) {
783 std::string mlist[200], stemp;
786 ifstream fin(fname.c_str());
788 for(
size_t i=0;i<nmods;i++) {
793 ofstream fouu(
"table.csv");
794 fouu <<
"Name, n0, B, K, ";
795 fouu <<
"S, L, Mmax, ";
796 fouu <<
"Rmax, nB_cent_max, R1.4, ";
797 fouu <<
"nB_cent_14, acausal, pressure_dec, ";
798 fouu <<
"neut_qual, max_mass_ok, ";
799 fouu <<
"inc_pressure, pos_neut, good_sat, ";
800 fouu <<
"pure_neut, other, success" << endl;
802 ofstream fout(
"table.html");
803 fout <<
"<html><body>" << endl;
804 fout <<
"<table border=0 cellspacing=0><tr bgcolor=\"#bbbbbb\">" << endl;
805 fout <<
"<td>Name </td>" << endl;
806 fout <<
"<td>n<sub>0</sub> (fm<sup>-3</sup>)"
807 <<
" </td>" << endl;
808 fout <<
"<td>B (MeV) </td>" << endl;
809 fout <<
"<td>K (MeV) </td>" << endl;
810 fout <<
"<td>S (MeV) </td>" << endl;
811 fout <<
"<td>L (MeV) </td>" << endl;
812 fout <<
"<td>M<sub>max</sub> (M<sub>sun</sub>)"
813 <<
" </td>" << endl;
814 fout <<
"<td>R<sub>max</sub> </td>" << endl;
815 fout <<
"<td>n<sub>B,cent,max</sub> </td>" << endl;
816 fout <<
"<td>R<sub>1.4</sub> </td>" << endl;
817 fout <<
"<td>n<sub>B,cent,1.4</sub> </td>" << endl;
818 fout <<
"<td>acausal </td>" << endl;
819 fout <<
"<td>pressure_dec </td>" << endl;
820 fout <<
"<td>neut_qual </td>" << endl;
821 fout <<
"<td>max_mass_ok </td>" << endl;
822 fout <<
"<td>inc_pressure </td>" << endl;
823 fout <<
"<td>pos_neut </td>" << endl;
824 fout <<
"<td>good_sat </td>" << endl;
825 fout <<
"<td>pure_neut </td>" << endl;
826 fout <<
"<td>other </td>" << endl;
827 fout <<
"<td>success </td>" << endl;
828 fout <<
"</tr>" << endl;
833 static const size_t N=66;
834 int list[N]={0,100,101,102,103,109,110,113,114,115,121,122,123,124,
835 125,126,127,128,129,130,131,132,133,134,135,145,147,17,1,
836 25,26,27,28,29,3,40,41,42,43,44,4,51,53,54,5,
837 62,63,64,64,66,67,68,69,6,71,73,75,76,77,7,81,82,84,
840 for(
size_t j=0;j<N;j++) {
846 cout <<
"Running model: " << i << endl;
849 tmp.push_back(mlist[i]);
867 ofstream fx(
"jim.dat",ios::app);
868 fx.setf(ios::scientific);
869 fx.setf(ios::showpos);
871 fx << mlist[i] <<
" ";
872 fx << res.
S <<
" " << res.
alt_S <<
" "
875 fx << sk.t0 <<
" " << sk.t1 <<
" " << sk.t2 <<
" " << sk.t3 <<
" ";
876 fx << sk.x0 <<
" " << sk.x1 <<
" " << sk.x2 <<
" " << sk.x3 <<
" ";
897 cout << sk.t0 <<
" " << sk.t1 <<
" " << sk.t2 <<
" " << sk.t3 << endl;
898 cout << sk.x0 <<
" " << sk.x1 <<
" " << sk.x2 <<
" " << sk.x3 << endl;
901 fout <<
"<tr bgcolor=\"#dddddd\">";
907 fout <<
"<td><a href=\"http://o2scl.svn.sourceforge.net/viewvc/"
908 <<
"o2scl/trunk/data/o2scl/skdata/" << res.
name
909 <<
"\">" << res.
name <<
"</a></td>";
910 fout <<
"<td>" << fd.
convert(res.
n0) <<
"</td>";
911 fout <<
"<td>" << fd.
convert(res.
B) <<
"</td>";
912 fout <<
"<td>" << fd.
convert(res.
K) <<
"</td>";
913 fout <<
"<td>" << fd.
convert(res.
S) <<
"</td>";
914 fout <<
"<td>" << fd.
convert(res.
L) <<
"</td>";
924 else fout <<
"<td>False</td>";
926 else fout <<
"<td>False</td>";
927 if (res.
pos_neut) fout <<
"<td>True</td>";
928 else fout <<
"<td>False</td>";
929 if (res.
good_sat) fout <<
"<td>True</td>";
930 else fout <<
"<td>False</td>";
931 if (res.
pure_neut) fout <<
"<td>True</td>";
932 else fout <<
"<td>False</td>";
933 fout <<
"<td>" << res.
other <<
"</td>";
934 if (res.
success) fout <<
"<td>True</td>";
935 else fout <<
"<td>False</td>";
936 fout <<
"</tr>" << endl;
939 fouu << res.
name <<
", ";
940 fouu << res.
n0 <<
", ";
941 fouu << res.
B <<
", ";
942 fouu << res.
K <<
", ";
943 fouu << res.
S <<
", ";
944 fouu << res.
L <<
", ";
945 fouu << res.
m_max <<
", ";
946 fouu << res.
r_max <<
", ";
947 fouu << res.
nb_max <<
", ";
948 fouu << res.
r_14 <<
", ";
949 fouu << res.
nb_14 <<
", ";
954 else fouu <<
"False, ";
956 else fouu <<
"False, ";
958 else fouu <<
"False, ";
960 else fouu <<
"False, ";
962 else fouu <<
"False, ";
963 fouu << res.
other <<
", ";
964 if (res.
success) fouu <<
"True ";
965 else fouu <<
"False ";
971 fout <<
"</table></body></html>" << endl;
983 int main(
int argc,
char *argv[]) {
985 cout.setf(ios::scientific);
993 cl.
prompt=
"ex_eos_had_skyrme> ";
995 int comm_option_cl_param=1;
996 int comm_option_both=2;
998 static const int narr=7;
1000 {0,
"run-all",
"Run all internally stored Skyrme models.",0,0,
"",
"",
1003 {
's',
"store",
"Store current model.",1,1,
"",
"",
1006 {
'l',
"load",
"Load internally stored model.",1,1,
"",
"",
1009 {0,
"unedf",
"Desc.",0,0,
"",
"",
1012 {
't',
"test",
"Test ex_eos_had_skyrme.",0,0,
"",
"",
1015 {
'u',
"summary",
"Summarize the properties of a Skyrme model.",
1019 {
'j',
"json",
"Summarize the properties of a Skyrme model.",
1033 p_verbose.
help=
"Verbose (default 1).";
1034 cl.
par_list.insert(make_pair(
"verbose",&p_verbose));
1038 p_output_files.
help=
"Output files (default 0).";
1039 cl.
par_list.insert(make_pair(
"output-files",&p_output_files));
1043 p_file_prefix.
help=
"File prefix (default \"\").";
1044 cl.
par_list.insert(make_pair(
"file-prefix",&p_file_prefix));
1048 p_name.
help=
"Model name (default \"\").";
1049 cl.
par_list.insert(make_pair(
"name",&p_name));
1053 p_reference.
help=
"Model reference (default \"\").";
1054 cl.
par_list.insert(make_pair(
"reference",&p_reference));
1058 p_t0hc.
help=
"Model parameter t0 in MeV.";
1059 cl.
par_list.insert(make_pair(
"t0hc",&p_t0hc));
1063 p_t1hc.
help=
"Model parameter t1 in MeV.";
1064 cl.
par_list.insert(make_pair(
"t1hc",&p_t1hc));
1068 p_t2hc.
help=
"Model parameter t2 in MeV.";
1069 cl.
par_list.insert(make_pair(
"t2hc",&p_t2hc));
1073 p_t3hc.
help=
"Model parameter t3 in MeV.";
1074 cl.
par_list.insert(make_pair(
"t3hc",&p_t3hc));
1078 p_x0.
help=
"Model parameter x0.";
1079 cl.
par_list.insert(make_pair(
"x0",&p_x0));
1083 p_x1.
help=
"Model parameter x1.";
1084 cl.
par_list.insert(make_pair(
"x1",&p_x1));
1088 p_x2.
help=
"Model parameter x2.";
1089 cl.
par_list.insert(make_pair(
"x2",&p_x2));
1093 p_x3.
help=
"Model parameter x3.";
1094 cl.
par_list.insert(make_pair(
"x3",&p_x3));
1098 p_a.
help=
"Model parameter a.";
1099 cl.
par_list.insert(make_pair(
"a",&p_a));
1103 p_b.
help=
"Model parameter b.";
1104 cl.
par_list.insert(make_pair(
"b",&p_b));
1108 p_W0hc.
help=
"Model parameter W0hc.";
1109 cl.
par_list.insert(make_pair(
"W0hc",&p_W0hc));
1112 p_alpha.
d=&se.
sk.alpha;
1113 p_alpha.
help=
"Model parameter alpha.";
1114 cl.
par_list.insert(make_pair(
"alpha",&p_alpha));
1118 p_b4.
help=
"Model parameter b4.";
1119 cl.
par_list.insert(make_pair(
"b4",&p_b4));
1123 p_b4p.
help=
"Model parameter b4p.";
1124 cl.
par_list.insert(make_pair(
"b4p",&p_b4p));