22 #ifndef O2SCL_FERMION_REL_H
23 #define O2SCL_FERMION_REL_H
33 #include <o2scl/constants.h>
34 #include <o2scl/mroot.h>
35 #include <o2scl/inte.h>
36 #include <o2scl/fermion.h>
37 #include <o2scl/root_brent_gsl.h>
38 #include <o2scl/inte_qagiu_gsl.h>
39 #include <o2scl/inte_qag_gsl.h>
41 #ifndef DOXYGEN_NO_O2NS
218 template<
class fp_t=
double>
284 calc_mu_tlate<fermion>(f,temper);
296 return calc_density_tlate<fermion>(f,temper);
303 pair_mu_tlate<fermion>(f,temper);
311 return pair_density_tlate<fermion>(f,temper);
317 return nu_from_n_tlate<fermion>(f,temper);
321 std::shared_ptr<inte<> >
nit;
324 std::shared_ptr<inte<> >
dit;
330 virtual const char *
type() {
return "fermion_rel"; }
390 template<
class fermion_t>
402 std::cout <<
"nu_from_n(): initial guess " << nex << std::endl;
407 if (temper>scale) scale=temper;
408 for(
size_t i=0;i<10;i++) {
409 if (nex<0.0) nex+=scale*1.0e5;
412 if (y<1.0-1.0e-6) i=10;
415 std::cout <<
"nu_from_n(): adjusted guess to " << nex << std::endl;
421 if (f.inc_rest_mass) {
424 nex=(f.ms-f.m)/temper;
428 std::cout <<
"nu_from_n(): adjusted guess (try 2) to "
434 if (y==1.0 || !std::isfinite(y)) {
441 funct mf=std::bind(std::mem_fn<fp_t(fp_t,
fermion &,fp_t)>
443 this,std::placeholders::_1,std::ref(f),temper);
457 std::cout <<
"nu_from_n(): density_root failed x="
458 << nex <<
" ." << std::endl;
459 std::cout <<
"\tTrying to make integrators more accurate."
464 fp_t tol1=
dit->tol_rel, tol2=
dit->tol_abs;
465 fp_t tol3=
nit->tol_rel, tol4=
nit->tol_abs;
475 std::cout <<
"nu_from_n(): density_root failed again x=" << nex
476 <<
" ." << std::endl;
477 std::cout <<
"Trying to bracket root." << std::endl;
480 fp_t lg=std::max(fabs(f.nu),f.ms);
481 fp_t bhigh=lg/temper, blow=-bhigh;
482 fp_t yhigh=mf(bhigh), ylow=mf(blow);
483 for(
size_t j=0;j<5 && yhigh>0.0;j++) {
487 for(
size_t j=0;j<5 && ylow<0.0;j++) {
491 if (yhigh<0.0 && ylow>0.0) {
501 std::cout <<
"nu_from_n(): density_root failed fourth solver "
502 << blow << std::endl;
506 std::cout <<
"nu_from_n(): Failed to bracket." << std::endl;
540 template<
class fermion_t>
558 if (f.non_interacting==
true) { f.nu=f.mu; f.ms=f.m; }
564 if (f.inc_rest_mass) {
565 psi=(f.nu-f.ms)/temper;
567 psi=(f.nu+(f.m-f.ms))/temper;
601 funct mfd=std::bind(std::mem_fn<fp_t(fp_t,
fermion &,fp_t)>
603 this,std::placeholders::_1,std::ref(f),temper);
604 funct mfe=std::bind(std::mem_fn<fp_t(fp_t,
fermion &,fp_t)>
606 this,std::placeholders::_1,std::ref(f),temper);
607 funct mfs=std::bind(std::mem_fn<fp_t(fp_t,
fermion &,fp_t)>
609 this,std::placeholders::_1,std::ref(f),temper);
615 f.n=
nit->integ(mfd,0.0,0.0);
621 f.ed=
nit->integ(mfe,0.0,0.0);
623 if (!f.inc_rest_mass) f.ed-=f.n*f.m;
624 unc.
ed=
nit->get_error()*prefac*temper;
628 f.en=
nit->integ(mfs,0.0,0.0);
639 funct mfd=std::bind(std::mem_fn<fp_t(fp_t,
fermion &,fp_t)>
641 this,std::placeholders::_1,std::ref(f),temper);
642 funct mfe=std::bind(std::mem_fn<fp_t(fp_t,
fermion &,fp_t)>
644 this,std::placeholders::_1,std::ref(f),temper);
645 funct mfs=std::bind(std::mem_fn<fp_t(fp_t,
fermion &,fp_t)>
647 this,std::placeholders::_1,std::ref(f),temper);
654 if (f.inc_rest_mass) {
671 O2SCL_ERR2(
"Zero density in degenerate limit in fermion_rel::",
672 "calc_mu(). Variable deg_limit set improperly?",
679 f.n=
dit->integ(mfd,0.0,ul);
685 f.ed=
dit->integ(mfe,0.0,ul);
692 if (f.inc_rest_mass) {
711 f.en=
dit->integ(mfs,ll,ul);
714 f.en=
dit->integ(mfs,0.0,ul);
724 f.pr=-f.ed+temper*f.en+f.nu*f.n;
737 template<
class fermion_t>
745 fp_t density_temp=f.n;
752 "fermion_rel::calc_density_tlate().",
761 #if !O2SCL_NO_RANGE_CHECK
766 if (!std::isfinite(f.n)) {
768 "fermion_rel::calc_density_tlate().",
exc_einval);
775 if (f.non_interacting==
true) { f.nu=f.mu; f.ms=f.m; }
784 if (f.non_interacting) { f.mu=f.nu; }
792 if (f.inc_rest_mass) {
793 psi=(f.nu-f.ms)/temper;
795 psi=(f.nu+(f.m-f.ms))/temper;
828 funct mfe=std::bind(std::mem_fn<fp_t(fp_t,
fermion &,fp_t)>
830 this,std::placeholders::_1,std::ref(f),temper);
831 funct mfs=std::bind(std::mem_fn<fp_t(fp_t,
fermion &,fp_t)>
833 this,std::placeholders::_1,std::ref(f),temper);
835 f.ed=
nit->integ(mfe,0.0,0.0);
837 if (!f.inc_rest_mass) f.ed-=f.n*f.m;
840 f.en=
nit->integ(mfs,0.0,0.0);
847 funct mfe=std::bind(std::mem_fn<fp_t(fp_t,
fermion &,fp_t)>
849 this,std::placeholders::_1,std::ref(f),temper);
850 funct mfs=std::bind(std::mem_fn<fp_t(fp_t,
fermion &,fp_t)>
852 this,std::placeholders::_1,std::ref(f),temper);
855 if (f.inc_rest_mass) {
866 if (f.inc_rest_mass) {
882 f.ed=
dit->integ(mfe,0.0,ul);
887 f.en=
dit->integ(mfs,ll,ul);
890 f.en=
dit->integ(mfs,0.0,ul);
902 O2SCL_ERR2(
"Zero density in degenerate limit in fermion_rel::",
903 "calc_mu(). Variable deg_limit set improperly?",
910 f.pr=-f.ed+temper*f.en+f.nu*f.n;
920 template<
class fermion_t>
925 if (f.non_interacting) { f.nu=f.mu; f.ms=f.m; }
959 if (f.inc_rest_mass) {
962 f.ed=f.ed+antip.
ed+2.0*antip.
n*f.m;
980 template<
class fermion_t>
995 fp_t density_temp=f.n;
997 if (f.non_interacting==
true) { f.nu=f.mu; f.ms=f.m; }
999 fp_t initial_guess=f.nu;
1001 fp_t nex=f.nu/temper;
1003 funct mf=std::bind(std::mem_fn<fp_t(fp_t,
fermion &,fp_t,
bool)>
1005 this,std::placeholders::_1,std::ref(f),temper,
false);
1016 fp_t lg=std::max(fabs(f.nu),f.ms);
1017 fp_t bhigh=lg/temper, blow=-bhigh;
1018 fp_t yhigh=mf(bhigh), ylow=mf(blow);
1019 for(
size_t j=0;j<5 && yhigh<0.0;j++) {
1023 for(
size_t j=0;j<5 && ylow>0.0;j++) {
1027 if (yhigh>0.0 && ylow<0.0) {
1041 fp_t tol1=
dit->tol_rel, tol2=
dit->tol_abs;
1042 fp_t tol3=
nit->tol_rel, tol4=
nit->tol_abs;
1043 dit->tol_rel/=1.0e2;
1044 dit->tol_abs/=1.0e2;
1045 nit->tol_rel/=1.0e2;
1046 nit->tol_abs/=1.0e2;
1052 if (nex<0.0) nex=1.0e-10;
1060 this,std::placeholders::_1,std::ref(f),
1075 ret=rbg.
solve(nex,lmf);
1091 std::cout.precision(14);
1092 std::cout <<
"m,ms,n,T: " << f.m <<
" " << f.ms <<
" "
1093 << f.n <<
" " << temper << std::endl;
1094 std::cout <<
"nu: " << initial_guess << std::endl;
1101 if (f.non_interacting==
true) { f.mu=f.nu; }
1120 #ifndef DOXYGEN_INTERNAL
1135 ret=(eta+u)*sqrt(u*u+2.0*eta*u);
1137 ret=(eta+u)*sqrt(u*u+2.0*eta*u)/(exp(eta+u-y)+1.0);
1139 ret=(eta+u)*sqrt(u*u+2.0*eta*u)*exp(y)/(exp(eta+u)+exp(y));
1142 if (!std::isfinite(ret)) {
1162 ret=(eta+u)*(eta+u)*sqrt(u*u+2.0*eta*u)/(exp(eta+u-y)+1.0);
1164 ret=(eta+u)*(eta+u)*sqrt(u*u+2.0*eta*u)*exp(y)/(exp(eta+u)+exp(y));
1167 if (!std::isfinite(ret)) {
1176 double ret, y, eta, term1, term2;
1185 term1=log(1.0+exp(y-eta-u))/(1.0+exp(y-eta-u));
1186 term2=log(1.0+exp(eta+u-y))/(1.0+exp(eta+u-y));
1187 ret=(eta+u)*sqrt(u*u+2.0*eta*u)*(term1+term2);
1189 if (!std::isfinite(ret)) {
1200 double E=gsl_hypot(k,f.
ms), ret;
1203 ret=k*k/(1.0+exp((E-f.
nu)/T));
1205 if (!std::isfinite(ret)) {
1207 "in fermion_rel::deg_density_fun().",
exc_einval);
1216 double E=gsl_hypot(k,f.
ms), ret;
1219 ret=k*k*E/(1.0+exp((E-f.
nu)/T));
1221 if (!std::isfinite(ret)) {
1223 "in fermion_rel::deg_energy_fun().",
exc_einval);
1232 double E=gsl_hypot(k,f.
ms), ret;
1243 double arg=E/T-f.
nu/T;
1244 ret=-k*k*(-1.0+arg)*exp(arg);
1246 double nx=1.0/(1.0+exp(E/T-f.
nu/T));
1247 ret=-k*k*(nx*log(nx)+(1.0-nx)*log(1.0-nx));
1250 if (!std::isfinite(ret)) {
1252 "in fermion_rel::deg_entropy_fun().",
exc_einval);
1281 yy=(ntemp-f.
n)/ntemp;
1294 yy=(ntemp-f.
n)/ntemp;
1304 funct mfe=std::bind(std::mem_fn<
double(
double,
fermion &,
double)>
1306 this,std::placeholders::_1,std::ref(f),T);
1308 nden=
nit->integ(mfe,0.0,0.0);
1316 funct mfe=std::bind(std::mem_fn<
double(
double,
fermion &,
double)>
1318 this,std::placeholders::_1,std::ref(f),T);
1333 nden=
dit->integ(mfe,0.0,ul);
1362 double nn_match=f.
n;
1364 double nden_p, nden_ap;
1375 if (!std::isfinite(f.
nu))
return 3;
1389 if (this->
calc_mu_ndeg(f,T,1.0e-8,
true) && std::isfinite(f.
n)) {
1390 double y1=f.
n/nn_match-1.0;
1391 if (!std::isfinite(y1)) {
1392 O2SCL_ERR(
"Value 'y1' not finite (10) in fermion_rel::pair_fun().",
1415 bool particles_done=
false;
1420 particles_done=
true;
1422 if (!std::isfinite(nden_p)) {
1423 O2SCL_ERR(
"Value 'nden_p' not finite (1) in fermion_rel::pair_fun().",
1431 if (this->
calc_mu_deg(f,T,1.0e-8) && std::isfinite(f.
n)) {
1432 particles_done=
true;
1434 if (!std::isfinite(nden_p)) {
1435 O2SCL_ERR(
"Value 'nden_p' not finite (2) in fermion_rel::pair_fun().",
1442 if (particles_done==
false) {
1448 funct mfe=std::bind(std::mem_fn<
double(
double,
fermion &,
double)>
1450 this,std::placeholders::_1,std::ref(f),T);
1452 nden_p=
nit->integ(mfe,0.0,0.0);
1454 if (!std::isfinite(nden_p)) {
1455 O2SCL_ERR(
"Value 'nden_p' not finite (3) in fermion_rel::pair_fun().",
1463 funct mfe=std::bind(std::mem_fn<
double(
double,
fermion &,
double)>
1465 this,std::placeholders::_1,std::ref(f),T);
1477 nden_p=
dit->integ(mfe,0.0,ul);
1483 if (!std::isfinite(nden_p)) {
1484 O2SCL_ERR(
"Value 'nden_p' not finite (4) in fermion_rel::pair_fun().",
1490 particles_done=
true;
1500 if (log_mode) f.
nu=-T*exp(x);
1503 if (log_mode) f.
nu=-T*exp(x)-2.0*f.
m;
1507 bool antiparticles_done=
false;
1521 antiparticles_done=
true;
1523 if (!std::isfinite(nden_ap)) {
1524 O2SCL_ERR2(
"Value 'nden_ap' not finite (5) in",
1525 "fermion_rel::pair_fun().",
1534 antiparticles_done=
true;
1536 if (!std::isfinite(nden_ap)) {
1537 O2SCL_ERR2(
"Value 'nden_ap' not finite (6) in",
1538 "fermion_rel::pair_fun().",
1545 if (antiparticles_done==
false) {
1551 funct mf=std::bind(std::mem_fn<
double(
double,
fermion &,
double)>
1553 this,std::placeholders::_1,std::ref(f),T);
1555 nden_ap=
nit->integ(mf,0.0,0.0);
1557 if (!std::isfinite(nden_ap)) {
1558 O2SCL_ERR2(
"Value 'nden_ap' not finite (7) in",
1559 "fermion_rel::pair_fun().",
1567 funct mf=std::bind(std::mem_fn<
double(
double,
fermion &,
double)>
1569 this,std::placeholders::_1,std::ref(f),T);
1581 nden_ap=
dit->integ(mf,0.0,ul);
1586 if (!std::isfinite(nden_ap)) {
1587 O2SCL_ERR2(
"Value 'nden_ap' not finite (8) in",
1588 "fermion_rel::pair_fun().",
1594 antiparticles_done=
true;
1599 if (nn_match==0.0) {
1600 y2=fabs(nden_p-nden_ap)/fabs(nden_p);
1602 y2=(nden_p-nden_ap)/nn_match-1.0;
1605 if (!std::isfinite(y2)) {
1606 O2SCL_ERR(
"Value 'y2' not finite (9) in fermion_rel::pair_fun().",
1626 #ifndef DOXYGEN_NO_O2NS