fermion_rel.h
Go to the documentation of this file.
1 /* -------------------------------------------------------------------
2 
3  Copyright (C) 2006-2020, Andrew W. Steiner
4 
5  This file is part of O2scl.
6 
7  O2scl is free software; you can redistribute it and/or modify
8  it under the terms of the GNU General Public License as published by
9  the Free Software Foundation; either version 3 of the License, or
10  (at your option) any later version.
11 
12  O2scl is distributed in the hope that it will be useful,
13  but WITHOUT ANY WARRANTY; without even the implied warranty of
14  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  GNU General Public License for more details.
16 
17  You should have received a copy of the GNU General Public License
18  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
19 
20  -------------------------------------------------------------------
21 */
22 #ifndef O2SCL_FERMION_REL_H
23 #define O2SCL_FERMION_REL_H
24 
25 /** \file fermion_rel.h
26  \brief File defining \ref o2scl::fermion_rel_tl
27 */
28 
29 #include <string>
30 #include <iostream>
31 #include <fstream>
32 #include <cmath>
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>
40 
41 #ifndef DOXYGEN_NO_O2NS
42 namespace o2scl {
43 #endif
44 
45  /** \brief Equation of state for a relativistic fermion
46 
47  This class computes the thermodynamics of a relativistic fermion
48  either as a function of the density or the chemical potential.
49  It employs direct integration, using two different integrators
50  for the degenerate and non-degenerate regimes. The default
51  integrators are inte_qag_gsl (for degenerate fermions) and
52  inte_qagiu_gsl (for non-degenerate fermions). For the functions
53  calc_mu() and calc_density(), if the temperature argument is
54  less than or equal to zero, the functions \ref
55  fermion_zerot::calc_mu_zerot() and \ref
56  fermion_zerot::calc_density_zerot() will be used to compute the
57  result.
58 
59  \hline
60  <b>Degeneracy parameter:</b>
61 
62  Define the degeneracy parameter
63  \f[
64  \psi=(\nu-m^{*})/T
65  \f]
66  where \f$ \nu \f$ is the effective chemical potential (including
67  the rest mass) and \f$
68  m^{*} \f$ is the effective mass. For \f$ \psi \f$ smaller than
69  \ref min_psi, the non-degenerate expansion in \ref
70  fermion_thermo::calc_mu_ndeg() is attempted first. If that
71  fails, then integration is used. For \f$ \psi \f$ greater than
72  \ref deg_limit (degenerate regime), a finite interval integrator
73  is used and for \f$ \psi \f$ less than \ref deg_limit
74  (non-degenerate regime), an integrator over the interval from
75  \f$ [0,\infty) \f$ is used. In the case where \ref
76  part::inc_rest_mass is false, the degeneracy parameter is
77  \f[
78  \psi=(\nu+m-m^{*})/T
79  \f]
80 
81  <b>Integration limits:</b>
82 
83  The upper limit on the degenerate integration is given by
84  \f[
85  \mathrm{upper~limit} = \sqrt{{\cal L}^2-m^{*,2}}
86  \f]
87  where \f$ {\cal L}\equiv u T+\nu \f$ and \f$ u \f$ is \ref
88  fermion_rel::upper_limit_fac . In the case where \ref
89  part::inc_rest_mass is false, the result is
90  \f[
91  \mathrm{upper~limit} = \sqrt{(m+{\cal L})^2-m^{*2}}
92  \f]
93 
94  The entropy is only significant at the Fermi surface, thus
95  in the degenerate case, the lower limit of the entropy
96  integral can be given be determined by the value of \f$ k \f$
97  which solves
98  \f[
99  - u = \frac{\sqrt{k^2+m^{* 2}}-\nu}{T}
100  \f]
101  The solution is
102  \f[
103  \mathrm{lower~limit} = \sqrt{(-u T+{\nu})^2-m^{*,2}}
104  \f]
105  but this solution is only valid if \f$ (m^{*}-\nu)/T < -u \f$.
106  In the case where part::inc_rest_mass is false, the result is
107  \f[
108  \mathrm{lower~limit} = \sqrt{(-u T + m +\nu)^2-m^{*,2}}
109  \f]
110  which is valid if \f$ (m^{*}-\nu - m)/T < -u \f$.
111 
112  <b>Entropy integrand:</b>
113 
114  In the degenerate regime, the entropy integrand
115  \f[
116  - k^2 \left[ f \log f + \left(1-f\right) \log
117  \left(1-f \right) \right]
118  \f]
119  where \f$ f \f$ is the fermionic distribution function can lose
120  precision when \f$ (E^{*} - \nu)/T \f$ is negative and
121  sufficiently large in absolute magnitude. Thus when \f$ (E^{*} -
122  \nu)/T < S \f$ where \f$ S \f$ is stored in \ref deg_entropy_fac
123  (default is -30), the integrand is written as
124  \f[
125  -k^2 \left( E/T-\nu/T \right) e^{E/T-\nu/T} \, .
126  \f]
127  If \f$ (E - \nu)/T < S \f$ is less than -1 times \ref exp_limit
128  (e.g. less than -200), then the entropy integrand is assumed
129  to be zero.
130 
131  <b>Non-degenerate integrands:</b>
132 
133  \comment
134  It's not at all clear that this dimensionless form is more
135  accurate than other potential alternatives. On the other hand,
136  it seems that the uncertainties in the integrations are larger
137  than the errors made by the integrand at present.
138  \endcomment
139  The integrands in the non-degenerate regime are written
140  in a dimensionless form, by defining \f$ u \f$ with
141  the relation
142  \f$ p = \sqrt{\left(T u + m^{*}\right)^2-m^{* 2}} \f$,
143  \f$ y \equiv \nu/ T \f$, and
144  \f$ \mathrm{mx} \equiv m^{*}/T \f$.
145  The density integrand is
146  \f[
147  \left(\mathrm{mx}+u\right) \sqrt{u^2+2 (\mathrm{mx}) u}
148  \left(\frac{e^{y}}{e^{\mathrm{mx}+u}+e^{y}}\right) \, ,
149  \f]
150  the energy integrand is
151  \f[
152  \left(\mathrm{mx}+u\right)^2 \sqrt{u^2+2 (\mathrm{mx}) u}
153  \left(\frac{e^{y}}{e^{\mathrm{mx}+u}+e^{y}}\right) \, ,
154  \f]
155  and the entropy integrand is
156  \f[
157  \left(\mathrm{mx}+u\right) \sqrt{u^2+2 (\mathrm{mx}) u}
158  \left(t_1+t_2\right) \, ,
159  \f]
160  where
161  \f{eqnarray*}
162  t_1 &=& \log \left(1+e^{y-\mathrm{mx}-u}\right)/
163  \left(1+e^{y-\mathrm{mx}-u}\right) \nonumber \\
164  t_2 &=& \log \left(1+e^{\mathrm{mx}+u-y}\right)/
165  \left(1+e^{\mathrm{mx}+u-y}\right) \, .
166  \f}
167 
168  \hline
169  <b>Accuracy:</b>
170 
171  The default settings for for this class give an accuracy of at
172  least 1 part in \f$ 10^6 \f$ (and frequently better than this).
173 
174  When the integrators provide numerical uncertainties, these
175  uncertainties are stored in \ref unc. In the case of
176  calc_density() and pair_density(), the uncertainty from the
177  numerical accuracy of the solver is not included. (There is also
178  a relatively small inaccuracy due to the mathematical evaluation
179  of the integrands which is not included in \ref unc.)
180 
181  One can improve the accuracy to within 1 part in \f$ 10^{10} \f$
182  using
183  \code
184  fermion_rel rf(1.0,2.0);
185  rf.upper_limit_fac=40.0;
186  rf.dit->tol_abs=1.0e-13;
187  rf.dit->tol_rel=1.0e-13;
188  rf.nit->tol_abs=1.0e-13;
189  rf.nit->tol_rel=1.0e-13;
190  rf.density_root->tol_rel=1.0e-10;
191  \endcode
192  which decreases the both the relative and absolute tolerances
193  for both the degenerate and non-degenerate integrators and
194  improves the accuracy of the solver which determines the
195  chemical potential from the density. Of course if these
196  tolerances are too small, the calculation may fail.
197 
198  \hline
199  <b>Todos:</b>
200 
201  \future The expressions which appear in in the integrand
202  functions density_fun(), etc. could likely be improved,
203  especially in the case where \ref o2scl::part::inc_rest_mass is
204  <tt>false</tt>. There should not be a need to check if
205  <tt>ret</tt> is finite.
206 
207  \future It appears this class doesn't compute the uncertainty in
208  the chemical potential or density with calc_density(). This
209  could be fixed.
210 
211  \future I'd like to change the lower limit on the entropy
212  integration, but the value in the code at the moment (stored
213  in <tt>ll</tt>) makes bm_part2.cpp worse.
214 
215  \future The function pair_mu() should set the antiparticle
216  integrators as done in fermion_deriv_rel.
217  */
218  template<class fp_t=double>
219  class fermion_rel_tl : public fermion_thermo_tl<fp_t> {
220 
221  public:
222 
223  /// \name Numerical parameters
224  //@{
225  /** \brief If true, call the error handler when convergence
226  fails (default true)
227  */
229 
230  /** \brief The smallest value of \f$ (\mu-m)/T \f$ for which
231  integration is used
232  */
233  fp_t min_psi;
234 
235  /** \brief The critical degeneracy at which to switch integration
236  techniques (default 2)
237  */
238  fp_t deg_limit;
239 
240  /** \brief The limit for exponentials to ensure integrals are finite
241  (default 200)
242  */
243  fp_t exp_limit;
244 
245  /// The factor for the degenerate upper limits (default 20)
247 
248  /// A factor for the degenerate entropy integration (default 30)
250 
251  /// Verbosity parameter (default 0)
252  int verbose;
253  //@}
254 
255  /// Storage for the uncertainty
257 
258  /// If true, use expansions for extreme conditions (default true)
260 
261  /// Create a fermion with mass \c m and degeneracy \c g
263  dit(new inte_qag_gsl<>),
264  density_root(new root_cern<>) {
265  deg_limit=2.0;
266 
267  exp_limit=200.0;
268  upper_limit_fac=20.0;
269  deg_entropy_fac=30.0;
270  min_psi=-4.0;
271  err_nonconv=true;
272  use_expansions=true;
273  density_root->tol_rel=4.0e-7;
274  verbose=0;
275  last_method=0;
276  }
277 
278  virtual ~fermion_rel_tl() {
279  }
280 
281  /** \brief Calculate properties as function of chemical potential
282  */
283  virtual void calc_mu(fermion &f, fp_t temper) {
284  calc_mu_tlate<fermion>(f,temper);
285  return;
286  }
287 
288  /** \brief Calculate properties as function of density
289 
290  This function uses the current value of \c nu (or \c mu if the
291  particle is non interacting) for an initial guess to solve for
292  the chemical potential. If this guess is too small, then this
293  function may fail.
294  */
295  virtual int calc_density(fermion &f, fp_t temper) {
296  return calc_density_tlate<fermion>(f,temper);
297  }
298 
299  /** \brief Calculate properties with antiparticles as function of
300  chemical potential
301  */
302  virtual void pair_mu(fermion &f, fp_t temper) {
303  pair_mu_tlate<fermion>(f,temper);
304  return;
305  }
306 
307  /** \brief Calculate properties with antiparticles as function of
308  density
309  */
310  virtual int pair_density(fermion &f, fp_t temper) {
311  return pair_density_tlate<fermion>(f,temper);
312  }
313 
314  /** \brief Calculate effective chemical potential from density
315  */
316  virtual int nu_from_n(fermion &f, fp_t temper) {
317  return nu_from_n_tlate<fermion>(f,temper);
318  }
319 
320  /// The non-degenerate integrator
321  std::shared_ptr<inte<> > nit;
322 
323  /// The degenerate integrator
324  std::shared_ptr<inte<> > dit;
325 
326  /// The solver for calc_density()
327  std::shared_ptr<root<> > density_root;
328 
329  /// Return string denoting type ("fermion_rel")
330  virtual const char *type() { return "fermion_rel"; }
331 
332  /** \brief An integer indicating the last numerical method used
333 
334  In all functions
335  - 0: no previous calculation or last calculation failed
336 
337  In \ref nu_from_n_tlate():
338  - 1: default solver
339  - 2: default solver with smaller tolerances
340  - 3: bracketing solver
341 
342  In \ref calc_mu_tlate():
343  - 4: non-degenerate expansion
344  - 5: degenerate expansion
345  - 6: exact integration, non-degenerate integrands
346  - 7: exact integration, degenerate integrands, lower limit
347  on entropy integration
348  - 8: exact integration, degenerate integrands, full
349  entropy integration
350  - 9: T=0 result
351 
352  In \ref calc_density_tlate(), the integer is a two-digit
353  number. The first digit (1 to 3) is the method used by \ref
354  nu_from_n_tlate() and the second digit is one of
355  - 1: nondegenerate expansion
356  - 2: degenerate expansion
357  - 3: exact integration, non-degenerate integrands
358  - 4: exact integration, degenerate integrands, lower limit
359  on entropy integration
360  - 5: exact integration, degenerate integrands, full
361  entropy integration
362  If \ref calc_density_tlate() uses the T=0 code, then
363  last_method is 40.
364 
365  In \ref pair_mu_tlate(), the integer is a three-digit number.
366  The third digit is always 0 (to ensure a value of last_method
367  which is unique from the other values reported from other
368  functions as described above). The first digit is the method
369  used for particles from \ref calc_mu_tlate() above and the
370  second digit is the method used for antiparticles.
371 
372  In \ref pair_density_tlate(), the integer is a four-digit
373  number. The first digit is from the list below and the
374  remaining three digits, if nonzero, are from \ref
375  pair_mu_tlate().
376  - 1: T=0 result
377  - 2: default solver
378  - 3: bracketing solver
379  - 4: default solver with smaller tolerances
380  - 5: default solver with smaller tolerances in log units
381  - 6: bracketing in log units
382  */
384 
385  /// \name Template versions of base functions
386  //@{
387  /** \brief Calculate the chemical potential from the density
388  (template version)
389  */
390  template<class fermion_t>
391  int nu_from_n_tlate(fermion_t &f, fp_t temper) {
392 
393  last_method=0;
394 
395  fp_t nex;
396 
397  // Try to ensure a good initial guess
398 
399  nex=f.nu/temper;
400  fp_t y=solve_fun(nex,f,temper);
401  if (verbose>1) {
402  std::cout << "nu_from_n(): initial guess " << nex << std::endl;
403  }
404 
405  if (y>1.0-1.0e-6) {
406  fp_t scale=f.ms;
407  if (temper>scale) scale=temper;
408  for(size_t i=0;i<10;i++) {
409  if (nex<0.0) nex+=scale*1.0e5;
410  else nex*=10.0;
411  y=solve_fun(nex,f,temper);
412  if (y<1.0-1.0e-6) i=10;
413  }
414  if (verbose>1) {
415  std::cout << "nu_from_n(): adjusted guess to " << nex << std::endl;
416  }
417  }
418 
419  // If that didn't work, try a different guess
420  if (y>1.0-1.0e-6) {
421  if (f.inc_rest_mass) {
422  nex=f.ms/temper;
423  } else {
424  nex=(f.ms-f.m)/temper;
425  }
426  y=solve_fun(nex,f,temper);
427  if (verbose>1) {
428  std::cout << "nu_from_n(): adjusted guess (try 2) to "
429  << nex << std::endl;
430  }
431  }
432 
433  // If neither worked, call the error handler
434  if (y==1.0 || !std::isfinite(y)) {
435  O2SCL_CONV2_RET("Couldn't find reasonable initial guess in ",
436  "fermion_rel::nu_from_n().",exc_einval,
437  this->err_nonconv);
438  }
439 
440  // Perform full solution
441  funct mf=std::bind(std::mem_fn<fp_t(fp_t,fermion &,fp_t)>
443  this,std::placeholders::_1,std::ref(f),temper);
444 
445  // The default o2scl::root object is of type root_cern,
446  // and this solver has problems when the root is near 0.
447  // Below, we switch to a root_brent_gsl object in the case
448  // that the default solver fails.
449 
450  bool drec=density_root->err_nonconv;
451  density_root->err_nonconv=false;
452  int ret=density_root->solve(nex,mf);
453 
454  if (ret!=0) {
455 
456  if (verbose>1) {
457  std::cout << "nu_from_n(): density_root failed x="
458  << nex << " ." << std::endl;
459  std::cout << "\tTrying to make integrators more accurate."
460  << std::endl;
461  }
462 
463  // If it fails, try to make the 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;
466  dit->tol_rel/=1.0e2;
467  dit->tol_abs/=1.0e2;
468  nit->tol_rel/=1.0e2;
469  nit->tol_abs/=1.0e2;
470  ret=density_root->solve(nex,mf);
471 
472  if (ret!=0) {
473 
474  if (verbose>1) {
475  std::cout << "nu_from_n(): density_root failed again x=" << nex
476  << " ." << std::endl;
477  std::cout << "Trying to bracket root." << std::endl;
478  }
479 
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++) {
484  bhigh*=1.0e2;
485  yhigh=mf(bhigh);
486  }
487  for(size_t j=0;j<5 && ylow<0.0;j++) {
488  blow*=1.0e2;
489  ylow=mf(blow);
490  }
491  if (yhigh<0.0 && ylow>0.0) {
493  rbg.err_nonconv=false;
494  ret=rbg.solve_bkt(blow,bhigh,mf);
495  if (ret==0) {
496  // Bracketing solver worked
497  last_method=3;
498  nex=blow;
499  } else {
500  if (verbose>1) {
501  std::cout << "nu_from_n(): density_root failed fourth solver "
502  << blow << std::endl;
503  }
504  }
505  } else if (verbose>1) {
506  std::cout << "nu_from_n(): Failed to bracket." << std::endl;
507  }
508  } else {
509  // Increasing tolerances worked
510  last_method=2;
511  }
512 
513  // Return tolerances to their original values
514  dit->tol_rel=tol1;
515  dit->tol_abs=tol2;
516  nit->tol_rel=tol3;
517  nit->tol_abs=tol4;
518 
519  } else {
520  // First solver worked
521  last_method=1;
522  }
523 
524  density_root->err_nonconv=drec;
525 
526  if (ret!=0) {
527  O2SCL_CONV2_RET("Density solver failed in ",
528  "fermion_rel::nu_from_n().",exc_efailed,
529  this->err_nonconv);
530  }
531 
532  f.nu=nex*temper;
533 
534  return success;
535  }
536 
537  /** \brief Calculate properties as function of chemical potential
538  (template version)
539  */
540  template<class fermion_t>
541  void calc_mu_tlate(fermion_t &f, fp_t temper) {
542 
543  last_method=0;
544 
545  // -----------------------------------------------------------------
546  // Handle T<=0
547 
548  if (temper<0.0) {
549  O2SCL_ERR2("Temperature less than zero in ",
550  "fermion_rel::calc_mu_tlate().",exc_einval);
551  }
552  if (temper==0.0) {
553  this->calc_mu_zerot(f);
554  last_method=9;
555  return;
556  }
557 
558  if (f.non_interacting==true) { f.nu=f.mu; f.ms=f.m; }
559 
560  // Compute the degeneracy parameter
561 
562  bool deg=true;
563  fp_t psi;
564  if (f.inc_rest_mass) {
565  psi=(f.nu-f.ms)/temper;
566  } else {
567  psi=(f.nu+(f.m-f.ms))/temper;
568  }
569  if (psi<deg_limit) deg=false;
570 
571  // Try the non-degenerate expansion if psi is small enough
572  if (use_expansions && psi<min_psi) {
573  bool acc=this->calc_mu_ndeg(f,temper,1.0e-14);
574  if (acc) {
575  unc.n=f.n*1.0e-14;
576  unc.ed=f.ed*1.0e-14;
577  unc.pr=f.pr*1.0e-14;
578  unc.en=f.en*1.0e-14;
579  last_method=4;
580  return;
581  }
582  }
583 
584  // Try the degenerate expansion if psi is large enough
585  if (use_expansions && psi>20.0) {
586  bool acc=this->calc_mu_deg(f,temper,1.0e-14);
587  if (acc) {
588  unc.n=f.n*1.0e-14;
589  unc.ed=f.ed*1.0e-14;
590  unc.pr=f.pr*1.0e-14;
591  unc.en=f.en*1.0e-14;
592  last_method=5;
593  return;
594  }
595  }
596 
597  if (!deg) {
598 
599  // If the temperature is large enough, perform the full integral
600 
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);
610 
611  fp_t prefac=f.g*pow(temper,3.0)/2.0/o2scl_const::pi2;
612 
613  // Compute the number density
614 
615  f.n=nit->integ(mfd,0.0,0.0);
616  f.n*=prefac;
617  unc.n=nit->get_error()*prefac;
618 
619  // Compute the energy density
620 
621  f.ed=nit->integ(mfe,0.0,0.0);
622  f.ed*=prefac*temper;
623  if (!f.inc_rest_mass) f.ed-=f.n*f.m;
624  unc.ed=nit->get_error()*prefac*temper;
625 
626  // Compute the entropy
627 
628  f.en=nit->integ(mfs,0.0,0.0);
629  f.en*=prefac;
630  unc.en=nit->get_error()*prefac;
631 
632  last_method=6;
633 
634  } else {
635 
636  // Otherwise, apply a degenerate approximation, by making the
637  // upper integration limit finite
638 
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);
648 
649  fp_t prefac=f.g/2.0/o2scl_const::pi2;
650 
651  // Compute the upper limit for degenerate integrals
652 
653  fp_t arg;
654  if (f.inc_rest_mass) {
655  arg=pow(upper_limit_fac*temper+f.nu,2.0)-f.ms*f.ms;
656  } else {
657  arg=pow(upper_limit_fac*temper+f.nu+f.m,2.0)-f.ms*f.ms;
658  }
659  fp_t ul;
660  if (arg>0.0) {
661  ul=sqrt(arg);
662  } else {
663  f.n=0.0;
664  f.ed=0.0;
665  f.pr=0.0;
666  f.en=0.0;
667  unc.n=0.0;
668  unc.ed=0.0;
669  unc.pr=0.0;
670  unc.en=0.0;
671  O2SCL_ERR2("Zero density in degenerate limit in fermion_rel::",
672  "calc_mu(). Variable deg_limit set improperly?",
673  exc_efailed);
674  return;
675  }
676 
677  // Compute the number density
678 
679  f.n=dit->integ(mfd,0.0,ul);
680  f.n*=prefac;
681  unc.n=dit->get_error()*prefac;
682 
683  // Compute the energy density
684 
685  f.ed=dit->integ(mfe,0.0,ul);
686  f.ed*=prefac;
687  unc.ed=dit->get_error()*prefac;
688 
689  // Compute the lower limit for the entropy integration
690 
691  fp_t ll;
692  if (f.inc_rest_mass) {
693  arg=pow(-upper_limit_fac*temper+f.nu,2.0)-f.ms*f.ms;
694  if (arg>0.0 && (f.ms-f.nu)/temper<-upper_limit_fac) {
695  ll=sqrt(arg);
696  } else {
697  ll=-1.0;
698  }
699  } else {
700  arg=pow(-upper_limit_fac*temper+f.nu+f.m,2.0)-f.ms*f.ms;
701  if (arg>0.0 && (f.ms-f.nu-f.m)/temper<-upper_limit_fac) {
702  ll=sqrt(arg);
703  } else {
704  ll=-1.0;
705  }
706  }
707 
708  // Compute the entropy
709 
710  if (ll>0.0) {
711  f.en=dit->integ(mfs,ll,ul);
712  last_method=7;
713  } else {
714  f.en=dit->integ(mfs,0.0,ul);
715  last_method=8;
716  }
717  f.en*=prefac;
718  unc.en=dit->get_error()*prefac;
719 
720  }
721 
722  // Compute the pressure using the thermodynamic identity
723 
724  f.pr=-f.ed+temper*f.en+f.nu*f.n;
725  unc.pr=sqrt(unc.ed*unc.ed+temper*unc.en*temper*unc.en+
726  f.nu*unc.n*f.nu*unc.n);
727 
728  return;
729  }
730 
731  /** \brief Calculate properties as function of density
732  (template version)
733 
734  \future There is still quite a bit of code duplication
735  between this function and \ref calc_mu_tlate() .
736  */
737  template<class fermion_t>
738  int calc_density_tlate(fermion_t &f, fp_t temper) {
739 
740  last_method=0;
741 
742  // The code below may modify the density, which is confusing to
743  // the user, so we store it here so we can restore it before
744  // this function returns
745  fp_t density_temp=f.n;
746 
747  // -----------------------------------------------------------------
748  // Handle T<=0
749 
750  if (temper<0.0) {
751  O2SCL_ERR2("Temperature less than zero in ",
752  "fermion_rel::calc_density_tlate().",
753  exc_einval);
754  }
755  if (temper==0.0) {
756  last_method=40;
757  this->calc_density_zerot(f);
758  return 0;
759  }
760 
761 #if !O2SCL_NO_RANGE_CHECK
762  // This may not be strictly necessary, because it should be clear
763  // that this function will produce gibberish if the density is not
764  // finite, but I've found this extra checking of the inputs useful
765  // for debugging.
766  if (!std::isfinite(f.n)) {
767  O2SCL_ERR2("Density not finite in ",
768  "fermion_rel::calc_density_tlate().",exc_einval);
769  }
770 #endif
771 
772  // -----------------------------------------------------------------
773  // First determine the chemical potential by solving for the density
774 
775  if (f.non_interacting==true) { f.nu=f.mu; f.ms=f.m; }
776 
777  int ret=nu_from_n(f,temper);
778  last_method*=10;
779  if (ret!=0) {
780  O2SCL_CONV2_RET("Function calc_density() failed in fermion_rel::",
781  "calc_density().",exc_efailed,this->err_nonconv);
782  }
783 
784  if (f.non_interacting) { f.mu=f.nu; }
785 
786  // -----------------------------------------------------------------
787  // Now use the chemical potential to compute the energy density,
788  // pressure, and entropy
789 
790  bool deg=true;
791  fp_t psi;
792  if (f.inc_rest_mass) {
793  psi=(f.nu-f.ms)/temper;
794  } else {
795  psi=(f.nu+(f.m-f.ms))/temper;
796  }
797  if (psi<deg_limit) deg=false;
798 
799  // Try the non-degenerate expansion if psi is small enough
800  if (use_expansions && psi<min_psi) {
801  bool acc=this->calc_mu_ndeg(f,temper,1.0e-14);
802  if (acc) {
803  unc.ed=f.ed*1.0e-14;
804  unc.pr=f.pr*1.0e-14;
805  unc.en=f.en*1.0e-14;
806  f.n=density_temp;
807  last_method+=1;
808  return 0;
809  }
810  }
811 
812  // Try the degenerate expansion if psi is large enough
813  if (use_expansions && psi>20.0) {
814  bool acc=this->calc_mu_deg(f,temper,1.0e-14);
815  if (acc) {
816  unc.n=f.n*1.0e-14;
817  unc.ed=f.ed*1.0e-14;
818  unc.pr=f.pr*1.0e-14;
819  unc.en=f.en*1.0e-14;
820  f.n=density_temp;
821  last_method+=2;
822  return 0;
823  }
824  }
825 
826  if (!deg) {
827 
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);
834 
835  f.ed=nit->integ(mfe,0.0,0.0);
836  f.ed*=f.g*pow(temper,4.0)/2.0/o2scl_const::pi2;
837  if (!f.inc_rest_mass) f.ed-=f.n*f.m;
838  unc.ed=nit->get_error()*f.g*pow(temper,4.0)/2.0/o2scl_const::pi2;
839 
840  f.en=nit->integ(mfs,0.0,0.0);
841  f.en*=f.g*pow(temper,3.0)/2.0/o2scl_const::pi2;
842  unc.en=nit->get_error()*f.g*pow(temper,3.0)/2.0/o2scl_const::pi2;
843  last_method+=3;
844 
845  } else {
846 
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);
853 
854  fp_t arg;
855  if (f.inc_rest_mass) {
856  arg=pow(upper_limit_fac*temper+f.nu,2.0)-f.ms*f.ms;
857  } else {
858  arg=pow(upper_limit_fac*temper+f.nu+f.m,2.0)-f.ms*f.ms;
859  }
860  fp_t ul;
861  if (arg>0.0) {
862 
863  ul=sqrt(arg);
864 
865  fp_t ll;
866  if (f.inc_rest_mass) {
867  arg=pow(-upper_limit_fac*temper+f.nu,2.0)-f.ms*f.ms;
868  if (arg>0.0 && (f.ms-f.nu)/temper<-upper_limit_fac) {
869  ll=sqrt(arg);
870  } else {
871  ll=-1.0;
872  }
873  } else {
874  arg=pow(-upper_limit_fac*temper+f.nu+f.m,2.0)-f.ms*f.ms;
875  if (arg>0.0 && (f.ms-f.nu-f.m)/temper<-upper_limit_fac) {
876  ll=sqrt(arg);
877  } else {
878  ll=-1.0;
879  }
880  }
881 
882  f.ed=dit->integ(mfe,0.0,ul);
883  f.ed*=f.g/2.0/o2scl_const::pi2;
884  unc.ed=dit->get_error()*f.g/2.0/o2scl_const::pi2;
885 
886  if (ll>0.0) {
887  f.en=dit->integ(mfs,ll,ul);
888  last_method+=4;
889  } else {
890  f.en=dit->integ(mfs,0.0,ul);
891  last_method+=5;
892  }
893  f.en*=f.g/2.0/o2scl_const::pi2;
894  unc.en=dit->get_error()*f.g/2.0/o2scl_const::pi2;
895 
896  } else {
897 
898  f.ed=0.0;
899  f.en=0.0;
900  unc.ed=0.0;
901  unc.en=0.0;
902  O2SCL_ERR2("Zero density in degenerate limit in fermion_rel::",
903  "calc_mu(). Variable deg_limit set improperly?",
904  exc_efailed);
905 
906  }
907  }
908 
909  f.n=density_temp;
910  f.pr=-f.ed+temper*f.en+f.nu*f.n;
911  unc.pr=sqrt(unc.ed*unc.ed+temper*unc.en*temper*unc.en+
912  f.nu*unc.n*f.nu*unc.n);
913 
914  return 0;
915  }
916 
917  /** \brief Calculate properties with antiparticles as function of
918  chemical potential (template version)
919  */
920  template<class fermion_t>
921  void pair_mu_tlate(fermion_t &f, fp_t temper) {
922 
923  last_method=0;
924 
925  if (f.non_interacting) { f.nu=f.mu; f.ms=f.m; }
926 
927  // AWS: 2/12/19: I'm taking this out, similar to the removal
928  // of the code in fermion_rel::pair_fun(). If I put it
929  // back in, I need to find a new value for last_method
930  // other than 9.
931  if (false && use_expansions) {
932  if (this->calc_mu_ndeg(f,temper,1.0e-14,true)) {
933  unc.n=1.0e-14*f.n;
934  unc.ed=1.0e-14*f.ed;
935  unc.en=1.0e-14*f.en;
936  unc.pr=1.0e-14*f.pr;
937  last_method=9;
938  return;
939  }
940  }
941 
942  fermion antip(f.m,f.g);
943  f.anti(antip);
944 
945  // Particles
946  calc_mu(f,temper);
947  fp_t unc_n=unc.n;
948  fp_t unc_pr=unc.pr;
949  fp_t unc_ed=unc.ed;
950  fp_t unc_en=unc.en;
951 
952  // Antiparticles
953  int lm=last_method*10;
954  calc_mu(antip,temper);
955  last_method+=lm;
956  last_method*=10;
957 
958  // Add up thermodynamic quantities
959  if (f.inc_rest_mass) {
960  f.ed+=antip.ed;
961  } else {
962  f.ed=f.ed+antip.ed+2.0*antip.n*f.m;
963  }
964  f.n-=antip.n;
965  f.pr+=antip.pr;
966  f.en+=antip.en;
967 
968  // Add up uncertainties
969  unc.n=gsl_hypot(unc.n,unc_n);
970  unc.ed=gsl_hypot(unc.ed,unc_ed);
971  unc.pr=gsl_hypot(unc.pr,unc_pr);
972  unc.en=gsl_hypot(unc.ed,unc_en);
973 
974  return;
975  }
976 
977  /** \brief Calculate thermodynamic properties with antiparticles
978  from the density (template version)
979  */
980  template<class fermion_t>
981  int pair_density_tlate(fermion_t &f, fp_t temper) {
982 
983  last_method=0;
984 
985  // -----------------------------------------------------------------
986  // Handle T<=0
987 
988  if (temper<=0.0) {
989  this->calc_density_zerot(f);
990  last_method=1000;
991  return success;
992  }
993 
994  // Storage the input density
995  fp_t density_temp=f.n;
996 
997  if (f.non_interacting==true) { f.nu=f.mu; f.ms=f.m; }
998 
999  fp_t initial_guess=f.nu;
1000 
1001  fp_t nex=f.nu/temper;
1002 
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);
1006 
1007  // Begin by trying the user-specified guess
1008  bool drec=density_root->err_nonconv;
1009  density_root->err_nonconv=false;
1010  int ret=density_root->solve(nex,mf);
1011 
1012  // If that doesn't work, try bracketing the root
1013  if (ret==0) {
1014  last_method=2000;
1015  } else {
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++) {
1020  bhigh*=1.0e2;
1021  yhigh=mf(bhigh);
1022  }
1023  for(size_t j=0;j<5 && ylow>0.0;j++) {
1024  blow*=1.0e2;
1025  ylow=mf(blow);
1026  }
1027  if (yhigh>0.0 && ylow<0.0) {
1028  root_brent_gsl<> rbg;
1029  rbg.err_nonconv=false;
1030  ret=rbg.solve_bkt(blow,bhigh,mf);
1031  if (ret==0) {
1032  nex=blow;
1033  last_method=3000;
1034  }
1035  }
1036  }
1037 
1038  if (ret!=0) {
1039 
1040  // If that fails, try to make the integrators more accurate
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;
1047  ret=density_root->solve(nex,mf);
1048  if (ret==0) last_method=4000;
1049 
1050  // AWS: 7/25/18: We work in log units below, so we ensure the
1051  // chemical potential is not negative
1052  if (nex<0.0) nex=1.0e-10;
1053 
1054  // If that failed, try working in log units
1055 
1056  // Function in log units
1057  funct lmf=std::bind(std::mem_fn<fp_t(fp_t,fermion &,
1058  fp_t,bool)>
1060  this,std::placeholders::_1,std::ref(f),
1061  temper,true);
1062 
1063  if (ret!=0) {
1064  nex=log(nex);
1065  ret=density_root->solve(nex,lmf);
1066  nex=exp(nex);
1067  if (ret==0) last_method=5000;
1068  }
1069 
1070  if (ret!=0) {
1071  // If that failed, try a different solver
1072  root_brent_gsl<> rbg;
1073  rbg.err_nonconv=false;
1074  nex=log(nex);
1075  ret=rbg.solve(nex,lmf);
1076  nex=exp(nex);
1077  if (ret==0) last_method=6000;
1078  }
1079 
1080  // Return tolerances to their original values
1081  dit->tol_rel=tol1;
1082  dit->tol_abs=tol2;
1083  nit->tol_rel=tol3;
1084  nit->tol_abs=tol4;
1085  }
1086 
1087  // Restore value of err_nonconv
1088  density_root->err_nonconv=drec;
1089 
1090  if (ret!=0) {
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;
1095  O2SCL_CONV2_RET("Density solver failed in fermion_rel::",
1096  "pair_density().",exc_efailed,this->err_nonconv);
1097  }
1098 
1099  f.nu=nex*temper;
1100 
1101  if (f.non_interacting==true) { f.mu=f.nu; }
1102 
1103  // Finally, now that we have the chemical potential, use pair_mu()
1104  // to evaluate the energy density, pressure, and entropy
1105  int lm=last_method;
1106  pair_mu(f,temper);
1107  last_method+=lm;
1108 
1109  // The function pair_mu() can modify the density, which would be
1110  // confusing to the user, so we return it to the user-specified
1111  // value.
1112  f.n=density_temp;
1113 
1114  return success;
1115  }
1116  //@}
1117 
1118  protected:
1119 
1120 #ifndef DOXYGEN_INTERNAL
1121 
1122  /// The integrand for the density for non-degenerate fermions
1123  fp_t density_fun(fp_t u, fermion &f, fp_t T) {
1124 
1125  double ret, y, eta;
1126 
1127  if (f.inc_rest_mass) {
1128  y=f.nu/T;
1129  } else {
1130  y=(f.nu+f.m)/T;
1131  }
1132  eta=f.ms/T;
1133 
1134  if (y-eta-u>exp_limit) {
1135  ret=(eta+u)*sqrt(u*u+2.0*eta*u);
1136  } else if (y>u+exp_limit && eta>u+exp_limit) {
1137  ret=(eta+u)*sqrt(u*u+2.0*eta*u)/(exp(eta+u-y)+1.0);
1138  } else {
1139  ret=(eta+u)*sqrt(u*u+2.0*eta*u)*exp(y)/(exp(eta+u)+exp(y));
1140  }
1141 
1142  if (!std::isfinite(ret)) {
1143  ret=0.0;
1144  }
1145 
1146  return ret;
1147  }
1148 
1149 
1150  /// The integrand for the energy density for non-degenerate fermions
1151  fp_t energy_fun(fp_t u, fermion &f, fp_t T) {
1152  double ret, y, eta;
1153 
1154  eta=f.ms/T;
1155 
1156  if (f.inc_rest_mass) {
1157  y=f.nu/T;
1158  } else {
1159  y=(f.nu+f.m)/T;
1160  }
1161  if (y>u+exp_limit && eta>u+exp_limit) {
1162  ret=(eta+u)*(eta+u)*sqrt(u*u+2.0*eta*u)/(exp(eta+u-y)+1.0);
1163  } else {
1164  ret=(eta+u)*(eta+u)*sqrt(u*u+2.0*eta*u)*exp(y)/(exp(eta+u)+exp(y));
1165  }
1166 
1167  if (!std::isfinite(ret)) {
1168  return 0.0;
1169  }
1170 
1171  return ret;
1172  }
1173 
1174  /// The integrand for the entropy density for non-degenerate fermions
1175  fp_t entropy_fun(fp_t u, fermion &f, fp_t T) {
1176  double ret, y, eta, term1, term2;
1177 
1178  if (f.inc_rest_mass) {
1179  y=f.nu/T;
1180  } else {
1181  y=(f.nu+f.m)/T;
1182  }
1183  eta=f.ms/T;
1184 
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);
1188 
1189  if (!std::isfinite(ret)) {
1190  return 0.0;
1191  }
1192 
1193  return ret;
1194  }
1195 
1196 
1197  /// The integrand for the density for degenerate fermions
1198  fp_t deg_density_fun(fp_t k, fermion &f, fp_t T) {
1199 
1200  double E=gsl_hypot(k,f.ms), ret;
1201  if (!f.inc_rest_mass) E-=f.m;
1202 
1203  ret=k*k/(1.0+exp((E-f.nu)/T));
1204 
1205  if (!std::isfinite(ret)) {
1206  O2SCL_ERR2("Returned not finite result ",
1207  "in fermion_rel::deg_density_fun().",exc_einval);
1208  }
1209 
1210  return ret;
1211  }
1212 
1213  /// The integrand for the energy density for degenerate fermions
1214  fp_t deg_energy_fun(fp_t k, fermion &f, fp_t T) {
1215 
1216  double E=gsl_hypot(k,f.ms), ret;
1217  if (!f.inc_rest_mass) E-=f.m;
1218 
1219  ret=k*k*E/(1.0+exp((E-f.nu)/T));
1220 
1221  if (!std::isfinite(ret)) {
1222  O2SCL_ERR2("Returned not finite result ",
1223  "in fermion_rel::deg_energy_fun().",exc_einval);
1224  }
1225 
1226  return ret;
1227  }
1228 
1229  /// The integrand for the entropy density for degenerate fermions
1230  fp_t deg_entropy_fun(fp_t k, fermion &f, fp_t T) {
1231 
1232  double E=gsl_hypot(k,f.ms), ret;
1233  if (!f.inc_rest_mass) E-=f.m;
1234 
1235  // If the argument to the exponential is really small, then the
1236  // value of the integrand is just zero
1237  if (((E-f.nu)/T)<-exp_limit) {
1238  ret=0.0;
1239  // Otherwise, if the argument to the exponential is still small,
1240  // then addition of 1 makes us lose precision, so we use an
1241  // alternative:
1242  } else if (((E-f.nu)/T)<-deg_entropy_fac) {
1243  double arg=E/T-f.nu/T;
1244  ret=-k*k*(-1.0+arg)*exp(arg);
1245  } else {
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));
1248  }
1249 
1250  if (!std::isfinite(ret)) {
1251  O2SCL_ERR2("Returned not finite result ",
1252  "in fermion_rel::deg_entropy_fun().",exc_einval);
1253  }
1254 
1255  return ret;
1256  }
1257 
1258  /// Solve for the chemical potential given the density
1259  fp_t solve_fun(fp_t x, fermion &f, fp_t T) {
1260  double nden, yy;
1261 
1262  f.nu=T*x;
1263 
1264  if (f.non_interacting) f.mu=f.nu;
1265 
1266  bool deg=true;
1267  double psi;
1268  if (f.inc_rest_mass) {
1269  psi=(f.nu-f.ms)/T;
1270  } else {
1271  psi=(f.nu+(f.m-f.ms))/T;
1272  }
1273  if (psi<deg_limit) deg=false;
1274 
1275  // Try the non-degenerate expansion if psi is small enough
1276  if (use_expansions && psi<min_psi) {
1277  double ntemp=f.n;
1278  bool acc=this->calc_mu_ndeg(f,T,1.0e-14);
1279  if (acc) {
1280  unc.n=f.n*1.0e-14;
1281  yy=(ntemp-f.n)/ntemp;
1282  f.n=ntemp;
1283  return yy;
1284  }
1285  f.n=ntemp;
1286  }
1287 
1288  // Try the degenerate expansion if psi is large enough
1289  if (use_expansions && psi>20.0) {
1290  double ntemp=f.n;
1291  bool acc=this->calc_mu_deg(f,T,1.0e-14);
1292  if (acc) {
1293  unc.n=f.n*1.0e-14;
1294  yy=(ntemp-f.n)/ntemp;
1295  f.n=ntemp;
1296  return yy;
1297  }
1298  f.n=ntemp;
1299  }
1300 
1301  // Otherwise, directly perform the integration
1302  if (!deg) {
1303 
1304  funct mfe=std::bind(std::mem_fn<double(double,fermion &,double)>
1306  this,std::placeholders::_1,std::ref(f),T);
1307 
1308  nden=nit->integ(mfe,0.0,0.0);
1309  nden*=f.g*pow(T,3.0)/2.0/o2scl_const::pi2;
1310  unc.n=nit->get_error()*f.g*pow(T,3.0)/2.0/o2scl_const::pi2;
1311 
1312  yy=(f.n-nden)/f.n;
1313 
1314  } else {
1315 
1316  funct mfe=std::bind(std::mem_fn<double(double,fermion &,double)>
1318  this,std::placeholders::_1,std::ref(f),T);
1319 
1320  double arg;
1321  if (f.inc_rest_mass) {
1322  arg=pow(upper_limit_fac*T+f.nu,2.0)-f.ms*f.ms;
1323  } else {
1324  arg=pow(upper_limit_fac*T+f.nu+f.m,2.0)-f.ms*f.ms;
1325  }
1326 
1327  double ul;
1328 
1329  if (arg>0.0) {
1330 
1331  ul=sqrt(arg);
1332 
1333  nden=dit->integ(mfe,0.0,ul);
1334  nden*=f.g/2.0/o2scl_const::pi2;
1335  unc.n=dit->get_error()*f.g/2.0/o2scl_const::pi2;
1336 
1337  } else {
1338 
1339  nden=0.0;
1340 
1341  }
1342 
1343  yy=(f.n-nden)/f.n;
1344  }
1345 
1346  return yy;
1347  }
1348 
1349 
1350  /** \brief Solve for the chemical potential given the density
1351  with antiparticles
1352 
1353  \future Particles and antiparticles have different degeneracy
1354  factors, so we separately use the expansions one at a time. It
1355  is probably better to separately generate a new expansion
1356  function which automatically handles the sum of particles and
1357  antiparticles.
1358  */
1359  fp_t pair_fun(fp_t x, fermion &f, fp_t T, bool log_mode) {
1360 
1361  // Temporary storage for density to match
1362  double nn_match=f.n;
1363  // Number density of particles and antiparticles
1364  double nden_p, nden_ap;
1365 
1366  // -----------------------------------------------------------------
1367 
1368  f.nu=T*x;
1369  if (log_mode) {
1370  f.nu=T*exp(x);
1371  }
1372 
1373  // Sometimes the exp() call above causes an overflow, so
1374  // we avoid extreme values
1375  if (!std::isfinite(f.nu)) return 3;
1376 
1377  if (f.non_interacting) f.mu=f.nu;
1378 
1379  // -----------------------------------------------------------------
1380  // First, try the non-degenerate expansion with both particles and
1381  // antiparticles together
1382 
1383  // AWS: 7/25/18: I'm taking this section out because it doesn't seem
1384  // to make sense to me, it apparently uses calc_mu_ndeg() which is
1385  // for particles only, and I'm not sure that's sufficient here. This
1386  // section also caused problems for the n=0, T!=0 case.
1387 
1388  if (false && use_expansions) {
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().",
1393  exc_einval);
1394  }
1395  // Make sure to restore the value of f.n to it's original value,
1396  // nn_match
1397  f.n=nn_match;
1398  return y1;
1399  }
1400  }
1401 
1402  // -----------------------------------------------------------------
1403  // If that doesn't work, evaluate particles and antiparticles
1404  // separately. This is the contribution for particles
1405 
1406  bool deg=true;
1407  double psi;
1408  if (f.inc_rest_mass) {
1409  psi=(f.nu-f.ms)/T;
1410  } else {
1411  psi=(f.nu+(f.m-f.ms))/T;
1412  }
1413  if (psi<deg_limit) deg=false;
1414 
1415  bool particles_done=false;
1416 
1417  // Try the non-degenerate expansion if psi is small enough
1418  if (use_expansions && psi<min_psi) {
1419  if (this->calc_mu_ndeg(f,T,1.0e-8) && std::isfinite(f.n)) {
1420  particles_done=true;
1421  nden_p=f.n;
1422  if (!std::isfinite(nden_p)) {
1423  O2SCL_ERR("Value 'nden_p' not finite (1) in fermion_rel::pair_fun().",
1424  exc_einval);
1425  }
1426  }
1427  }
1428 
1429  // Try the degenerate expansion if psi is large enough
1430  if (use_expansions && particles_done==false && psi>20.0) {
1431  if (this->calc_mu_deg(f,T,1.0e-8) && std::isfinite(f.n)) {
1432  particles_done=true;
1433  nden_p=f.n;
1434  if (!std::isfinite(nden_p)) {
1435  O2SCL_ERR("Value 'nden_p' not finite (2) in fermion_rel::pair_fun().",
1436  exc_einval);
1437  }
1438  }
1439  }
1440 
1441  // If neither expansion worked, use direct integration
1442  if (particles_done==false) {
1443 
1444  if (!deg) {
1445 
1446  // Nondegenerate case
1447 
1448  funct mfe=std::bind(std::mem_fn<double(double,fermion &,double)>
1450  this,std::placeholders::_1,std::ref(f),T);
1451 
1452  nden_p=nit->integ(mfe,0.0,0.0);
1453  nden_p*=f.g*pow(T,3.0)/2.0/o2scl_const::pi2;
1454  if (!std::isfinite(nden_p)) {
1455  O2SCL_ERR("Value 'nden_p' not finite (3) in fermion_rel::pair_fun().",
1456  exc_einval);
1457  }
1458 
1459  } else {
1460 
1461  // Degenerate case
1462 
1463  funct mfe=std::bind(std::mem_fn<double(double,fermion &,double)>
1465  this,std::placeholders::_1,std::ref(f),T);
1466 
1467  double arg;
1468  if (f.inc_rest_mass) {
1469  arg=pow(upper_limit_fac*T+f.nu,2.0)-f.ms*f.ms;
1470  } else {
1471  arg=pow(upper_limit_fac*T+f.nu+f.m,2.0)-f.ms*f.ms;
1472  }
1473 
1474  double ul;
1475  if (arg>0.0) {
1476  ul=sqrt(arg);
1477  nden_p=dit->integ(mfe,0.0,ul);
1478  nden_p*=f.g/2.0/o2scl_const::pi2;
1479  } else {
1480  nden_p=0.0;
1481  }
1482 
1483  if (!std::isfinite(nden_p)) {
1484  O2SCL_ERR("Value 'nden_p' not finite (4) in fermion_rel::pair_fun().",
1485  exc_einval);
1486  }
1487 
1488  }
1489 
1490  particles_done=true;
1491 
1492  // End of 'if (particles_done==false)'
1493  }
1494 
1495  // -----------------------------------------------------------------
1496  // Compute the contribution from the antiparticles
1497 
1498  if (f.inc_rest_mass) {
1499  f.nu=-T*x;
1500  if (log_mode) f.nu=-T*exp(x);
1501  } else {
1502  f.nu=-T*x-2.0*f.m;
1503  if (log_mode) f.nu=-T*exp(x)-2.0*f.m;
1504  }
1505  if (f.non_interacting) f.mu=f.nu;
1506 
1507  bool antiparticles_done=false;
1508 
1509  // Evaluate the degeneracy parameter
1510  deg=true;
1511  if (f.inc_rest_mass) {
1512  psi=(f.nu-f.ms)/T;
1513  } else {
1514  psi=(f.nu+f.m-f.ms)/T;
1515  }
1516  if (psi<deg_limit) deg=false;
1517 
1518  // Try the non-degenerate expansion if psi is small enough
1519  if (use_expansions && psi<min_psi) {
1520  if (this->calc_mu_ndeg(f,T,1.0e-8)) {
1521  antiparticles_done=true;
1522  nden_ap=f.n;
1523  if (!std::isfinite(nden_ap)) {
1524  O2SCL_ERR2("Value 'nden_ap' not finite (5) in",
1525  "fermion_rel::pair_fun().",
1526  exc_einval);
1527  }
1528  }
1529  }
1530 
1531  // Try the degenerate expansion if psi is large enough
1532  if (use_expansions && antiparticles_done==false && psi>20.0) {
1533  if (this->calc_mu_deg(f,T,1.0e-8)) {
1534  antiparticles_done=true;
1535  nden_ap=f.n;
1536  if (!std::isfinite(nden_ap)) {
1537  O2SCL_ERR2("Value 'nden_ap' not finite (6) in",
1538  "fermion_rel::pair_fun().",
1539  exc_einval);
1540  }
1541  }
1542  }
1543 
1544  // If neither expansion worked, use direct integration
1545  if (antiparticles_done==false) {
1546 
1547  if (!deg) {
1548 
1549  // Nondegenerate case
1550 
1551  funct mf=std::bind(std::mem_fn<double(double,fermion &,double)>
1553  this,std::placeholders::_1,std::ref(f),T);
1554 
1555  nden_ap=nit->integ(mf,0.0,0.0);
1556  nden_ap*=f.g*pow(T,3.0)/2.0/o2scl_const::pi2;
1557  if (!std::isfinite(nden_ap)) {
1558  O2SCL_ERR2("Value 'nden_ap' not finite (7) in",
1559  "fermion_rel::pair_fun().",
1560  exc_einval);
1561  }
1562 
1563  } else {
1564 
1565  // Degenerate case
1566 
1567  funct mf=std::bind(std::mem_fn<double(double,fermion &,double)>
1569  this,std::placeholders::_1,std::ref(f),T);
1570 
1571  double arg;
1572  if (f.inc_rest_mass) {
1573  arg=pow(upper_limit_fac*T+f.nu,2.0)-f.ms*f.ms;
1574  } else {
1575  arg=pow(upper_limit_fac*T+f.nu+f.m,2.0)-f.ms*f.ms;
1576  }
1577 
1578  double ul;
1579  if (arg>0.0) {
1580  ul=sqrt(arg);
1581  nden_ap=dit->integ(mf,0.0,ul);
1582  nden_ap*=f.g/2.0/o2scl_const::pi2;
1583  } else {
1584  nden_ap=0.0;
1585  }
1586  if (!std::isfinite(nden_ap)) {
1587  O2SCL_ERR2("Value 'nden_ap' not finite (8) in",
1588  "fermion_rel::pair_fun().",
1589  exc_einval);
1590  }
1591 
1592  }
1593 
1594  antiparticles_done=true;
1595  }
1596 
1597  double y2;
1598  // Finish computing the function value
1599  if (nn_match==0.0) {
1600  y2=fabs(nden_p-nden_ap)/fabs(nden_p);
1601  } else {
1602  y2=(nden_p-nden_ap)/nn_match-1.0;
1603  }
1604 
1605  if (!std::isfinite(y2)) {
1606  O2SCL_ERR("Value 'y2' not finite (9) in fermion_rel::pair_fun().",
1607  exc_einval);
1608  }
1609 
1610  // Make sure to restore the value of f.n to it's original value,
1611  // nn_match
1612  f.n=nn_match;
1613  return y2;
1614  }
1615 
1616 
1617 #endif
1618 
1619  };
1620 
1621  /** \brief Double-precision version of
1622  \ref o2scl::fermion_rel_tl
1623  */
1625 
1626 #ifndef DOXYGEN_NO_O2NS
1627 }
1628 #endif
1629 
1630 #endif
o2scl::fermion_rel_tl::deg_entropy_fun
fp_t deg_entropy_fun(fp_t k, fermion &f, fp_t T)
The integrand for the entropy density for degenerate fermions.
Definition: fermion_rel.h:1230
o2scl::fermion_rel_tl::calc_density
virtual int calc_density(fermion &f, fp_t temper)
Calculate properties as function of density.
Definition: fermion_rel.h:295
o2scl::fermion_rel_tl::calc_mu_tlate
void calc_mu_tlate(fermion_t &f, fp_t temper)
Calculate properties as function of chemical potential (template version)
Definition: fermion_rel.h:541
O2SCL_CONV2_RET
#define O2SCL_CONV2_RET(d, d2, n, b)
o2scl::fermion_rel_tl::dit
std::shared_ptr< inte<> > dit
The degenerate integrator.
Definition: fermion_rel.h:324
o2scl::inte_qag_gsl
o2scl_const::pi2
const double pi2
o2scl::fermion_rel_tl::unc
fermion unc
Storage for the uncertainty.
Definition: fermion_rel.h:256
o2scl::part_tl::mu
fp_t mu
Chemical potential.
Definition: part.h:118
o2scl::part_tl::ms
fp_t ms
Effective mass (Dirac unless otherwise specified)
Definition: part.h:122
exc_efailed
exc_efailed
o2scl::fermion_rel_tl::deg_energy_fun
fp_t deg_energy_fun(fp_t k, fermion &f, fp_t T)
The integrand for the energy density for degenerate fermions.
Definition: fermion_rel.h:1214
o2scl::fermion_rel_tl::pair_fun
fp_t pair_fun(fp_t x, fermion &f, fp_t T, bool log_mode)
Solve for the chemical potential given the density with antiparticles.
Definition: fermion_rel.h:1359
O2SCL_ERR2
#define O2SCL_ERR2(d, d2, n)
o2scl::fermion_rel_tl::solve_fun
fp_t solve_fun(fp_t x, fermion &f, fp_t T)
Solve for the chemical potential given the density.
Definition: fermion_rel.h:1259
o2scl::part_tl::n
fp_t n
Number density.
Definition: part.h:112
o2scl::fermion_zerot_tl< double >::calc_density_zerot
virtual void calc_density_zerot(fermion_tl< double > &f)
Zero temperature fermions from o2scl::part_tl::n and o2scl::part_tl::ms.
Definition: fermion.h:255
o2scl::root_brent_gsl::solve_bkt
virtual int solve_bkt(fp_t &x1, fp_t x2, func_t &f)
o2scl::fermion_rel_tl::pair_mu_tlate
void pair_mu_tlate(fermion_t &f, fp_t temper)
Calculate properties with antiparticles as function of chemical potential (template version)
Definition: fermion_rel.h:921
o2scl::part_tl::non_interacting
bool non_interacting
True if the particle is non-interacting (default true)
Definition: part.h:130
o2scl::fermion_rel_tl::pair_density_tlate
int pair_density_tlate(fermion_t &f, fp_t temper)
Calculate thermodynamic properties with antiparticles from the density (template version)
Definition: fermion_rel.h:981
o2scl::fermion_rel_tl::density_fun
fp_t density_fun(fp_t u, fermion &f, fp_t T)
The integrand for the density for non-degenerate fermions.
Definition: fermion_rel.h:1123
o2scl::part_tl::en
fp_t en
Entropy density.
Definition: part.h:120
o2scl::fermion_rel_tl::fermion_rel_tl
fermion_rel_tl()
Create a fermion with mass m and degeneracy g.
Definition: fermion_rel.h:262
o2scl::inte_qagiu_gsl
o2scl::part_tl::g
fp_t g
Degeneracy (e.g. spin and color if applicable)
Definition: part.h:108
o2scl::fermion_thermo_tl< double >::calc_mu_ndeg
virtual bool calc_mu_ndeg(fermion &f, double temper, double prec=1.0e-18, bool inc_antip=false)
Non-degenerate expansion for fermions.
Definition: fermion.h:599
o2scl::part_tl::ed
fp_t ed
Energy density.
Definition: part.h:114
root_bkt< funct, funct, double >::solve
virtual int solve(double &x, funct &func)
o2scl::fermion_rel_tl::min_psi
fp_t min_psi
The smallest value of for which integration is used.
Definition: fermion_rel.h:233
o2scl::fermion_rel_tl::err_nonconv
bool err_nonconv
If true, call the error handler when convergence fails (default true)
Definition: fermion_rel.h:228
o2scl::fermion_rel_tl
Equation of state for a relativistic fermion.
Definition: fermion_rel.h:219
success
success
o2scl::fermion_rel_tl::deg_limit
fp_t deg_limit
The critical degeneracy at which to switch integration techniques (default 2)
Definition: fermion_rel.h:238
o2scl::fermion_rel_tl::calc_mu
virtual void calc_mu(fermion &f, fp_t temper)
Calculate properties as function of chemical potential.
Definition: fermion_rel.h:283
o2scl::root_brent_gsl
o2scl::fermion_tl< double >
exc_einval
exc_einval
o2scl::part_tl::inc_rest_mass
bool inc_rest_mass
If true, include the mass in the energy density and chemical potential (default true)
Definition: part.h:128
o2scl::fermion_rel_tl::use_expansions
bool use_expansions
If true, use expansions for extreme conditions (default true)
Definition: fermion_rel.h:259
o2scl::fermion_rel_tl::deg_density_fun
fp_t deg_density_fun(fp_t k, fermion &f, fp_t T)
The integrand for the density for degenerate fermions.
Definition: fermion_rel.h:1198
o2scl::part_tl::m
fp_t m
Mass.
Definition: part.h:110
o2scl::part_tl::pr
fp_t pr
Pressure.
Definition: part.h:116
o2scl::fermion_rel_tl::last_method
int last_method
An integer indicating the last numerical method used.
Definition: fermion_rel.h:383
o2scl::fermion_rel_tl::nu_from_n_tlate
int nu_from_n_tlate(fermion_t &f, fp_t temper)
Calculate the chemical potential from the density (template version)
Definition: fermion_rel.h:391
o2scl::fermion_rel_tl::type
virtual const char * type()
Return string denoting type ("fermion_rel")
Definition: fermion_rel.h:330
o2scl::root_cern
o2scl::fermion_rel_tl::pair_density
virtual int pair_density(fermion &f, fp_t temper)
Calculate properties with antiparticles as function of density.
Definition: fermion_rel.h:310
o2scl::fermion_thermo_tl< double >::calc_mu_deg
virtual bool calc_mu_deg(fermion &f, double temper, double prec=1.0e-18)
Degenerate expansion for fermions.
Definition: fermion.h:621
o2scl::fermion_rel_tl::energy_fun
fp_t energy_fun(fp_t u, fermion &f, fp_t T)
The integrand for the energy density for non-degenerate fermions.
Definition: fermion_rel.h:1151
o2scl::fermion_rel_tl::deg_entropy_fac
fp_t deg_entropy_fac
A factor for the degenerate entropy integration (default 30)
Definition: fermion_rel.h:249
O2SCL_ERR
#define O2SCL_ERR(d, n)
o2scl::funct
std::function< double(double)> funct
o2scl::fermion_rel_tl::exp_limit
fp_t exp_limit
The limit for exponentials to ensure integrals are finite (default 200)
Definition: fermion_rel.h:243
o2scl::fermion_thermo_tl
Fermion with finite-temperature thermodynamics [abstract base].
Definition: fermion.h:309
o2scl::fermion_zerot_tl< double >::calc_mu_zerot
virtual void calc_mu_zerot(fermion_tl< double > &f)
Zero temperature fermions from o2scl::part_tl::mu or o2scl::part_tl::nu and o2scl::part_tl::ms.
Definition: fermion.h:217
o2scl::fermion_rel_tl::entropy_fun
fp_t entropy_fun(fp_t u, fermion &f, fp_t T)
The integrand for the entropy density for non-degenerate fermions.
Definition: fermion_rel.h:1175
root< funct, funct, double >::err_nonconv
bool err_nonconv
o2scl::fermion_rel_tl::nu_from_n
virtual int nu_from_n(fermion &f, fp_t temper)
Calculate effective chemical potential from density.
Definition: fermion_rel.h:316
o2scl::part_tl::nu
fp_t nu
Effective chemical potential.
Definition: part.h:124
o2scl::fermion_rel_tl::nit
std::shared_ptr< inte<> > nit
The non-degenerate integrator.
Definition: fermion_rel.h:321
o2scl::fermion_rel_tl::density_root
std::shared_ptr< root<> > density_root
The solver for calc_density()
Definition: fermion_rel.h:327
psi
const double psi
o2scl::fermion_rel_tl::upper_limit_fac
fp_t upper_limit_fac
The factor for the degenerate upper limits (default 20)
Definition: fermion_rel.h:246
o2scl::fermion_rel_tl::pair_mu
virtual void pair_mu(fermion &f, fp_t temper)
Calculate properties with antiparticles as function of chemical potential.
Definition: fermion_rel.h:302
o2scl::fermion_rel_tl::calc_density_tlate
int calc_density_tlate(fermion_t &f, fp_t temper)
Calculate properties as function of density (template version)
Definition: fermion_rel.h:738
o2scl::fermion_rel_tl::verbose
int verbose
Verbosity parameter (default 0)
Definition: fermion_rel.h:252

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).