23 #ifndef O2SCL_PART_DERIV_H
24 #define O2SCL_PART_DERIV_H
34 #include <o2scl/part.h>
35 #include <o2scl/fermion.h>
36 #include <o2scl/fermion_rel.h>
38 #ifndef DOXYGEN_NO_O2NS
45 template<
class fp_t=
double>
70 template<
class fp_t=
double>
123 template<
class fp_t=
double>
161 void deriv_f(fp_t &dmudn, fp_t &dmudT, fp_t &dsdT_n) {
175 template<
class fp_t=
double>
265 typedef fermion_deriv_tl<double> fermion_deriv;
269 template<
class fp_t=
double>
320 typedef part_deriv_tl<double> part_deriv;
415 template<
class part_deriv_t>
417 return (p.dsdT-p.dndT*p.dndT/p.dndmu)*temper/p.n;
464 template<
class part_deriv_t>
466 return temper/p.n*p.dsdT+p.en*p.en*temper/p.n/p.n/p.n*p.dndmu-
467 2.0*p.en*temper/p.n/p.n*p.dndT;
512 template<
class part_deriv_t>
514 return (p.dndT*p.dndT-p.dndmu*p.dsdT)/
515 (p.n*p.n*p.dsdT-2.0*p.n*p.en*p.dndT+p.en*p.en*p.dndmu);
548 template<
class part_deriv_t>
550 return p.dndmu/p.n/p.n;
589 template<
class part_deriv_t>
591 return p.en/p.n/p.n*p.dndmu-p.dndT/p.n;
682 template<
class part_deriv_t>
685 if (p.inc_rest_mass) {
690 return (p.n*p.n*p.dsdT-2.0*p.n*p.en*p.dndT+p.en*p.en*p.dndmu)/
691 (edt+p.pr)/(p.dndmu*p.dsdT-p.dndT*p.dndT);
696 typedef deriv_thermo_base_tl<double> deriv_thermo_base;
711 template<
class fp_t=
double>
731 pi=boost::math::constants::pi<fp_t>();
732 pi2=boost::math::constants::pi_sqr<fp_t>();
740 virtual int calc_mu(fermion_deriv &f, fp_t temper)=0;
744 virtual int calc_density(fermion_deriv &f, fp_t temper)=0;
749 virtual int pair_mu(fermion_deriv &f, fp_t temper)=0;
754 virtual int pair_density(fermion_deriv &f, fp_t temper)=0;
757 virtual int nu_from_n(fermion_deriv &f, fp_t temper)=0;
786 fp_t s2x=sqrt(2.0+x);
792 fp_t dndmu_term1=sx*s2x*(1.0+x)/f.
ms/f.
ms;
796 (-1.0+2.0*x*(2.0+x))/
797 f.
ms/f.
ms/sx/s2x/x/(2.0+x);
805 (1.0+x)/sx/s2x/x3/(x+2.0)/(x+2.0)/(x+2.0)/f.
ms/f.
ms;
807 f.
ms/f.
ms/pow(x*(2.0+x),2.5);
814 (3.0+2.0*x*(2.0+x))/f.
ms/f.
ms/pow(x*(2.0+x),5.5);
816 (7.0+6.0*x*(2.0+x))/f.
ms/f.
ms/pow(x*(2.0+x),4.5);
819 (1.0+x)/f.
ms/f.
ms/pow(x*(2.0+x),3.5);
822 f.
dndmu=prefac*(dndmu_term1+dndmu_term2+dndmu_term3+dndmu_term4);
823 f.
dndT=prefac*(dndT_term2+dndT_term3+dndT_term4);
824 f.
dsdT=prefac*(dsdT_term2+dsdT_term3+dsdT_term4);
838 fp_t prec,
bool inc_antip=
false) {
850 psi_num=f.
nu+f.
m-f.
ms;
861 static const size_t max_term=200;
865 fp_t first_dndmu=0.0;
874 for(
size_t j=1;j<=max_term;j++) {
882 fp_t pterm, nterm, enterm;
883 fp_t dndmu_term, dndT_term, dsdT_term;
888 if (inc_antip==
false) {
889 dndmu_term=nterm*jot;
890 dndT_term=jot*enterm-nterm/tt;
891 dsdT_term=(3.0*tt-2.0*dj*xx-2.0*dj)/tt/tt*enterm+
892 (5.0*dj*tt-2.0*dj*dj*xx+5.0*dj*tt*xx-dj*dj*xx*xx)/
895 dndmu_term=nterm*jot;
896 dndT_term=jot*enterm*tanh(jot*(xx+1.0))-
897 (tt+2.0*dj*(1.0+xx))/sinh(jot*(xx+1.0))*nterm/tt/tt;
898 dsdT_term=(2.0*dj*(1.0+xx)*tanh(jot*(xx+1.0))-3.0*tt)*enterm/tt/tt+
899 (2.0*pow(dj*1.0+xx,2.0)*tanh(jot*(xx+1.0))-
900 dj*dj*(2.0+2.0*xx+xx*xx)*cosh(jot*(xx+1.0))-
901 5.0*dj*(1.0+xx)*tt)*nterm/dj/tt/tt/tt;
908 first_dndT=dndT_term;
909 first_dsdT=dsdT_term;
910 first_dndmu=dndmu_term;
925 if (first_dndT==0.0 && first_dndmu==0.0 && first_dsdT==0.0) {
935 fabs(dndT_term)<prec*fabs(first_dndT) &&
936 fabs(dndmu_term)<prec*fabs(first_dndmu) &&
937 fabs(dsdT_term)<prec*fabs(first_dsdT)) {
953 typedef fermion_deriv_thermo_tl<double> fermion_deriv_thermo;
974 template<
class part_t,
class thermo_t>
976 std::string file,
bool nr_mode=
false,
977 int verbose=0,
bool external=
false) {
990 if (external==
false) {
997 std::cout <<
"In part_calibrate(), loading file named\n\t'"
998 << fname <<
"'.\n" << std::endl;
1004 #ifndef O2SCL_NO_HDF_INPUT
1010 std::string str=
"Failed to load data from file '"+fname+
1011 "' in part_calibrate(). Bad filename?";
1027 part_t bad, dev, exact;
1028 double m_bad=0.0, mu_bad=0.0, T_bad=0.0, mot_bad=0.0, psi_bad=0.0;
1029 p.non_interacting=
true;
1036 for(
size_t k=0;k<4;k++) {
1038 double ret_local=0.0;
1041 dev.n=0.0; dev.dndmu=0.0; dev.dndT=0.0; dev.dsdT=0.0;
1042 bad.n=0.0; bad.dndmu=0.0; bad.dndT=0.0; bad.dsdT=0.0;
1045 for(
double T=1.0e-2;T<=1.001e2;T*=1.0e2) {
1050 double mot=tab.
get(
"mot",i);
1051 double psi=tab.
get(
"psi",i);
1052 exact.n=tab.
get(
"n",i);
1053 exact.dndmu=tab.
get(
"dndmu",i);
1054 exact.dndT=tab.
get(
"dndT",i);
1055 exact.dsdT=tab.
get(
"dsdT",i);
1062 exact.n*=pow(T,3.0);
1063 exact.dndmu*=pow(T,2.0);
1064 exact.dndT*=pow(T,2.0);
1065 exact.dsdT*=pow(T,2.0);
1067 dev.dndmu+=fabs((p.dndmu-exact.dndmu)/exact.dndmu);
1068 dev.dndT+=fabs((p.dndT-exact.dndT)/exact.dndT);
1069 dev.dsdT+=fabs((p.dsdT-exact.dsdT)/exact.dsdT);
1073 check_derivs<part_t>(p,exact,bad,k,T,mot,
psi,mu_bad,m_bad,T_bad,
1074 mot_bad,psi_bad,ret_local);
1077 std::cout.precision(5);
1079 std::cout <<
"T,ms,nu,psi,mot: " << T <<
" "
1080 << p.ms <<
" " << p.nu
1081 <<
" " <<
psi <<
" " << mot << std::endl;
1083 std::cout <<
"T,m,mu,psi,mot: " << T <<
" "
1084 << p.m <<
" " << p.mu
1085 <<
" " <<
psi <<
" " << mot << std::endl;
1087 std::cout.precision(5);
1088 std::cout <<
"n,dndmu,dndT,dsdT: " << std::endl;
1089 std::cout <<
"approx: " << p.n <<
" " << p.dndmu <<
" "
1091 << p.dsdT << std::endl;
1092 std::cout <<
"exact : " << exact.n <<
" " << exact.dndmu <<
" "
1093 << exact.dndT <<
" " << exact.dsdT << std::endl;
1094 std::cout <<
"bad : " << bad.n <<
" " << bad.dndmu <<
" "
1095 << bad.dndT <<
" " << bad.dsdT << std::endl;
1096 std::cout << std::endl;
1103 if (ret_local>ret) {
1119 std::cout <<
"Function calc_mu(), include rest mass:" << std::endl;
1121 std::cout <<
"Function calc_mu(), without rest mass:" << std::endl;
1123 std::cout <<
"Function calc_mu(), include rest mass, "
1124 <<
"interacting:" << std::endl;
1126 std::cout <<
"Function calc_mu(), without rest mass, "
1127 <<
"interacting:" << std::endl;
1130 std::cout <<
"Average performance: " << std::endl;
1131 std::cout <<
"dndmu: " << dev.dndmu <<
" dndT: "
1132 << dev.dndT <<
" dsdT: " << dev.dsdT << std::endl;
1133 std::cout <<
"Worst case: " << std::endl;
1134 std::cout <<
"dndmu: " << bad.dndmu <<
" dndT: "
1135 << bad.dndT <<
" dsdT: " << bad.dsdT << std::endl;
1136 std::cout <<
"mu: " << mu_bad <<
" m: " << m_bad <<
" T: " << T_bad
1137 <<
" mot: " << mot_bad
1138 <<
"\n\tpsi: " << psi_bad << std::endl;
1139 std::cout << std::endl;
1147 p.non_interacting=
true;
1157 for(
size_t k=0;k<4;k++) {
1159 double ret_local=0.0;
1162 dev.mu=0.0; dev.dndmu=0.0; dev.dndT=0.0; dev.dsdT=0.0;
1163 bad.mu=0.0; bad.dndmu=0.0; bad.dndT=0.0; bad.dsdT=0.0;
1166 for(
double T=1.0e-2;T<=1.001e2;T*=1.0e2) {
1171 double mot=tab.
get(
"mot",i);
1172 double psi=tab.
get(
"psi",i);
1174 exact.dndmu=tab.
get(
"dndmu",i);
1175 exact.dndT=tab.
get(
"dndT",i);
1176 exact.dsdT=tab.
get(
"dsdT",i);
1182 exact.dndmu*=pow(T,2.0);
1183 exact.dndT*=pow(T,2.0);
1184 exact.dsdT*=pow(T,2.0);
1193 th.calc_density(p,T);
1195 dev.dndmu+=fabs((p.dndmu-exact.dndmu)/exact.dndmu);
1196 dev.dndT+=fabs((p.dndT-exact.dndT)/exact.dndT);
1197 dev.dsdT+=fabs((p.dsdT-exact.dsdT)/exact.dsdT);
1201 check_derivs<part_t>(p,exact,bad,k,T,mot,
psi,mu_bad,m_bad,
1202 T_bad,mot_bad,psi_bad,ret_local);
1205 std::cout.precision(5);
1207 std::cout <<
"T,ms,n,psi,mot: " << T <<
" "
1208 << p.ms <<
" " << p.n <<
" "
1209 <<
psi <<
" " << mot << std::endl;
1211 std::cout <<
"T,m,n,psi,mot: " << T <<
" " << p.m <<
" "
1212 << p.n <<
" " <<
psi <<
" " << mot << std::endl;
1214 std::cout.precision(6);
1215 std::cout <<
"mu,dndmu,dndT,dsdT: " << std::endl;
1216 std::cout <<
"approx: " << p.mu <<
" "
1217 << p.dndmu <<
" " << p.dndT <<
" "
1218 << p.dsdT << std::endl;
1219 std::cout <<
"exact : " << exact.mu <<
" "
1220 << exact.dndmu <<
" " << exact.dndT <<
" "
1221 << exact.dsdT << std::endl;
1222 std::cout <<
"bad : " << bad.mu <<
" " << bad.dndmu <<
" "
1223 << bad.dndT <<
" " << bad.dsdT << std::endl;
1224 std::cout << std::endl;
1231 if (ret_local>ret) {
1247 std::cout <<
"Function calc_density(), include rest mass:"
1250 std::cout <<
"Function calc_density(), without rest mass:"
1253 std::cout <<
"Function calc_density(), include "
1254 <<
"rest mass, interacting:"
1257 std::cout <<
"Function calc_density(), without rest mass, "
1262 std::cout <<
"Average performance: " << std::endl;
1263 std::cout <<
"dndmu: " << dev.dndmu <<
" dndT: "
1264 << dev.dndT <<
" dsdT: " << dev.dsdT << std::endl;
1265 std::cout <<
"Worst case: " << std::endl;
1266 std::cout <<
"dndmu: " << bad.dndmu <<
" dndT: "
1267 << bad.dndT <<
" dsdT: " << bad.dsdT << std::endl;
1268 std::cout <<
"mu: " << mu_bad <<
" m: " << m_bad
1269 <<
" T: " << T_bad <<
" mot: " << mot_bad
1270 <<
"\n\tpsi: " << psi_bad << std::endl;
1271 std::cout << std::endl;
1288 for(
size_t k=0;k<4;k++) {
1290 double ret_local=0.0;
1293 dev.n=0.0; dev.dndmu=0.0; dev.dndT=0.0; dev.dsdT=0.0;
1294 bad.n=0.0; bad.dndmu=0.0; bad.dndT=0.0; bad.dsdT=0.0;
1297 for(
double T=1.0e-2;T<=1.001e2;T*=1.0e2) {
1302 double mot=tab.
get(
"mot",i);
1303 double psi=tab.
get(
"psi",i);
1304 exact.n=tab.
get(
"pair_n",i);
1305 exact.dndmu=tab.
get(
"pair_dndmu",i);
1306 exact.dndT=tab.
get(
"pair_dndT",i);
1307 exact.dsdT=tab.
get(
"pair_dsdT",i);
1314 exact.n*=pow(T,3.0);
1315 exact.dndmu*=pow(T,2.0);
1316 exact.dndT*=pow(T,2.0);
1317 exact.dsdT*=pow(T,2.0);
1319 dev.dndmu+=fabs((p.dndmu-exact.dndmu)/exact.dndmu);
1320 dev.dndT+=fabs((p.dndT-exact.dndT)/exact.dndT);
1321 dev.dsdT+=fabs((p.dsdT-exact.dsdT)/exact.dsdT);
1325 check_derivs<part_t>(p,exact,bad,k,T,mot,
psi,mu_bad,m_bad,
1326 T_bad,mot_bad,psi_bad,ret_local);
1329 std::cout.precision(5);
1330 std::cout <<
"T,m,mu,psi,mot: " << T <<
" " << p.m
1332 <<
" " <<
psi <<
" " << mot << std::endl;
1333 std::cout.precision(6);
1334 std::cout <<
"n,dndmu,dndT,dsdT: " << std::endl;
1335 std::cout <<
"approx: " << p.n <<
" " << p.dndmu <<
" "
1337 << p.dsdT << std::endl;
1338 std::cout <<
"exact : " << exact.n <<
" "
1339 << exact.dndmu <<
" " << exact.dndT <<
" "
1340 << exact.dsdT << std::endl;
1341 std::cout <<
"bad : " << bad.n <<
" " << bad.dndmu <<
" "
1342 << bad.dndT <<
" " << bad.dsdT << std::endl;
1343 std::cout << std::endl;
1350 if (ret_local>ret) {
1366 std::cout <<
"Function pair_mu(), include rest mass"
1369 std::cout <<
"Function pair_mu(), without rest mass"
1372 std::cout <<
"Function pair_mu(), include rest mass, "
1373 <<
"interacting" << std::endl;
1375 std::cout <<
"Function pair_mu(), without rest mass, "
1376 <<
"interacting" << std::endl;
1379 std::cout <<
"Average performance: " << std::endl;
1380 std::cout <<
"dndmu: " << dev.dndmu <<
" dndT: "
1381 << dev.dndT <<
" dsdT: " << dev.dsdT << std::endl;
1382 std::cout <<
"Worst case: " << std::endl;
1383 std::cout <<
"dndmu: " << bad.dndmu <<
" dndT: "
1384 << bad.dndT <<
" dsdT: " << bad.dsdT << std::endl;
1385 std::cout <<
"mu: " << mu_bad <<
" m: " << m_bad
1386 <<
" T: " << T_bad <<
" mot: " << mot_bad
1387 <<
"\n\tpsi: " << psi_bad << std::endl;
1388 std::cout << std::endl;
1403 for(
size_t k=0;k<4;k++) {
1405 double ret_local=0.0;
1408 dev.mu=0.0; dev.dndmu=0.0; dev.dndT=0.0; dev.dsdT=0.0;
1409 bad.mu=0.0; bad.dndmu=0.0; bad.dndT=0.0; bad.dsdT=0.0;
1412 for(
double T=1.0e-2;T<=1.001e2;T*=1.0e2) {
1417 double mot=tab.
get(
"mot",i);
1418 double psi=tab.
get(
"psi",i);
1419 p.n=tab.
get(
"pair_n",i);
1420 exact.dndmu=tab.
get(
"pair_dndmu",i);
1421 exact.dndT=tab.
get(
"pair_dndT",i);
1422 exact.dsdT=tab.
get(
"pair_dsdT",i);
1428 exact.dndmu*=pow(T,2.0);
1429 exact.dndT*=pow(T,2.0);
1430 exact.dsdT*=pow(T,2.0);
1435 th.pair_density(p,T);
1437 dev.dndmu+=fabs((p.dndmu-exact.dndmu)/exact.dndmu);
1438 dev.dndT+=fabs((p.dndT-exact.dndT)/exact.dndT);
1439 dev.dsdT+=fabs((p.dsdT-exact.dsdT)/exact.dsdT);
1443 check_derivs<part_t>(p,exact,bad,k,T,mot,
psi,mu_bad,m_bad,
1444 T_bad,mot_bad,psi_bad,ret_local);
1447 std::cout.precision(5);
1448 std::cout <<
"T,m,n,psi,mot: " << T <<
" " << p.m <<
" "
1449 << p.n <<
" " <<
psi <<
" " << mot << std::endl;
1450 std::cout.precision(6);
1451 std::cout <<
"mu,dndmu,dndT,dsdT: " << std::endl;
1452 std::cout <<
"approx: " << p.mu <<
" " << p.dndmu <<
" "
1453 << p.dndT <<
" " << p.dsdT << std::endl;
1454 std::cout <<
"exact : " << exact.mu <<
" "
1455 << exact.dndmu <<
" " << exact.dndT <<
" "
1456 << exact.dsdT << std::endl;
1457 std::cout <<
"bad : " << bad.mu <<
" " << bad.dndmu <<
" "
1458 << bad.dndT <<
" " << bad.dsdT << std::endl;
1459 std::cout << std::endl;
1466 if (ret_local>ret) {
1482 std::cout <<
"Function pair_density(), include rest mass"
1485 std::cout <<
"Function pair_density(), without rest mass"
1488 std::cout <<
"Function pair_density(), include rest mass, "
1489 <<
"interacting" << std::endl;
1491 std::cout <<
"Function pair_density(), without rest mass, "
1492 <<
"interacting" << std::endl;
1495 std::cout <<
"Average performance: " << std::endl;
1496 std::cout <<
"dndmu: "
1497 << dev.dndmu <<
" dndT: "
1498 << dev.dndT <<
" dsdT: " << dev.dsdT << std::endl;
1499 std::cout <<
"Worst case: " << std::endl;
1500 std::cout <<
"dndmu: " << bad.dndmu
1501 <<
" dndT: " << bad.dndT
1502 <<
" dsdT: " << bad.dsdT << std::endl;
1503 std::cout <<
"mu: " << mu_bad <<
" m: " << m_bad
1504 <<
" T: " << T_bad <<
" mot: " << mot_bad
1505 <<
"\n\tpsi: " << psi_bad << std::endl;
1506 std::cout << std::endl;
1529 #ifndef DOXYGEN_NO_O2NS