26 #ifndef O2SCL_SKYRME_EOS_H
27 #define O2SCL_SKYRME_EOS_H
33 #include <o2scl/constants.h>
34 #include <o2scl/mroot.h>
35 #include <o2scl/eos_had_base.h>
36 #include <o2scl/part.h>
37 #include <o2scl/fermion_nonrel.h>
38 #include <o2scl/fermion_deriv_nr.h>
40 #ifndef DOXYGEN_NO_O2NS
236 double &ham3,
double &ham4,
237 double &ham5,
double &ham6);
244 template<
class fermion_t>
246 (fermion_t &ne, fermion_t &pr,
double ltemper,
thermo &locth,
247 double term,
double term2,
double ham1,
double ham2,
248 double ham3,
double ham4,
double ham5,
double ham6) {
251 double na=pow(nb,alpha);
252 double npa=pow(pr.n,alpha);
253 double nna=pow(ne.n,alpha);
255 double ham=ne.ed+pr.ed+ham1*nb*nb+ham2*(ne.n*ne.n+pr.n*pr.n)+
256 ham3*na*ne.n*pr.n+ham4*(nna*ne.n*ne.n+npa*pr.n*pr.n)+
257 ham5*nb*nb*na+ham6*(ne.n*ne.n+pr.n*pr.n)*na;
260 if (ne.inc_rest_mass) {
261 gn=2.0*ne.ms*(ne.ed-ne.n*ne.m);
265 if (pr.inc_rest_mass) {
266 gp=2.0*pr.ms*(pr.ed-pr.n*pr.m);
273 double common=2.0*ham1*nb+ham5*(2.0+alpha)*nb*na;
274 double dhdnn=common+2.0*ham2*ne.n+ham3*na*pr.n*(alpha*ne.n/nb+1.0)+
275 ham4*(nna*ne.n*(2.0+alpha))+
276 ham6*(2.0*ne.n*na+(ne.n*ne.n+pr.n*pr.n)*alpha*na/nb);
277 double dhdnp=common+2.0*ham2*pr.n+ham3*na*ne.n*(alpha*pr.n/nb+1.0)+
278 ham4*(npa*pr.n*(2.0+alpha))+
279 ham6*(2.0*pr.n*na+(ne.n*ne.n+pr.n*pr.n)*alpha*na/nb);
282 ne.mu=ne.nu+dhdnn+(gn+gp)*term+gn*term2;
283 pr.mu=pr.nu+dhdnp+(gn+gp)*term+gp*term2;
287 locth.
en=ne.en+pr.en;
288 locth.
pr=ltemper*locth.
en+ne.mu*ne.n+pr.mu*pr.n-locth.
ed;
295 template<
class fermion_t>
297 (fermion_t &ne, fermion_t &pr,
double ltemper,
thermo &locth,
298 thermo_np_deriv_helm &locthd,
double term,
double term2,
299 double ham1,
double ham2,
double ham3,
300 double ham4,
double ham5,
double ham6) {
303 double na=pow(nb,alpha);
304 double npa=pow(pr.n,alpha);
305 double nna=pow(ne.n,alpha);
309 double opatpa=(1.0+alpha)*(2.0+alpha);
310 double common2=2.0*ham1+2.0*ham2;
311 double dhdnn2=common2+ham4*nna*opatpa+ham5*na*opatpa+
312 ham3*(2.0*pr.n/nb*na*alpha+
313 ne.n*pr.n*na/nb/nb*(-1.0+alpha)*alpha)+
314 ham6*(2.0*na+4.0*ne.n/nb*na*alpha+na/nb/nb*
315 (ne.n*ne.n+pr.n*pr.n)*(-1.0+alpha)*alpha);
316 double dhdnp2=common2+ham4*npa*opatpa+ham5*na*opatpa+
317 ham3*(2.0*ne.n/nb*na*alpha+
318 ne.n*pr.n*na/nb/nb*(-1.0+alpha)*alpha)+
319 ham6*(2.0*na+4.0*pr.n/nb*na*alpha+na/nb/nb*
320 (ne.n*ne.n+pr.n*pr.n)*(-1.0+alpha)*alpha);
321 double dhdnndnp=2.0*ham1+na/nb/nb*
322 (ham5*nb*nb*opatpa+ham6*alpha*
323 (4.0*ne.n*pr.n+ne.n*ne.n*(1.0+alpha)+pr.n*pr.n*(1.0+alpha))+
324 ham3*(ne.n*ne.n*(1.0+alpha)+pr.n*pr.n*(1.0+alpha)+
325 ne.n*pr.n*(2.0+alpha+alpha*alpha)));
330 double n_dsdT_f=0.0, p_dsdT_f=0.0;
331 double n_dmudT_f=0.0, p_dmudT_f=0.0;
332 double n_dmudn_f=0.0, p_dmudn_f=0.0;
333 ne.deriv_f(n_dmudn_f,n_dmudT_f,n_dsdT_f);
334 pr.deriv_f(p_dmudn_f,p_dmudT_f,p_dsdT_f);
339 if (ne.inc_rest_mass) {
340 hn=10.0*ne.ms*(ne.ed-ne.n*ne.m)*ne.ms;
342 hn=10.0*ne.ms*ne.ed*ne.ms;
344 if (pr.inc_rest_mass) {
345 hp=10.0*pr.ms*(pr.ed-pr.n*pr.m)*pr.ms;
347 hp=10.0*pr.ms*pr.ed*pr.ms;
352 locthd.dsdT=n_dsdT_f+p_dsdT_f;
354 locthd.dmundT=2.0*ltemper*ne.ms*(term+term2)*n_dsdT_f+
355 2.0*ltemper*pr.ms*term*p_dsdT_f+n_dmudT_f;
357 locthd.dmupdT=2.0*ltemper*pr.ms*(term+term2)*p_dsdT_f+
358 2.0*ltemper*ne.ms*term*n_dsdT_f+p_dmudT_f;
360 locthd.dmundnn=-(term+term2)*(term+term2)*hn-term*term*hp+
361 pow(1.0+3.0*(term+term2)*ne.n*ne.ms,2.0)*n_dmudn_f+
362 9.0*term*term*pr.n*pr.ms*pr.n*pr.ms*p_dmudn_f+dhdnn2;
364 locthd.dmupdnp=-(term+term2)*(term+term2)*hp-term*term*hn+
365 pow(1.0+3.0*(term+term2)*pr.n*pr.ms,2.0)*p_dmudn_f+
366 9.0*term*term*ne.n*ne.ms*ne.n*ne.ms*n_dmudn_f+dhdnp2;
368 locthd.dmudn_mixed=-term*(term+term2)*(hn+hp)+
369 3.0*term*ne.ms*ne.n*(1.0+3.0*(term+term2)*ne.n*ne.ms)*n_dmudn_f+
370 3.0*term*pr.ms*pr.n*(1.0+3.0*(term+term2)*pr.n*pr.ms)*p_dmudn_f+
386 template<
class fermion_t>
389 if (!std::isfinite(ne.n) || !std::isfinite(pr.n) ||
391 O2SCL_ERR2(
"Nucleon densities or temperature not finite in ",
394 if (ne.n<0.0 || pr.n<0.0) {
395 std::string str=((std::string)
"Nucleon densities negative, n_n=")+
396 std::to_string(ne.n)+
", n_p="+std::to_string(pr.n)+
", in "+
397 "eos_had_skyrme::check_input().";
400 if (fabs(ne.g-2.0)>1.0e-10 || fabs(pr.g-2.0)>1.0e-10) {
401 O2SCL_ERR((((std::string)
"Neutron (")+std::to_string(ne.g)+
402 ") or proton ("+std::to_string(pr.g)+
") spin deg"+
403 "eneracies wrong in "+
404 "eos_had_skyrme::check_input().").c_str(),
407 if (fabs(ne.m-4.5)>1.0 || fabs(pr.m-4.5)>1.0) {
408 O2SCL_ERR((((std::string)
"Neutron (")+std::to_string(ne.m)+
409 ") or proton ("+std::to_string(pr.m)+
") masses wrong "+
410 "in eos_had_skyrme::check_input().").c_str(),
413 if (ne.non_interacting==
true || pr.non_interacting==
true) {
414 O2SCL_ERR2(
"Neutron or protons non-interacting in ",
443 virtual int calc_temp_e(fermion &ne, fermion &pr,
double temper,
450 double temper,
thermo &th,
451 thermo_np_deriv_helm &thd);
456 virtual int calc_e(fermion &ne, fermion &pr,
thermo <);
459 virtual const char *
type() {
return "eos_had_skyrme"; }
464 double t0, t1, t2, t3, x0, x1, x2, x3, alpha, a, b;
508 virtual double feoa(
double nb);
516 virtual double fmsom(
double nb);
544 virtual double fesym(
double nb,
double alpha=0.0);
575 double &f0,
double &g0,
double &f0p,
576 double &g0p,
double &f1,
double &g1,
577 double &f1p,
double &g1p);
589 double &f1,
double &g1);
598 template<
class fermion_t>
599 void eff_mass(fermion_t &ne, fermion_t &pr,
double &term,
604 term=0.25*(t1*(1.0+x1/2.0)+t2*(1.0+x2/2.0));
605 term2=0.25*(t2*(0.5+x2)-t1*(0.5+x1));
606 ne.ms=ne.m/(1.0+2.0*(nb*term+ne.n*term2)*ne.m);
607 pr.ms=pr.m/(1.0+2.0*(nb*term+pr.n*term2)*pr.m);
641 int calpar(
double gt0=-10.0,
double gt3=70.0,
double galpha=0.2,
642 double gt1=2.0,
double gt2=-1.0);
690 (
double Crr00,
double Crr10,
double Crr0D,
double Crr1D,
691 double Crt0,
double Crt1,
double CrDr0,
double CrDr1,
692 double CrnJ0,
double CrnJ1,
double alpha2);
712 (
double &Crr00,
double &Crr10,
double &Crr0D,
double &Crr1D,
713 double &Crt0,
double &Crt1,
double &CrDr0,
double &CrDr1,
714 double &CrnJ0,
double &CrnJ1,
double &alpha2);
743 (
double n0,
double EoA,
double K,
double Ms_star,
double a,
double L,
744 double Mv_star,
double CrDr0,
double CrDr1,
double CrnJ0,
double CrnJ1);
756 #ifndef DOXYGEN_NO_O2NS
764 double fixn0, fixeoa, fixesym, fixcomp, fixmsom;
771 #ifndef DOXYGEN_NO_O2NS