29 #include <o2scl/constants.h>
30 #include <o2scl/inte.h>
31 #include <o2scl/funct.h>
32 #include <o2scl/mroot.h>
35 #include <o2scl/lib_settings.h>
37 #include <o2scl/table.h>
38 #include <o2scl/hdf_file.h>
39 #include <o2scl/hdf_io.h>
45 #ifndef DOXYGEN_NO_O2NS
64 const char *
type() {
return "thermo"; }
179 virtual void init(fp_t mass, fp_t dof) {
222 virtual const char *
type() {
return "part_tl"; }
247 template<
class part_t>
250 p.non_interacting=
false;
254 p.non_interacting=
true;
259 p.inc_rest_mass=
true;
261 p.inc_rest_mass=
false;
269 template<
class part_t>
303 template<
class part1_t,
class part2_t,
class part3_t>
305 double T,
double mot,
double psi,
306 double &mu_bad,
double &m_bad,
307 double &T_bad,
double &mot_bad,
double &psi_bad,
309 if (fabs((p.n-exact.n)/exact.n)>bad.n) {
310 bad.n=fabs((p.n-exact.n)/exact.n);
311 if (bad.n>ret_local) {
331 template<
class part1_t,
class part2_t,
class part3_t>
333 double T,
double mot,
double psi,
334 double &mu_bad,
double &m_bad,
335 double &T_bad,
double &mot_bad,
double &psi_bad,
338 if (fabs((p.nu-exact.nu)/exact.nu)>bad.mu) {
339 bad.mu=fabs((p.nu-exact.nu)/exact.nu);
340 if (bad.mu>ret_local) {
350 if (fabs((p.mu-exact.mu)/exact.mu)>bad.mu) {
351 bad.mu=fabs((p.mu-exact.mu)/exact.mu);
352 if (bad.mu>ret_local) {
368 template<
class part1_t,
class part2_t,
class part3_t>
369 void check_eps(part1_t &p, part2_t &exact, part3_t &bad,
size_t k,
370 double T,
double mot,
double psi,
371 double &mu_bad,
double &m_bad,
372 double &T_bad,
double &mot_bad,
double &psi_bad,
374 if (fabs((p.ed-exact.ed)/exact.ed)>bad.ed) {
375 bad.ed=fabs((p.ed-exact.ed)/exact.ed);
376 if (bad.ed>ret_local) {
390 if (fabs((p.pr-exact.pr)/exact.pr)>bad.pr) {
391 bad.pr=fabs((p.pr-exact.pr)/exact.pr);
392 if (bad.pr>ret_local) {
406 if (fabs((p.en-exact.en)/exact.en)>bad.en) {
407 bad.en=fabs((p.en-exact.en)/exact.en);
408 if (bad.en>ret_local) {
428 template<
class part1_t,
class part2_t,
class part3_t>
430 double T,
double mot,
double psi,
431 double &mu_bad,
double &m_bad,
432 double &T_bad,
double &mot_bad,
double &psi_bad,
434 if (fabs((p.dndT-exact.dndT)/exact.dndT)>bad.dndT) {
435 bad.dndT=fabs((p.dndT-exact.dndT)/exact.dndT);
436 if (bad.dndT>ret_local) {
450 if (fabs((p.dndmu-exact.dndmu)/exact.dndmu)>bad.dndmu) {
451 bad.dndmu=fabs((p.dndmu-exact.dndmu)/exact.dndmu);
452 if (bad.dndmu>ret_local) {
466 if (fabs((p.dsdT-exact.dsdT)/exact.dsdT)>bad.dsdT) {
467 bad.dsdT=fabs((p.dsdT-exact.dsdT)/exact.dsdT);
468 if (bad.dsdT>ret_local) {
499 template<
class part_t,
class thermo_t>
501 std::string file,
bool nr_mode=
false,
502 int verbose=0,
bool external=
false) {
515 if (external==
false) {
522 std::cout <<
"In part_calibrate(), loading file named\n\t'"
523 << fname <<
"'.\n" << std::endl;
529 #ifndef O2SCL_NO_HDF_INPUT
535 std::string str=
"Failed to load data from file '"+fname+
536 "' in part_calibrate(). Bad filename?";
552 part bad, dev, exact;
553 double m_bad=0.0, mu_bad=0.0, T_bad=0.0, mot_bad=0.0, psi_bad=0.0;
561 for(
size_t k=0;k<4;k++) {
563 double ret_local=0.0;
566 dev.
n=0.0; dev.
ed=0.0; dev.
pr=0.0; dev.
en=0.0;
567 bad.
n=0.0; bad.
ed=0.0; bad.
pr=0.0; bad.
en=0.0;
570 for(
double T=1.0e-2;T<=1.001e2;T*=1.0e2) {
575 double mot=tab.
get(
"mot",i);
576 double psi=tab.
get(
"psi",i);
577 exact.
n=tab.
get(
"n",i);
578 exact.
ed=tab.
get(
"ed",i);
579 exact.
pr=tab.
get(
"pr",i);
580 exact.
en=tab.
get(
"en",i);
590 exact.
ed=exact.
ed*pow(T,4.0)+exact.
n*p.m;
592 exact.
ed=exact.
ed*pow(T,4.0);
596 exact.
ed*=pow(T,4.0);
598 exact.
ed=exact.
ed*pow(T,4.0)-exact.
n*p.m;
601 exact.
pr*=pow(T,4.0);
602 exact.
en*=pow(T,3.0);
604 dev.
n+=fabs((p.n-exact.
n)/exact.
n);
605 dev.
ed+=fabs((p.ed-exact.
ed)/exact.
ed);
606 dev.
pr+=fabs((p.pr-exact.
pr)/exact.
pr);
607 dev.
en+=fabs((p.en-exact.
en)/exact.
en);
611 check_density<part_t>(p,exact,bad,k,T,mot,
psi,mu_bad,m_bad,T_bad,
612 mot_bad,psi_bad,ret_local);
614 check_eps<part_t>(p,exact,bad,k,T,mot,
psi,mu_bad,m_bad,T_bad,
615 mot_bad,psi_bad,ret_local);
618 std::cout.precision(5);
620 std::cout <<
"T,ms,nu,psi,mot: " << T <<
" "
621 << p.ms <<
" " << p.nu <<
" "
622 <<
psi <<
" " << mot << std::endl;
624 std::cout <<
"T,m,mu,psi,mot: " << T <<
" "
625 << p.m <<
" " << p.mu <<
" "
626 <<
psi <<
" " << mot << std::endl;
628 std::cout.precision(5);
629 std::cout <<
"n,ed,pr,en: " << std::endl;
630 std::cout <<
"approx: " << p.n <<
" " << p.ed <<
" "
631 << p.pr <<
" " << p.en << std::endl;
632 std::cout <<
"exact : " << exact.
n <<
" " << exact.
ed <<
" "
633 << exact.
pr <<
" " << exact.
en << std::endl;
634 std::cout <<
"bad : " << bad.
n <<
" " << bad.
ed <<
" "
635 << bad.
pr <<
" " << bad.
en << std::endl;
636 std::cout <<
"ret_local,ret: " << ret_local <<
" "
638 std::cout << std::endl;
661 std::cout <<
"Function calc_mu(), include rest mass"
664 std::cout <<
"Function calc_mu(), without rest mass"
667 std::cout <<
"Function calc_mu(), include rest mass, "
668 <<
"interacting" << std::endl;
670 std::cout <<
"Function calc_mu(), without rest mass, "
671 <<
"interacting" << std::endl;
674 std::cout <<
"Average performance: " << std::endl;
675 std::cout <<
"n: " << dev.
n <<
" ed: " << dev.
ed <<
" pr: "
676 << dev.
pr <<
" en: " << dev.
en << std::endl;
677 std::cout <<
"Worst case: " << std::endl;
678 std::cout <<
"n: " << bad.
n <<
" ed: " << bad.
ed <<
" pr: "
679 << bad.
pr <<
" en: " << bad.
en << std::endl;
680 std::cout <<
"mu: " << mu_bad <<
" m: " << m_bad
681 <<
" T: " << T_bad <<
" mot: " << mot_bad
682 <<
"\n\tpsi: " << psi_bad << std::endl;
683 std::cout <<
"ret_local,ret: " << ret_local <<
" "
685 std::cout << std::endl;
693 p.non_interacting=
true;
703 for(
size_t k=0;k<4;k++) {
705 double ret_local=0.0;
708 dev.
mu=0.0; dev.
ed=0.0; dev.
pr=0.0; dev.
en=0.0;
709 bad.
mu=0.0; bad.
ed=0.0; bad.
pr=0.0; bad.
en=0.0;
712 for(
double T=1.0e-2;T<=1.001e2;T*=1.0e2) {
717 double mot=tab.
get(
"mot",i);
718 double psi=tab.
get(
"psi",i);
719 p.n=tab.
get(
"n",i)*pow(T,3.0);
720 exact.
ed=tab.
get(
"ed",i);
721 exact.
pr=tab.
get(
"pr",i);
722 exact.
en=tab.
get(
"en",i);
730 exact.
ed=exact.
ed*pow(T,4.0)+exact.
n*p.m;
732 exact.
ed=exact.
ed*pow(T,4.0);
736 exact.
ed*=pow(T,4.0);
738 exact.
ed=exact.
ed*pow(T,4.0)-exact.
n*p.m;
741 exact.
pr*=pow(T,4.0);
742 exact.
en*=pow(T,3.0);
751 th.calc_density(p,T);
754 dev.
nu+=fabs((p.nu-exact.
nu)/exact.
nu);
756 dev.
mu+=fabs((p.mu-exact.
mu)/exact.
mu);
758 dev.
ed+=fabs((p.ed-exact.
ed)/exact.
ed);
759 dev.
pr+=fabs((p.pr-exact.
pr)/exact.
pr);
760 dev.
en+=fabs((p.en-exact.
en)/exact.
en);
764 check_chem_pot<part_t>(p,exact,bad,k,T,mot,
psi,mu_bad,m_bad,T_bad,
765 mot_bad,psi_bad,ret_local);
767 check_eps<part_t>(p,exact,bad,k,T,mot,
psi,mu_bad,m_bad,T_bad,
768 mot_bad,psi_bad,ret_local);
771 std::cout.precision(5);
773 std::cout <<
"T,ms,n,psi,mot: " << T <<
" "
774 << p.ms <<
" " << p.n <<
" "
775 <<
psi <<
" " << mot << std::endl;
777 std::cout <<
"T,m,n,psi,mot: " << T <<
" "
778 << p.m <<
" " << p.n <<
" "
779 <<
psi <<
" " << mot << std::endl;
781 std::cout.precision(6);
783 std::cout <<
"nu,ed,pr,en: " << std::endl;
784 std::cout <<
"approx: " << p.nu <<
" " << p.ed <<
" "
785 << p.pr <<
" " << p.en << std::endl;
786 std::cout <<
"exact : " << exact.
nu <<
" " << exact.
ed <<
" "
787 << exact.
pr <<
" " << exact.
en << std::endl;
789 std::cout <<
"mu,ed,pr,en: " << std::endl;
790 std::cout <<
"approx: " << p.mu <<
" " << p.ed <<
" "
791 << p.pr <<
" " << p.en << std::endl;
792 std::cout <<
"exact : " << exact.
mu <<
" " << exact.
ed <<
" "
793 << exact.
pr <<
" " << exact.
en << std::endl;
795 std::cout <<
"bad : " << bad.
mu <<
" " << bad.
ed <<
" "
796 << bad.
pr <<
" " << bad.
en << std::endl;
797 std::cout <<
"ret_local,ret: " << ret_local <<
" "
799 std::cout << std::endl;
822 std::cout <<
"Function calc_density(), include rest mass"
825 std::cout <<
"Function calc_density(), without rest mass"
828 std::cout <<
"Function calc_density(), include rest mass, "
829 <<
"interacting" << std::endl;
831 std::cout <<
"Function calc_density(), without rest mass, "
832 <<
"interacting" << std::endl;
835 std::cout <<
"Average performance: " << std::endl;
836 std::cout <<
"mu: " << dev.
mu <<
" ed: " << dev.
ed <<
" pr: "
837 << dev.
pr <<
" en: " << dev.
en << std::endl;
838 std::cout <<
"Worst case: " << std::endl;
839 std::cout <<
"mu: " << bad.
mu <<
" ed: " << bad.
ed <<
" pr: "
840 << bad.
pr <<
" en: " << bad.
en << std::endl;
841 std::cout <<
"mu: " << mu_bad <<
" m: " << m_bad
842 <<
" T: " << T_bad <<
" mot: " << mot_bad
843 <<
"\n\tpsi: " << psi_bad << std::endl;
844 std::cout <<
"ret_local,ret: " << ret_local <<
" "
846 std::cout << std::endl;
863 for(
size_t k=0;k<4;k++) {
865 double ret_local=0.0;
868 dev.
n=0.0; dev.
ed=0.0; dev.
pr=0.0; dev.
en=0.0;
869 bad.
n=0.0; bad.
ed=0.0; bad.
pr=0.0; bad.
en=0.0;
872 for(
double T=1.0e-2;T<=1.001e2;T*=1.0e2) {
877 double mot=tab.
get(
"mot",i);
878 double psi=tab.
get(
"psi",i);
879 exact.
n=tab.
get(
"pair_n",i);
880 exact.
ed=tab.
get(
"pair_ed",i);
881 exact.
pr=tab.
get(
"pair_pr",i);
882 exact.
en=tab.
get(
"pair_en",i);
891 exact.
ed*=pow(T,4.0);
893 exact.
ed=exact.
ed*pow(T,4.0)-exact.
n*p.m;
895 exact.
pr*=pow(T,4.0);
896 exact.
en*=pow(T,3.0);
898 dev.
n+=fabs((p.n-exact.
n)/exact.
n);
899 dev.
ed+=fabs((p.ed-exact.
ed)/exact.
ed);
900 dev.
pr+=fabs((p.pr-exact.
pr)/exact.
pr);
901 dev.
en+=fabs((p.en-exact.
en)/exact.
en);
905 check_density<part_t>(p,exact,bad,k,T,mot,
psi,mu_bad,m_bad,T_bad,
906 mot_bad,psi_bad,ret_local);
908 check_eps<part_t>(p,exact,bad,k,T,mot,
psi,mu_bad,m_bad,T_bad,
909 mot_bad,psi_bad,ret_local);
912 std::cout.precision(5);
913 std::cout <<
"T,m,mu,psi,mot: " << T <<
" " << p.m <<
" "
914 << p.mu <<
" " <<
psi <<
" " << mot << std::endl;
915 std::cout.precision(6);
916 std::cout << i <<
" " << exact.
m <<
" " << exact.
ms <<
" "
917 << exact.
mu <<
" " << exact.
nu << std::endl;
918 std::cout <<
"n,ed,pr,en: " << std::endl;
919 std::cout <<
"approx: " << p.n <<
" " << p.ed <<
" "
920 << p.pr <<
" " << p.en << std::endl;
921 std::cout <<
"exact : " << exact.
n <<
" " << exact.
ed <<
" "
922 << exact.
pr <<
" " << exact.
en << std::endl;
923 std::cout <<
"bad : " << bad.
n <<
" " << bad.
ed <<
" "
924 << bad.
pr <<
" " << bad.
en << std::endl;
925 std::cout <<
"ret_local,ret: " << ret_local <<
" "
927 std::cout << std::endl;
950 std::cout <<
"Function pair_mu(), include rest mass"
953 std::cout <<
"Function pair_mu(), without rest mass"
956 std::cout <<
"Function pair_mu(), include rest mass, "
957 <<
"interacting" << std::endl;
959 std::cout <<
"Function pair_mu(), without rest mass, "
960 <<
"interacting" << std::endl;
963 std::cout <<
"Average performance: " << std::endl;
964 std::cout <<
"n: " << dev.
n <<
" ed: " << dev.
ed <<
" pr: "
965 << dev.
pr <<
" en: " << dev.
en << std::endl;
966 std::cout <<
"Worst case: " << std::endl;
967 std::cout <<
"n: " << bad.
n <<
" ed: " << bad.
ed <<
" pr: "
968 << bad.
pr <<
" en: " << bad.
en << std::endl;
969 std::cout <<
"mu: " << mu_bad <<
" m: " << m_bad
970 <<
" T: " << T_bad <<
" mot: " << mot_bad
971 <<
"\n\tpsi: " << psi_bad << std::endl;
972 std::cout <<
"ret_local,ret: " << ret_local <<
" "
974 std::cout << std::endl;
989 for(
size_t k=0;k<4;k++) {
991 double ret_local=0.0;
994 dev.
mu=0.0; dev.
ed=0.0; dev.
pr=0.0; dev.
en=0.0;
995 bad.
mu=0.0; bad.
ed=0.0; bad.
pr=0.0; bad.
en=0.0;
998 for(
double T=1.0e-2;T<=1.001e2;T*=1.0e2) {
1003 double mot=tab.
get(
"mot",i);
1004 double psi=tab.
get(
"psi",i);
1005 p.n=tab.
get(
"pair_n",i)*pow(T,3.0);
1006 exact.
ed=tab.
get(
"pair_ed",i);
1007 exact.
pr=tab.
get(
"pair_pr",i);
1008 exact.
en=tab.
get(
"pair_en",i);
1015 exact.
ed*=pow(T,4.0);
1017 exact.
ed=exact.
ed*pow(T,4.0)-exact.
n*p.m;
1019 exact.
pr*=pow(T,4.0);
1020 exact.
en*=pow(T,3.0);
1025 th.pair_density(p,T);
1028 dev.
nu+=fabs((p.nu-exact.
nu)/exact.
nu);
1030 dev.
mu+=fabs((p.mu-exact.
mu)/exact.
mu);
1032 dev.
ed+=fabs((p.ed-exact.
ed)/exact.
ed);
1033 dev.
pr+=fabs((p.pr-exact.
pr)/exact.
pr);
1034 dev.
en+=fabs((p.en-exact.
en)/exact.
en);
1038 check_chem_pot<part_t>(p,exact,bad,k,T,mot,
psi,mu_bad,m_bad,T_bad,
1039 mot_bad,psi_bad,ret_local);
1041 check_eps<part_t>(p,exact,bad,k,T,mot,
psi,mu_bad,m_bad,T_bad,
1042 mot_bad,psi_bad,ret_local);
1046 std::cout.precision(5);
1048 std::cout <<
"T,ms,n,psi,mot: " << T <<
" " << p.ms <<
" "
1049 << p.n <<
" " <<
psi <<
" " << mot << std::endl;
1051 std::cout <<
"T,m,n,psi,mot: " << T <<
" " << p.m <<
" "
1052 << p.n <<
" " <<
psi <<
" " << mot << std::endl;
1054 std::cout.precision(6);
1056 std::cout <<
"nu,ed,pr,en: " << std::endl;
1057 std::cout <<
"approx: " << p.nu <<
" " << p.ed <<
" "
1058 << p.pr <<
" " << p.en << std::endl;
1059 std::cout <<
"exact : " << exact.
nu <<
" " << exact.
ed <<
" "
1060 << exact.
pr <<
" " << exact.
en << std::endl;
1062 std::cout <<
"mu,ed,pr,en: " << std::endl;
1063 std::cout <<
"approx: " << p.mu <<
" " << p.ed <<
" "
1064 << p.pr <<
" " << p.en << std::endl;
1065 std::cout <<
"exact : " << exact.
mu <<
" " << exact.
ed <<
" "
1066 << exact.
pr <<
" " << exact.
en << std::endl;
1068 std::cout <<
"bad : " << bad.
mu <<
" " << bad.
ed <<
" "
1069 << bad.
pr <<
" " << bad.
en << std::endl;
1070 std::cout <<
"ret_local,ret: " << ret_local <<
" "
1071 << ret << std::endl;
1072 std::cout << std::endl;
1079 if (ret_local>ret) {
1095 std::cout <<
"Function pair_density(), include rest mass"
1098 std::cout <<
"Function pair_density(), without rest mass"
1101 std::cout <<
"Function pair_density(), include rest mass, "
1102 <<
"interacting" << std::endl;
1104 std::cout <<
"Function pair_density(), without rest mass, "
1105 <<
"interacting" << std::endl;
1108 std::cout <<
"Average performance: " << std::endl;
1109 std::cout <<
"mu: " << dev.
mu <<
" ed: " << dev.
ed <<
" pr: "
1110 << dev.
pr <<
" en: " << dev.
en << std::endl;
1111 std::cout <<
"Worst case: " << std::endl;
1112 std::cout <<
"mu: " << bad.
mu <<
" ed: " << bad.
ed <<
" pr: "
1113 << bad.
pr <<
" en: " << bad.
en << std::endl;
1114 std::cout <<
"mu: " << mu_bad <<
" m: " << m_bad
1115 <<
" T: " << T_bad <<
" mot: " << mot_bad
1116 <<
"\n\tpsi: " << psi_bad << std::endl;
1117 std::cout <<
"ret_local,ret: " << ret_local <<
" "
1118 << ret << std::endl;
1119 std::cout << std::endl;
1126 if (ret_local>ret) {
1146 #ifndef DOXYGEN_NO_O2NS