eos_had_base.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2018, Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 /** \file eos_had_base.h
24  \brief File defining \ref o2scl::eos_had_base
25 */
26 #ifndef O2SCL_HADRONIC_EOS_H
27 #define O2SCL_HADRONIC_EOS_H
28 
29 #include <iostream>
30 #include <string>
31 
32 #include <boost/numeric/ublas/vector.hpp>
33 
34 #include <o2scl/deriv_gsl.h>
35 #include <o2scl/mroot.h>
36 #include <o2scl/mroot_hybrids.h>
37 #include <o2scl/mm_funct.h>
38 #include <o2scl/eos_base.h>
39 #include <o2scl/fermion_eff.h>
40 #include <o2scl/part.h>
41 #include <o2scl/lib_settings.h>
42 
43 #ifndef DOXYGEN_NO_O2NS
44 namespace o2scl {
45 #endif
46 
47  /** \brief Hadronic equation of state [abstract base]
48 
49  \note This class and all of its children expect the neutron and
50  proton properties to be specified in powers of \f$ \mathrm{fm}
51  \f$ so that masses and chemical potentials are in powers of \f$
52  \mathrm{fm}^{-1} \f$ and energy densities and pressures are in
53  powers of \f$ \mathrm{fm}^{-4} \f$ .
54 
55  Denote the number density of neutrons as \f$ n_n \f$, the number
56  density of protons as \f$ n_p \f$, the total baryon density \f$
57  n_B = n_n + n_p \f$, the asymmetry \f$ \delta \equiv
58  (n_n-n_p)/n_B \f$, the nuclear saturation density as \f$ n_0
59  \approx 0.16~\mathrm{fm}^{-3} \f$, and the quantity \f$ \epsilon
60  \equiv (n_B-n_0)/3n_0 \f$. (Note that some authors define
61  \f$ \delta \f$ as \f$ n_n - n_p \f$, which is not the same as
62  the definition above.) Then the energy per baryon of nucleonic
63  matter can be written as an expansion around
64  \f$ \epsilon =\delta = 0 \f$
65  \f[
66  E(n_B,\delta) = -B + \frac{\tilde{K}}{2!} {\epsilon}^2 +
67  \frac{Q}{3!} {\epsilon}^3 + \delta^2 \left(S + L \epsilon +
68  \frac{K_{\mathrm{sym}}}{2!} {\epsilon}^2
69  + \frac{Q_{\mathrm{sym}}}{3!} \epsilon^3 \right) +
70  E_4(n_B,\delta) + {\cal O}(\delta^6)
71  \qquad \left(\mathrm{Eq.}~1\right)
72  \f]
73  where \f$ E_4 \f$ represents the quartic terms
74  \f[
75  E_4(n_B,\delta)
76  = \delta^4 \left(S_4 + L_4 \epsilon + \frac{K_4}{2!} {\epsilon}^2
77  + \frac{Q_4}{3!} \epsilon^3 \right) \qquad
78  \left(\mathrm{Eq.}~2\right)
79  \f]
80  (Adapted slightly from \ref Piekarewicz09). From this, one can
81  compute the energy density of nuclear matter \f$
82  \varepsilon(n_B,\delta) = n_B E(n_B,\delta) \f$, the chemical
83  potentials \f$ \mu_i \equiv (\partial \varepsilon) / (\partial
84  n_i ) \f$ and the pressure \f$ P = -\varepsilon + \mu_n n_n +
85  \mu_p n_p \f$. This expansion motivates the definition of
86  several separate terms. The binding energy \f$ B \f$ of
87  symmetric nuclear matter (\f$ \delta = 0 \f$) is around 16 MeV.
88 
89  The compression modulus is usually defined by \f$ \chi = -1/V
90  (dV/dP) = 1/n (dP/dn)^{-1} \f$ . In nuclear physics it has
91  become common to use the incompressibility (or bulk) modulus
92  with an extra factor of 9 (motivated by the 3 in the denominator
93  in the definition of \f$ \epsilon \f$), \f$ K=9/(n \chi) \f$ and
94  refer to \f$ K \f$ simply as the incompressibility. Here, we
95  define the function
96  \f[
97  K(n_B,\delta) \equiv 9 \left( \frac{\partial P}{\partial n_B}
98  \right) = 9 n_B \left(\frac{\partial^2 \varepsilon}
99  {\partial n_B^2}\right)
100  \f]
101  This quantity is computed by the function \ref fcomp() by
102  computing the first derivative of the pressure, which is more
103  numerically stable than the second derivative of the energy
104  density (since most \o2 EOSs compute the pressure exactly). This
105  function is typically evaluated at the point \f$
106  (n_B=n_0,\delta=0) \f$ and is stored in \ref comp by the
107  function \ref saturation(). This quantity is not always the same
108  as \f$ \tilde{K} \f$, defined here as
109  \f[
110  \tilde{K}(n_B,\delta) = 9 n_B^2 \left(\frac{\partial^2 E}{\partial
111  n_B^2}\right) = K(n_B,\delta) - \frac{1}{n_B} 18 P(n_B,\delta)
112  \f]
113  We denote \f$ K \equiv K(n_B=n_0,\delta=0) \f$ and similarly for
114  \f$ \tilde{K} \f$, the quantity in Eq. 1 above. In nuclear
115  matter at saturation, the pressure is zero and \f$ K = \tilde{K}
116  \f$. See \ref Chabanat97 for further discussion of the distinction
117  between \f$ K \f$ and \f$ \tilde{K} \f$.
118 
119  The symmetry energy, \f$ S(n_B,\delta), \f$ can be defined as
120  \f[
121  S(n_B,\delta) \equiv \frac{1}{2 n_B}\frac{\partial^2 \varepsilon}
122  {\partial \delta^2}
123  \f]
124  and the parameter \f$ S \f$ in Eq. 1 is just \f$ S(n_0,0) \f$.
125  Using
126  \f[
127  \left(\frac{\partial \varepsilon}{\partial \delta}\right)_{n_B} =
128  \frac{\partial \varepsilon}{\partial n_n}
129  \left(\frac{\partial n_n}{\partial \delta}\right)_{n_B} +
130  \frac{\partial \varepsilon}{\partial n_p}
131  \left(\frac{\partial n_p}{\partial \delta}\right)_{n_B}
132  = \frac{n_B}{2} \left(\mu_n - \mu_p \right)
133  \f]
134  this can be rewritten
135  \f[
136  S(n_B,\delta) = \frac{1}{4} \frac{\partial}{\partial \delta}
137  \left(\mu_n - \mu_p\right)
138  \f]
139  where the dependence of the chemical potentials on \f$ n_B \f$
140  and \f$ \delta \f$ is not written explicitly. This quantity is
141  computed by function \ref fesym(). Note that many of the
142  functions in this class are written in terms of the proton
143  fraction \f$ x_p = (1-\delta)/2 \f$ denoted as <tt>'pf'</tt>
144  instead of as functions of \f$ \delta \f$. Frequently, \f$
145  S(n_B,\delta) \f$ is evaluated at \f$ \delta=0 \f$ to give a
146  univariate function of the baryon density. It is sometimes also
147  evaluated at the point \f$ (n_B=n_0, \delta=0) \f$, and this
148  value is denoted by \f$ S \f$ above and is typically stored in
149  \ref esym.
150  Alternatively, one can define the symmetry energy by
151  \f[
152  \tilde{S}(n_B) = E(n_B,\delta=1)-E(n_B,\delta=0)
153  \f]
154  which is computed by function \ref fesym_diff() . The
155  functions \f$ S(n_B,\delta=0) \f$ and \f$ \tilde{S}(n_B) \f$
156  are equal when \f$ {\cal O}(\delta^4) \f$ terms are zero.
157  In this case, \f$ \mu_n - \mu_p \f$ is proportional to
158  \f$ \delta \f$ and so
159  \f[
160  S(n_B) = \tilde{S}(n_B) = \frac{1}{4}
161  \frac{(\mu_n-\mu_p)}{\delta} \, .
162  \f]
163 
164  These functions can also be generalized to finite temperature
165  \f[
166  S(n_B,\delta,T) = \frac{1}{4} \frac{\partial}{\partial \delta}
167  \left[\mu_n(n_B,\delta,T) - \mu_p(n_B,\delta,T)\right] \, ,
168  \f]
169  and
170  \f[
171  \tilde{S}(n_B,T) = F(n_B,\delta=1,T)-F(n_B,\delta=0,T) \, .
172  \f]
173 
174  The symmetry energy slope parameter \f$ L \f$, can be defined
175  by
176  \f[
177  L(n_B,\delta) \equiv 3 n_B \frac{\partial S(n_B,\delta)}
178  {\partial n_B} = 3 n_B \frac{\partial}{\partial n_B} \left[
179  \frac{1}{2 n_B} \frac{\partial^2 \varepsilon}{\partial \delta^2}
180  \right]
181  \f]
182  This can be rewritten as
183  \f[
184  L(n_B,\delta) = \frac{3 n_B}{4} \frac{\partial}{\partial n_B}
185  \frac{\partial}{\partial \delta} \left(\mu_n - \mu_p\right)
186  \f]
187  (where the derivatives can be evaluated in either order)
188  and this is the method used to compute this function
189  in \ref fesym_slope(). Alternatively, using
190  \f[
191  \left(\frac{\partial \varepsilon}{\partial n_B}\right)_{\delta} =
192  \frac{\partial \varepsilon}{\partial n_n}
193  \left(\frac{\partial n_n}{\partial n_B}\right)_{\delta} +
194  \frac{\partial \varepsilon}{\partial n_p}
195  \left(\frac{\partial n_p}{\partial n_B}\right)_{\delta}
196  = \frac{1}{2} \left(\mu_n + \mu_p \right)
197  \f]
198  \f$ L \f$ can be rewritten
199  \f{eqnarray*}
200  L(n_B,\delta) &=& 3 n_B \left[\frac{-1}{2 n_B^2}
201  \frac{\partial^2 \varepsilon}{\partial \delta^2} +
202  \frac{1}{4 n_B} \frac{\partial^2}{\partial \delta^2}
203  \left(\mu_n + \mu_p\right)\right] \\
204  &=& \frac{3}{4}\frac{\partial^2}{\partial \delta^2}
205  \left(\mu_n + \mu_p\right) - 3 S(n_B,\delta) \, .
206  \f}
207 
208  The third derivative with respect to the baryon density is
209  sometimes called the skewness. Here, we define
210  \f[
211  Q(n_B,\delta) = 27 n_B^3 \frac{\partial^3}{\partial n_B^3}
212  \left(\frac{\varepsilon}{n_B}\right) =
213  27 n_B^3 \frac{\partial^2}{\partial n_B^2}
214  \left(\frac{P}{n_B^2}\right)
215  \f]
216  and this function is computed in \ref fkprime() by computing
217  the second derivative of the pressure.
218 
219  The second derivative of the symmetry energy with respect
220  to the baryon density is
221  \f[
222  K_{\mathrm{sym}}(n_B,\delta) = 9 n_B^2
223  \frac{\partial^2}{\partial n_B^2} S(n_B,\delta)
224  \f]
225  and this function is computed in \ref fesym_curve().
226 
227  The third derivative of the symmetry energy with respect
228  to the baryon density is
229  \f[
230  Q_{\mathrm{sym}}(n_B,\delta) = 27 n_B^3
231  \frac{\partial^3}{\partial n_B^3} S(n_B,\delta)
232  \f]
233  and this function is computed in \ref fesym_skew(). Note that
234  the numerical evaluation of higher derivatives can make \ref
235  eos_had_base::fesym_curve() and \ref eos_had_base::fesym_skew()
236  inaccurate.
237 
238  Note that assuming terms of order \f$ \epsilon^3 \f$ and higher
239  are zero and solving for the baryon density for which \f$ P=0 \f$
240  gives, to order \f$ \delta^2 \f$ (\ref Piekarewicz09),
241  \f[
242  n_{B,\mathrm{sat}} = n_0 \left[ 1 - \frac{3 L \delta^2}{K} \right]
243  \f]
244  this implies a new 'incompressibility' around the saturation
245  point, i.e.
246  \f[
247  K(n_B=n_{B,\mathrm{sat}},\delta)=
248  K+\delta^2 \left( K_{\mathrm{sym}}-6 L- \frac{L Q}{K} \right)
249  + {\cal O}\left(\delta^4\right)
250  \f]
251  The quantity in parenthesis is referred to by some authors as
252  \f$ K_{\tau} \f$. Note that, because one is evaluating this at
253  \f$ n_B=n_{B,\mathrm{sat}} \f$, this is distinct from
254  \f[
255  \tilde{K}_{\tau} \equiv \frac{1}{2} \frac{\partial^2 K(n_B,\delta)}
256  {\partial \delta^2}
257  \f]
258  which is equal to \f$ K_{\mathrm{sym}} + 6 L \f$ to lowest order
259  in \f$ \delta \f$ at \f$ n_B = n_0 \f$.
260 
261  The quartic symmetry energy \f$ S_4(n_B,\delta) \f$ can be defined as
262  \f[
263  S_4(n_B,\delta) \equiv \frac{1}{24 n_B}\frac{\partial^4 \varepsilon}
264  {\partial \delta^4}
265  \f]
266  However, fourth derivatives are difficult numerically, and so an
267  alternative quantity is preferable. Instead, one can evaluate
268  the extent to which \f$ {\cal O}(\delta^4) \f$ terms are
269  important from
270  \f[
271  \eta(n_B) \equiv \frac{E(n_B,1)-E(n_B,1/2)}
272  {3 \left[E(n_B,1/2)-E(n_B,0)\right]}
273  \f]
274  as described in \ref Steiner06 . This function can be expressed
275  either in terms of \f$ \tilde{S} \f$ or \f$ S_4 \f$
276  \f[
277  \eta(n_B) = \frac{5 \tilde{S}(n_B) - S(n_B,0)}
278  {\tilde{S}(n_B) + 3 S(n_B,0)} =
279  \frac{5 S_4(n_B,0) + 4 S(n_B,0)}
280  {S_4(n_B,0) + 4 S(n_B,0)}
281  \f]
282  Alternatively, \f$ S_4 \f$ can be written
283  \f[
284  4 S(n_B) \left[ \frac{1-\eta(n_B)}{\eta(n_B)-5} \right] \, .
285  \f]
286 
287  Evaluating this function at the saturation density gives
288  \f[
289  \eta(n_0) = \frac{4 S + 5 S_4}{4 S + S_4}
290  \f]
291  (Note that \f$ S_4 \f$ is referred to as \f$ Q \f$ in
292  \ref Steiner06). Sometimes it is useful to separate out
293  the kinetic and potential parts of the energy density when
294  computing \f$ \eta \f$, and the class \ref eos_had_sym4_base
295  is useful for this purpose.
296 
297  The function \f$ L_4 \f$ can also be rewritten in
298  \f$ \eta^{\prime} \f$ (now suppressing the dependence
299  on \f$ n_B \f$),
300  \f[
301  \eta^{\prime} = \frac{16 \left( L_4 S - L S_4 \right)}
302  {3 n_B \left(4 S +S_4 \right)^2}
303  \f]
304  then using the expression for \f$ S_4 \f$,
305  \f[
306  \eta^{\prime} = \frac{\left(\eta -5\right)}{48 n_B S }
307  \left[ L_4 \left(\eta -5\right) + 4 L
308  \left(\eta -1\right)\right]
309  \f]
310 
311  \future Replace fmsom() with f_effm_scalar(). This has to wait
312  until f_effm_scalar() has a sensible definition when mn is
313  not equal to mp
314  \future Could write a function to compute the "symmetry free energy"
315  or the "symmetry entropy"
316  \future Compute the speed of sound or the number susceptibilities?
317  \future A lot of the numerical derivatives here might possibly
318  request negative number densities for the nucleons, which
319  may cause exceptions, espescially at very low densities.
320  Since the default EOS objects are GSL derivatives, one can
321  get around these issues by setting the GSL derivative object
322  step size, but this is a temporary solution.
323  */
324  class eos_had_base : public eos_base {
325 
326  public:
327 
329 
330  eos_had_base();
331 
332  virtual ~eos_had_base() {};
333 
334  /** \brief Binding energy (without the rest mass) in
335  \f$ \mathrm{fm}^{-1} \f$
336  */
337  double eoa;
338 
339  /// Compression modulus in \f$ \mathrm{fm}^{-1} \f$
340  double comp;
341 
342  /// Symmetry energy in \f$ \mathrm{fm}^{-1} \f$
343  double esym;
344 
345  /// Saturation density in \f$ \mathrm{fm}^{-3} \f$
346  double n0;
347 
348  /// Effective mass (neutron)
349  double msom;
350 
351  /// Skewness in \f$ \mathrm{fm}^{-1} \f$
352  double kprime;
353 
354  /** \brief If true, call the error handler if msolve() or
355  msolve_de() does not converge (default true)
356  */
358 
359  /// \name Equation of state
360  //@{
361  /** \brief Equation of state as a function of the chemical potentials
362  */
363  virtual int calc_p(fermion &n, fermion &p, thermo &th)=0;
364 
365  /** \brief Equation of state as a function of density
366  */
367  virtual int calc_e(fermion &n, fermion &p, thermo &th)=0;
368  //@}
369 
370  /// \name EOS properties
371  //@{
372  /** \brief Calculate the incompressibility in \f$ \mathrm{fm}^{-1} \f$
373  using calc_e()
374 
375  This function computes \f$ K (n_B,\delta) = 9 n_B \partial^2
376  \varepsilon /(\partial n_B^2) = 9 \partial P / (\partial n_B)
377  \f$ . The value of \f$ K(n_0,0) \f$, often referred to as the
378  "compressibility", is stored in \ref comp by \ref saturation()
379  and is about 240 MeV at saturation density.
380  */
381  virtual double fcomp(double nb, double delta=0.0);
382 
383  /** \brief Compute the incompressibility and its uncertainty
384 
385  This function works like \ref fcomp(), except that it also
386  returns the uncertainty in \c unc.
387  */
388  virtual double fcomp_err(double nb, double delta, double &unc);
389 
390  /** \brief Calculate the energy per baryon in \f$ \mathrm{fm}^{-1} \f$
391  using calc_e()
392 
393  This function computes the energy per baryon of matter
394  without the nucleon rest masses at the specified baryon
395  density, \c nb, and isospin asymmetry \c delta.
396  */
397  virtual double feoa(double nb, double delta=0.0);
398 
399  /** \brief Calculate symmetry energy of matter in
400  \f$ \mathrm{fm}^{-1} \f$ using \ref calc_dmu_delta() .
401 
402  This function computes the symmetry energy,
403  \f[
404  \left(\frac{1}{2 n_B}\frac{d^2 \varepsilon}{d \delta^2}
405  \right) = \frac{1}{4} \frac{\partial}{\partial \delta}
406  \left(\mu_n - \mu_p \right)
407  \f]
408  at the value of \f$ n_B \f$ given in \c nb and \f$ \delta \f$
409  given in \c delta. The symmetry energy at \f$ \delta=0 \f$ at
410  the saturation density and is stored in \ref esym by
411  \ref saturation().
412  */
413  virtual double fesym(double nb, double delta=0.0);
414 
415  /** \brief Calculate symmetry energy of matter and its
416  uncertainty in \f$ \mathrm{fm}^{-1} \f$
417 
418  This estimates the uncertainty due to the numerical
419  differentiation, assuming that difference betwen the neutron
420  and proton chemical potentials is computed exactly by \ref
421  calc_dmu_delta() .
422  */
423  virtual double fesym_err(double nb, double delta, double &unc);
424 
425  /** \brief The symmetry energy slope parameter
426  in \f$ \mathrm{fm}^{-1} \f$
427 
428  This returns the value of the "slope parameter" of the
429  symmetry energy as a function of baryon density \c nb and
430  isospin asymmetry \c delta. This ranges between about zero and
431  200 MeV for most equations of state.
432  */
433  virtual double fesym_slope(double nb, double delta=0.0);
434 
435  /** \brief The curvature of the symmetry energy
436  in \f$ \mathrm{fm}^{-1} \f$
437  */
438  virtual double fesym_curve(double nb, double delta=0.0);
439 
440  /** \brief The skewness of the symmetry energy
441  in \f$ \mathrm{fm}^{-1} \f$
442  */
443  virtual double fesym_skew(double nb, double delta=0.0);
444 
445  /** \brief Calculate symmetry energy of matter as energy of
446  neutron matter minus the energy of nuclear matter
447  in \f$ \mathrm{fm}^{-1} \f$
448 
449  This function returns the energy per baryon of neutron matter
450  minus the energy per baryon of nuclear matter. This will
451  deviate significantly from the results from fesym() only if
452  the dependence of the symmetry energy on \f$ \delta \f$ is not
453  quadratic.
454  */
455  virtual double fesym_diff(double nb);
456 
457  /** \brief The strength parameter for quartic terms in the
458  symmetry energy
459  */
460  virtual double feta(double nb);
461 
462  /** \brief The derivative of the strength parameter for quartic
463  terms in the symmetry energy
464  */
465  virtual double feta_prime(double nb);
466 
467  /** \brief Calculate skewness of nuclear matter
468  in \f$ \mathrm{fm}^{-1} \f$ using calc_e()
469 
470  The skewness is defined to be
471  \f$ 27 n^3 d^3 (\varepsilon/n)/(d n^3) =
472  27 n^3 d^2 (P/n^2) / (d n^2) \f$
473  and is denoted 'kprime'. This definition seems to be ambiguous
474  for densities other than the saturation density and is not
475  quite analogous to the compression modulus.
476  */
477  virtual double fkprime(double nb, double delta=0.0);
478 
479  /** \brief Calculate reduced neutron effective mass using calc_e()
480 
481  Neutron effective mass (as stored in <tt>part::ms</tt>)
482  divided by vacuum mass (as stored in <tt>part::m</tt>) in
483  nuclear matter at saturation density. Note that this simply
484  uses the value of n.ms from calc_e(), so that this effective
485  mass could be either the Landau or Dirac mass depending on the
486  context. Note that this may not be equal to the reduced proton
487  effective mass.
488  */
489  virtual double fmsom(double nb, double delta=0.0);
490 
491  /** \brief Neutron (reduced) effective mass
492  */
493  virtual double f_effm_neut(double nb, double delta=0.0);
494 
495  /** \brief Proton (reduced) effective mass
496  */
497  virtual double f_effm_prot(double nb, double delta=0.0);
498 
499  /** \brief Scalar effective mass
500 
501  Given the reduced nucleon effective masses, \f$ m_n^{*} \f$
502  and \f$ m_p^{*} \f$, the scalar and vector effective masses
503  are defined by (see e.g. \ref Farine01)
504  \f[
505  \frac{1}{m^{*}_n} = (1+\delta) \frac{1}{m^{*}_s} -
506  \delta \frac{1}{m^{*}_v}
507  \f]
508  \f[
509  \frac{1}{m^{*}_p} = (1-\delta) \frac{1}{m^{*}_s} +
510  \delta \frac{1}{m^{*}_v}
511  \f]
512  this implies
513  \f[
514  m_{\mathrm{scalar}}^{*} =
515  \frac{2 m^{*}_n m^{*}_p}{m^{*}_n+m^{*}_p}
516  \f]
517  and
518  \f[
519  m_{\mathrm{vector}}^{*} =
520  \frac{2 m^{*}_n m^{*}_p \delta}{m^{*}_n - m^{*}_p
521  + \delta(m^{*}_n + m^{*}_p)}
522  \f]
523  */
524  virtual double f_effm_scalar(double nb, double delta=0.0);
525 
526  /** \brief Vector effective mass
527 
528  See documentation for \ref eos_had_base::f_effm_scalar().
529 
530  Note that the vector effective mass diverges when \f$ m^{*}_n
531  = m^{*}_p \f$ and \f$ \delta = 0 \f$, but many models have
532  vector effective masses which are independent of \f$ \delta
533  \f$. For now, we set \f$ \delta =1 \f$ to be the default
534  value, corresponding to neutron matter.
535  */
536  virtual double f_effm_vector(double nb, double delta=1.0);
537 
538  /** \brief Calculate saturation density using calc_e()
539 
540  This function finds the baryon density for which the pressure
541  vanishes.
542  */
543  virtual double fn0(double delta, double &leoa);
544 
545  /** \brief Compute the number susceptibilities as a function of
546  the chemical potentials, \f$ \partial^2 P / \partial \mu_i
547  \mu_j \f$
548  */
549  virtual void f_number_suscept(double mun, double mup, double &dPdnn,
550  double &dPdnp, double &dPdpp);
551 
552  /** \brief Compute the 'inverse' number susceptibilities as a
553  function of the densities, \f$ \partial^2 \varepsilon /
554  \partial n_i n_j \f$
555  */
556  virtual void f_inv_number_suscept(double mun, double mup, double &dednn,
557  double &dednp, double &dedpp);
558 
559  /** \brief Calculates some of the EOS properties at the saturation
560  density
561 
562  This computes the saturation density, and the
563  incompressibility, the symmetry energy, the binding energy,
564  the reduced neutron effective mass at the saturation density,
565  and the skewness in isospin-symmetric matter. The results are
566  stored in \ref n0, \ref comp, \ref esym, \ref eoa, \ref msom,
567  and \ref kprime, respectively.
568 
569  \future It would be great to provide numerical uncertainties
570  in the saturation properties.
571  */
572  virtual int saturation();
573  //@}
574 
575  /// \name Functions for calculating physical properties
576  //@{
577  /** \brief Compute the neutron chemical potential at fixed
578  density
579 
580  This function uses \ref neutron, \ref proton, \ref
581  eos_base::eos_thermo, and \ref calc_e() .
582  */
583  double calc_mun_e(double nn, double np);
584 
585  /** \brief Compute the energy density as a function of the nucleon
586  densities
587  */
588  double calc_ed(double nn, double np);
589 
590  /** \brief Compute the pressure as a function of the nucleon
591  chemical potentials
592  */
593  double calc_pr(double nn, double np);
594 
595  /** \brief Compute the proton chemical potential at fixed
596  density
597 
598  This function uses \ref neutron, \ref proton, \ref
599  eos_base::eos_thermo, and \ref calc_e() .
600  */
601  double calc_mup_e(double nn, double np);
602 
603  /** \brief Compute the neutron density at fixed
604  chemical potential
605 
606  This function uses \ref neutron, \ref proton, \ref
607  eos_base::eos_thermo, and \ref calc_e() .
608  */
609  double calc_nn_p(double mun, double mup);
610 
611  /** \brief Compute the proton density at fixed
612  chemical potential
613 
614  This function uses \ref neutron, \ref proton, \ref
615  eos_base::eos_thermo, and \ref calc_e() .
616  */
617  double calc_np_p(double mun, double mup);
618 
619  /** \brief Compute the difference between neutron and proton chemical
620  potentials as a function of the isospin asymmetry
621 
622  This function uses \ref neutron, \ref proton, \ref
623  eos_base::eos_thermo, and \ref calc_e() .
624  */
625  double calc_dmu_delta(double delta, double nb);
626 
627  /** \brief Compute the sum of the neutron and proton chemical
628  potentials as a function of the isospin asymmetry
629 
630  This uses \ref neutron, \ref proton, \ref eos_base::eos_thermo,
631  and \ref calc_e() .
632  */
633  double calc_musum_delta(double delta, double nb);
634 
635  /** \brief Compute the pressure as a function of baryon density
636  at fixed isospin asymmetry
637 
638  Used by fcomp().
639  */
640  double calc_pressure_nb(double nb, double delta=0.0);
641 
642  /** \brief Compute the energy density as a function of baryon density
643  at fixed isospin asymmetry
644 
645  This uses \ref neutron, \ref proton, \ref eos_base::eos_thermo,
646  and \ref calc_e() .
647  */
648  double calc_edensity_nb(double nb, double delta=0.0);
649 
650  /** \brief Compute derivatives at constant proton fraction
651  */
652  void const_pf_derivs(double nb, double pf,
653  double &dednb_pf, double &dPdnb_pf);
654 
655  /** \brief Calculate pressure / baryon density squared in nuclear
656  matter as a function of baryon density at fixed isospin asymmetry
657 
658  Used by fkprime().
659 
660  This uses \ref neutron, \ref proton, \ref eos_base::eos_thermo,
661  and \ref calc_e() .
662  */
663  double calc_press_over_den2(double nb, double delta=0.0);
664 
665  /** \brief Calculate energy density as a function of the
666  isospin asymmetry at fixed baryon density
667 
668  Used by fesym().
669 
670  This function calls \ref eos_had_base::calc_e() with the
671  internally stored neutron and proton objects.
672  */
673  double calc_edensity_delta(double delta, double nb);
674  //@}
675 
676  /// \name Nuclear matter functions
677  //@{
678  /** \brief Solve for the chemical potentials given the densities
679 
680  The neutron and proton chemical potentials should be stored in
681  <tt>x[0]</tt> and <tt>x[1]</tt> and the neutron and proton
682  densities should be stored in <tt>pa[0]</tt> and
683  <tt>pa[1]</tt>.
684 
685  Because this function is designed to be used in a solver,
686  it returns <tt>exc_efailed</tt> without calling the error
687  handler if the densities are not finite.
688 
689  This function is used by \ref eos_had_pres_base::calc_e().
690  */
691  int nuc_matter_p(size_t nv, const ubvector &x, ubvector &y,
692  double nn0, double np0);
693 
694  /** \brief Solve for the densities given the chemical potentials
695 
696  The neutron and proton densities should be stored in
697  <tt>x[0]</tt> and <tt>x[1]</tt> and the neutron and proton
698  chemical potentials should be stored in <tt>pa[0]</tt> and
699  <tt>pa[1]</tt>.
700 
701  Because this function is designed to be used in a solver,
702  it returns <tt>exc_efailed</tt> without calling the error
703  handler if the chemical potentials are not finite.
704 
705  This function is used by \ref eos_had_eden_base::calc_p().
706  */
707  int nuc_matter_e(size_t nv, const ubvector &x, ubvector &y,
708  double mun0, double mup0);
709  //@}
710 
711  /// \name Set auxilliary objects
712  //@{
713  /** \brief Set class mroot object for use in calculating chemical
714  potentials from densities
715 
716  \note While in principle this allows one to use any \ref mroot
717  object, in practice some of the current EOSs require \ref
718  mroot_hybrids because it automatically avoids regions
719  where the equations are undefined.
720  */
721  virtual void set_mroot(mroot<> &mr);
722 
723  /** \brief Set class mroot object for use calculating saturation density
724 
725  \note While in principle this allows one to use any \ref mroot
726  object, in practice some of the current EOSs require \ref
727  mroot_hybrids because it automatically avoids regions
728  where the equations are undefined.
729  */
730  virtual void set_sat_root(root<> &mr);
731 
732  /// Set \ref deriv_base object to use to find saturation properties
733  virtual void set_sat_deriv(deriv_base<> &de);
734 
735  /** \brief Set the second \ref deriv_base object to use to find
736  saturation properties
737 
738  Computing the slope of the symmetry energy at the saturation
739  density requires two derivative objects, because it has to
740  take an isospin derivative and a density derivative. Thus this
741  second \ref deriv_base object is used in the function
742  fesym_slope().
743  */
744  virtual void set_sat_deriv2(deriv_base<> &de);
745 
746  /// Set neutron and proton
747  virtual void set_n_and_p(fermion &n, fermion &p);
748  //@}
749 
750  /** \brief The defaut neutron
751 
752  By default this has a spin degeneracy of 2 and a mass of \ref
753  o2scl_mks::mass_neutron . Also the value of
754  <tt>part::non_interacting</tt> is set to <tt>false</tt>.
755  */
757 
758  /** \brief The defaut proton
759 
760  By default this has a spin degeneracy of 2 and a mass of \ref
761  o2scl_mks::mass_proton . Also the value of
762  <tt>part::non_interacting</tt> is set to <tt>false</tt>.
763  */
765 
766  /// \name Default solvers and derivative classes
767  //@{
768  /** \brief The default object for derivatives
769 
770  The value of deriv_gsl::h is set to \f$ 10^{-3} \f$ in
771  the eos_had_base constructor.
772  */
774 
775  /** \brief The second default object for derivatives
776 
777  The value of deriv_gsl::h is set to \f$ 10^{-3} \f$ in
778  the eos_had_base constructor.
779  */
781 
782  /** \brief The default solver
783 
784  Used by calc_e() to solve nuc_matter_p() (2 variables) and by
785  calc_p() to solve nuc_matter_e() (2 variables).
786  */
788 
789  /** \brief The default solver for calculating the saturation
790  density
791 
792  Used by fn0() (which is called by saturation()) to solve
793  saturation_matter_e() (1 variable).
794  */
796  //@}
797 
798  /// \name Other functions
799  //@{
800  /** \brief Calculate coefficients for \gradient \part of Hamiltonian
801 
802  \note This is still somewhat experimental.
803 
804  We want the \gradient \part of the Hamiltonian in the form
805  \f[
806  {\cal H}_{\mathrm{grad}} = \frac{1}{2} \sum_{i=n,p}
807  \sum_{j=n,p} Q_{ij}
808  \vec{\nabla} n_i \cdot \vec{\nabla} n_j
809  \f]
810 
811  The expression for the \gradient terms from \ref Pethick95 is
812  \f{eqnarray*}
813  {\cal H}_{\mathrm{grad}}&=&-\frac{1}{4}
814  \left(2 P_1 + P_{1;f}-P_{2;f}\right)
815  \nonumber \\
816  && +\frac{1}{2} \left( Q_1+Q_2 \right)
817  \left(n_n \nabla^2 n_n+n_p \nabla^2 n_p\right) \nonumber \\
818  && + \frac{1}{4}\left( Q_1-Q_2 \right)
819  \left[\left(\nabla n_n\right)^2+\left(\nabla n_p\right)^2
820  \right] \nonumber \\
821  && + \frac{1}{2} \frac{d Q_2}{d n}
822  \left( n_n \nabla n_n + n_p \nabla n_p \right) \cdot \nabla n
823  \f}
824  This can be rewritten
825  \f{eqnarray*}
826  {\cal H}_{\mathrm{grad}}&=&\frac{1}{2}\left(\nabla n\right)^2
827  \left[ \frac{3}{2} P_1+n \frac{d P_1}{d n}\right] \nonumber \\
828  && - \frac{3}{4} \left[ \left( \nabla n_n\right)^2 +
829  \left( \nabla n_p \right)^2 \right] \nonumber \\
830  && -\frac{1}{2} \left[ \right] \cdot \nabla n \frac{d Q_1}{d n}
831  \nonumber \\ && - \frac{1}{4} \left( \nabla n\right)^2 P_2
832  - \frac{1}{4} \left[ \left( \nabla n_n\right)^2 +
833  \left( \nabla n_p \right)^2 \right] Q_2
834  \f}
835  or
836  \f{eqnarray*}
837  {\cal H}_{\mathrm{grad}}&=&\frac{1}{4} \left( \nabla n\right)^2
838  \left[3 P_1 + 2 n \frac{d P_1}{d n}-P_2\right] \nonumber \\
839  && - \frac{1}{4} \left( 3 Q_1+Q_2 \right)
840  \left[ \left( \nabla n_n\right)^2 +
841  \left( \nabla n_p \right)^2 \right] \nonumber \\
842  && - \frac{1}{2} \frac{d Q_1}{d n}
843  \left[ n_n \nabla n_n + n_p \nabla n_p \right]
844  \cdot \nabla n
845  \f}
846  or
847  \f{eqnarray*}
848  {\cal H}_{\mathrm{grad}}&=&\frac{1}{4} \left( \nabla n\right)^2
849  \left[3 P_1 + 2 n \frac{d P_1}{d n}-P_2\right] \nonumber \\
850  && - \frac{1}{4} \left( 3 Q_1+Q_2 \right)
851  \left[ \left( \nabla n_n\right)^2 +
852  \left( \nabla n_p \right)^2 \right] \nonumber \\
853  && - \frac{1}{2} \frac{d Q_1}{d n}
854  \left[ n_n \left( \nabla n_n \right)^2 +
855  n_p \left( \nabla n_p \right)^2 + n \nabla n_n \cdot
856  \nabla n_p \right]
857  \f}
858 
859  Generally, for Skyrme-like interactions
860  \f{eqnarray*}
861  P_i &=& \frac{1}{4} t_i \left(1+\frac{1}{2} x_i \right) \nonumber \\
862  Q_i &=& \frac{1}{4} t_i \left(\frac{1}{2} + x_i \right) \, .
863  \f}
864  for \f$ i=1,2 \f$ .
865 
866  This function uses the assumption \f$ x_1=x_2=0 \f$ to
867  calculate \f$ t_1 \f$ and \f$ t_2 \f$ from the neutron
868  and proton effective masses assuming the Skyrme form. The
869  values of \f$ Q_{ij} \f$ and their derivatives are then computed.
870 
871  The functions set_n_and_p() and set_thermo() will be called by
872  gradient_qij(), to facilitate the use of the \c n, \c p, and
873  \c th parameters.
874 
875  */
876  void gradient_qij(fermion &n, fermion &p, thermo &th,
877  double &qnn, double &qnp, double &qpp,
878  double &dqnndnn, double &dqnndnp,
879  double &dqnpdnn, double &dqnpdnp,
880  double &dqppdnn, double &dqppdnp);
881 
882  /// Return string denoting type ("eos_had_base")
883  virtual const char *type() { return "eos_had_base"; }
884  //@}
885 
886  /// \name Consistency checks
887  //@{
888  /** \brief Check the chemical potentials by computing the
889  derivatives numerically
890  */
891  void check_mu(fermion &n, fermion &p, thermo &th,
892  double &mun_deriv,
893  double &mup_deriv,
894  double &mun_err, double &mup_err);
895 
896  /** \brief Check the densities by computing the
897  derivatives numerically
898  */
899  void check_den(fermion &n, fermion &p, thermo &th,
900  double &nn_deriv, double &np_deriv,
901  double &nn_err, double &np_err);
902  //@}
903 
904 #ifndef DOXYGEN_INTERNAL
905 
906  protected:
907 
908  /// Compute t1 for gradient_qij().
909  double t1_fun(double barn);
910 
911  /// Compute t2 for gradient_qij().
912  double t2_fun(double barn);
913 
914  /// The EOS solver
916 
917  /// The solver to compute saturation properties
919 
920  /// The derivative object for saturation properties
922 
923  /// The second derivative object for saturation properties
925 
926  /// The neutron object
928 
929  /// The proton object
931 
932 #endif
933 
934  };
935 
936  /// A hadronic EOS based on a function of the densities [abstract base]
938  public:
939 
940  /** \brief Equation of state as a function of density
941  */
942  virtual int calc_e(fermion &n, fermion &p, thermo &th)=0;
943 
944  /** \brief Equation of state as a function of the chemical potentials
945  */
946  virtual int calc_p(fermion &n, fermion &p, thermo &th);
947 
948  };
949 
950  /** \brief A hadronic EOS based on a function of the chemical
951  potentials [abstract base]
952  */
954  public:
955 
956  /** \brief Equation of state as a function of the chemical potentials
957  */
958  virtual int calc_p(fermion &n, fermion &p, thermo &th)=0;
959 
960  /** \brief Equation of state as a function of density
961  */
962  virtual int calc_e(fermion &n, fermion &p, thermo &th);
963 
964  };
965 
966  /// A finite temperature hadronic EOS [abstract base]
968 
969 #ifndef DOXYGEN_INTERNAL
970 
971  protected:
972 
973  /// Fermion thermodynamics (default is \ref def_fet)
975 
976  /// Solve for nuclear matter at finite temperature given density
977  int nuc_matter_temp_e(size_t nv, const ubvector &x,
978  ubvector &y, double nn0, double np0, double T);
979 
980  /// Solve for nuclear matter at finite temperature given mu
981  int nuc_matter_temp_p(size_t nv, const ubvector &x,
982  ubvector &y, double mun0, double mup0, double T);
983 
984  /** \brief Solve for the liquid gas phase transition as
985  a function of the densities
986  */
987  int liqgas_dens_solve(size_t nv, const ubvector &x,
988  ubvector &y, fermion &n1, fermion &p1,
989  fermion &n2, fermion &p2, double T,
990  thermo &th1, thermo &th2);
991 
992  /** \brief Solve for the liquid-gas phase transition at fixed
993  baryon density and electron fraction
994  */
995  int liqgas_solve(size_t nv, const ubvector &x,
996  ubvector &y, fermion &n1, fermion &p1,
997  fermion &n2, fermion &p2, double nB0,
998  double Ye0, double T,
999  thermo &th1, thermo &th2);
1000 
1001  /** \brief Solve for the liquid-gas phase transition in
1002  beta-equilibrium
1003  */
1004  int liqgas_beta_solve(size_t nv, const ubvector &x,
1005  ubvector &y, fermion &n1, fermion &p1,
1006  fermion &n2, fermion &p2,
1007  double nB0, double T,
1008  thermo &th1, thermo &th2, fermion &e);
1009 
1010  /** \brief Compute the entropy
1011  */
1012  double calc_entropy_delta(double delta, double nb, double T);
1013 
1014  /** \brief Compute the difference between the neutron and
1015  proton chemical potentials
1016  */
1017  double calc_dmu_delta_T(double delta, double nb, double T);
1018 
1019 #endif
1020 
1021  public:
1022 
1023  eos_had_temp_base() {
1024  fet=&def_fet;
1025  }
1026 
1027  virtual ~eos_had_temp_base() {}
1028 
1029  /// \name Basic usage
1030  //@{
1031  /** \brief Equation of state as a function of density
1032  */
1033  virtual int calc_e(fermion &n, fermion &p, thermo &th)=0;
1034 
1035  /** \brief Equation of state as a function of densities at
1036  finite temperature
1037  */
1038  virtual int calc_temp_e(fermion &n, fermion &p, double T,
1039  thermo &th)=0;
1040 
1041  /** \brief Equation of state as a function of the chemical potentials
1042  */
1043  virtual int calc_p(fermion &n, fermion &p, thermo &th)=0;
1044 
1045  /** \brief Equation of state as a function of the chemical potentials
1046  at finite temperature
1047  */
1048  virtual int calc_temp_p(fermion &n, fermion &p, double T,
1049  thermo &th)=0;
1050  //@}
1051 
1052  /// Computing finite-temperature integrals
1053  //@{
1054  /** \brief Set the object for computing finite-temperature fermions
1055  (default is \ref def_fet)
1056  */
1058  fet=&f;
1059  }
1060 
1061  /// Default fermion thermodynamics object
1063  //@}
1064 
1065  /// \name Liquid-gas transition functions
1066  //@{
1067  /** \brief Compute liquid-gas phase transition densities using
1068  \ref eos_had_temp_base::calc_temp_e() .
1069 
1070  At fixed baryon number density for \c n1, this determines the
1071  baryon number densities for \c p1, \c n2, and \c p2 which give
1072  chemical and mechanical equilibrium at a fixed temperature
1073  \c T. The thermodynamic quantities assuming bulk matter for
1074  each set is stored in \c th1 and \c th2.
1075  */
1076  virtual int calc_liqgas_dens_temp_e
1077  (fermion &n1, fermion &p1, fermion &n2, fermion &p2,
1078  double T, thermo &th1, thermo &th2);
1079 
1080  /** \brief Compute the liquid-gas phase transition using
1081  \ref eos_had_temp_base::calc_temp_e() .
1082 
1083  At fixed baryon number density, \c nB, fixed electron
1084  fraction, \c Ye, and fixed temperature, \c T, this function
1085  determines the baryon number densities for \c n1, \c p1, \c
1086  n2, and \c p2 which give chemical and mechanical equilibrium.
1087  The thermodynamic quantities assuming bulk matter for each set
1088  is stored in \c th1 and \c th2, and the volume fraction of
1089  phase 1 is stored in \c chi.
1090  */
1091  virtual int calc_liqgas_temp_e
1092  (fermion &n1, fermion &p1, fermion &n2, fermion &p2,
1093  double nB, double Ye, double T, thermo &th1, thermo &th2,
1094  double &chi);
1095 
1096  /** \brief Compute the liquid-gas phase transition in
1097  beta-equilibrium using \ref eos_had_temp_base::calc_temp_e() .
1098 
1099  At fixed baryon number density, \c nB, and fixed temperature,
1100  \c T, this function determines the baryon number densities for
1101  \c n1, \c p1, \c n2, and \c p2 which give chemical and
1102  mechanical equilibrium assuming beta-equilibrium with
1103  electrons. The thermodynamic quantities assuming bulk matter
1104  for each set is stored in \c th1 and \c th2, the electron
1105  fraction is stored in \c Ye, and the volume fraction of phase
1106  1 is stored in \c chi.
1107  */
1108  virtual int calc_liqgas_beta_temp_e
1109  (fermion &n1, fermion &p1, fermion &n2, fermion &p2,
1110  double nB, double T, thermo &th1, thermo &th2,
1111  double &Ye, double &chi);
1112  //@}
1113 
1114  /// \name Functions related to the symmetry energy
1115  //@{
1116  /** \brief Compute the symmetry energy at finite temperature
1117  */
1118  virtual double fesym_T(double nb, double T, double delta=0.0);
1119 
1120  /** \brief Compute the symmetry entropy at finite temperature
1121  */
1122  virtual double fsyment_T(double nb, double T, double delta=0.0);
1123  //@}
1124 
1125  /// \name Helper functions
1126  //@{
1127  /// Neutron chemical potential as a function of the densities
1128  virtual double calc_temp_mun_e(double nn, double np, double T);
1129  /// Proton chemical potential as a function of the densities
1130  virtual double calc_temp_mup_e(double nn, double np, double T);
1131  /// Neutron density as a function of the chemical potentials
1132  virtual double calc_temp_nn_p(double mun, double mup, double T);
1133  /// Proton density as a function of the chemical potentials
1134  virtual double calc_temp_np_p(double mun, double mup, double T);
1135 
1136  /** \brief Compute the free energy as a function of the temperature
1137  and the densities
1138  */
1139  double calc_fr(double nn, double np, double T);
1140  //@}
1141 
1142  /// \name Susceptibilities
1143  //@{
1144  /** \brief Compute the number susceptibilities as a function of
1145  the chemical potentials, \f$ \partial^2 P / \partial \mu_i
1146  \mu_j \f$ at a fixed temperature
1147  */
1148  virtual void f_number_suscept_T
1149  (double mun, double mup, double T, double &dPdnn,
1150  double &dPdnp, double &dPdpp);
1151 
1152  /** \brief Compute the 'inverse' number susceptibilities as a
1153  function of the densities, \f$ \partial^2 \varepsilon /
1154  \partial n_i n_j \f$ at a fixed temperature
1155  */
1156  virtual void f_inv_number_suscept_T
1157  (double mun, double mup, double T, double &dednn,
1158  double &dednp, double &dedpp);
1159  //@}
1160 
1161  /// \name Consistency check
1162  //@{
1163  /** \brief Check the entropy by computing the
1164  derivative numerically
1165  */
1166  void check_en(fermion &n, fermion &p, double T, thermo &th,
1167  double &en_deriv, double &en_err);
1168 
1169  /** \brief Check the chemical potentials at finite temperature
1170  by computing the derivative numerically
1171  */
1172  void check_mu_T(fermion &n, fermion &p, double T, thermo &th,
1173  double &mun_deriv, double &mup_deriv,
1174  double &mun_err, double &mup_err);
1175  //@}
1176 
1177  };
1178 
1179  /** \brief A hadronic EOS at finite temperature
1180  based on a function of the densities [abstract base]
1181 
1182  This class provides automatic computation of \ref calc_e() and
1183  \ref calc_temp_e() assuming that \ref calc_p() and \ref
1184  calc_temp_p() are specified.
1185  */
1187  public:
1188 
1189  /** \brief Equation of state as a function of density
1190  */
1191  virtual int calc_e(fermion &n, fermion &p, thermo &th)=0;
1192 
1193  /** \brief Equation of state as a function of densities at
1194  finite temperature
1195  */
1196  virtual int calc_temp_e(fermion &n, fermion &p, double T,
1197  thermo &th)=0;
1198 
1199  /** \brief Equation of state as a function of the chemical potentials
1200  */
1201  virtual int calc_p(fermion &n, fermion &p, thermo &th);
1202 
1203  /** \brief Equation of state as a function of the chemical potentials
1204  at finite temperature
1205  */
1206  virtual int calc_temp_p(fermion &n, fermion &p, double T,
1207  thermo &th);
1208 
1209  };
1210 
1211  /** \brief A hadronic EOS at finite temperature based on a function
1212  of the chemical potentials [abstract base]
1213 
1214  This class provides automatic computation of \ref calc_p() and
1215  \ref calc_temp_p() assuming that \ref calc_e() and \ref
1216  calc_temp_e() are specified.
1217  */
1219  public:
1220 
1221  /** \brief Equation of state as a function of the chemical potentials
1222  */
1223  virtual int calc_p(fermion &n, fermion &p, thermo &th)=0;
1224 
1225  /** \brief Equation of state as a function of the chemical potentials
1226  at finite temperature
1227  */
1228  virtual int calc_temp_p(fermion &n, fermion &p, double T,
1229  thermo &th)=0;
1230 
1231  /** \brief Equation of state as a function of density
1232  */
1233  virtual int calc_e(fermion &n, fermion &p, thermo &th);
1234 
1235  /** \brief Equation of state as a function of densities at
1236  finite temperature
1237  */
1238  virtual int calc_temp_e(fermion &n, fermion &p, double T,
1239  thermo &th);
1240  };
1241 
1242 #ifndef DOXYGEN_NO_O2NS
1243 }
1244 #endif
1245 
1246 #endif
1247 
root_cern def_sat_root
The default solver for calculating the saturation density.
Definition: eos_had_base.h:795
virtual double fesym_skew(double nb, double delta=0.0)
The skewness of the symmetry energy in .
virtual void set_sat_deriv(deriv_base<> &de)
Set deriv_base object to use to find saturation properties.
virtual void set_fermion_eval_thermo(fermion_eval_thermo &f)
Computing finite-temperature integrals.
virtual double feta_prime(double nb)
The derivative of the strength parameter for quartic terms in the symmetry energy.
double calc_ed(double nn, double np)
Compute the energy density as a function of the nucleon densities.
fermion_eval_thermo * fet
Fermion thermodynamics (default is def_fet)
Definition: eos_had_base.h:974
Equation of state base class.
Definition: eos_base.h:39
double comp
Compression modulus in .
Definition: eos_had_base.h:340
double calc_np_p(double mun, double mup)
Compute the proton density at fixed chemical potential.
virtual double fkprime(double nb, double delta=0.0)
Calculate skewness of nuclear matter in using calc_e()
deriv_gsl def_deriv
The default object for derivatives.
Definition: eos_had_base.h:773
const double barn
virtual double f_effm_neut(double nb, double delta=0.0)
Neutron (reduced) effective mass.
double t1_fun(double barn)
Compute t1 for gradient_qij().
void const_pf_derivs(double nb, double pf, double &dednb_pf, double &dPdnb_pf)
Compute derivatives at constant proton fraction.
void check_mu(fermion &n, fermion &p, thermo &th, double &mun_deriv, double &mup_deriv, double &mun_err, double &mup_err)
Check the chemical potentials by computing the derivatives numerically.
virtual double fmsom(double nb, double delta=0.0)
Calculate reduced neutron effective mass using calc_e()
double calc_mun_e(double nn, double np)
Compute the neutron chemical potential at fixed density.
double calc_pressure_nb(double nb, double delta=0.0)
Compute the pressure as a function of baryon density at fixed isospin asymmetry.
mroot * eos_mroot
The EOS solver.
Definition: eos_had_base.h:915
virtual double fcomp(double nb, double delta=0.0)
Calculate the incompressibility in using calc_e()
virtual double f_effm_scalar(double nb, double delta=0.0)
Scalar effective mass.
virtual double fn0(double delta, double &leoa)
Calculate saturation density using calc_e()
double esym
Symmetry energy in .
Definition: eos_had_base.h:343
virtual void f_inv_number_suscept(double mun, double mup, double &dednn, double &dednp, double &dedpp)
Compute the &#39;inverse&#39; number susceptibilities as a function of the densities, .
virtual double fesym_diff(double nb)
Calculate symmetry energy of matter as energy of neutron matter minus the energy of nuclear matter in...
void check_den(fermion &n, fermion &p, thermo &th, double &nn_deriv, double &np_deriv, double &nn_err, double &np_err)
Check the densities by computing the derivatives numerically.
A hadronic EOS at finite temperature based on a function of the densities [abstract base]...
A finite temperature hadronic EOS [abstract base].
Definition: eos_had_base.h:967
virtual const char * type()
Return string denoting type ("eos_had_base")
Definition: eos_had_base.h:883
deriv_base * sat_deriv
The derivative object for saturation properties.
Definition: eos_had_base.h:921
double n0
Saturation density in .
Definition: eos_had_base.h:346
virtual int saturation()
Calculates some of the EOS properties at the saturation density.
double calc_press_over_den2(double nb, double delta=0.0)
Calculate pressure / baryon density squared in nuclear matter as a function of baryon density at fixe...
double calc_nn_p(double mun, double mup)
Compute the neutron density at fixed chemical potential.
root * sat_root
The solver to compute saturation properties.
Definition: eos_had_base.h:918
virtual int calc_p(fermion &n, fermion &p, thermo &th)=0
Equation of state as a function of the chemical potentials.
virtual void set_sat_root(root<> &mr)
Set class mroot object for use calculating saturation density.
A hadronic EOS at finite temperature based on a function of the chemical potentials [abstract base]...
void gradient_qij(fermion &n, fermion &p, thermo &th, double &qnn, double &qnp, double &qpp, double &dqnndnn, double &dqnndnp, double &dqnpdnn, double &dqnpdnp, double &dqppdnn, double &dqppdnp)
Calculate coefficients for gradient gradient part part of Hamiltonian.
double t2_fun(double barn)
Compute t2 for gradient_qij().
double calc_mup_e(double nn, double np)
Compute the proton chemical potential at fixed density.
fermion * proton
The proton object.
Definition: eos_had_base.h:930
deriv_gsl def_deriv2
The second default object for derivatives.
Definition: eos_had_base.h:780
double calc_edensity_delta(double delta, double nb)
Calculate energy density as a function of the isospin asymmetry at fixed baryon density.
A hadronic EOS based on a function of the densities [abstract base].
Definition: eos_had_base.h:937
virtual double f_effm_vector(double nb, double delta=1.0)
Vector effective mass.
double calc_musum_delta(double delta, double nb)
Compute the sum of the neutron and proton chemical potentials as a function of the isospin asymmetry...
deriv_base * sat_deriv2
The second derivative object for saturation properties.
Definition: eos_had_base.h:924
int nuc_matter_e(size_t nv, const ubvector &x, ubvector &y, double mun0, double mup0)
Solve for the densities given the chemical potentials.
Hadronic equation of state [abstract base].
Definition: eos_had_base.h:324
A hadronic EOS based on a function of the chemical potentials [abstract base].
Definition: eos_had_base.h:953
virtual double fcomp_err(double nb, double delta, double &unc)
Compute the incompressibility and its uncertainty.
double calc_edensity_nb(double nb, double delta=0.0)
Compute the energy density as a function of baryon density at fixed isospin asymmetry.
virtual int calc_e(fermion &n, fermion &p, thermo &th)=0
Equation of state as a function of density.
virtual void f_number_suscept(double mun, double mup, double &dPdnn, double &dPdnp, double &dPdpp)
Compute the number susceptibilities as a function of the chemical potentials, .
fermion def_proton
The defaut proton.
Definition: eos_had_base.h:764
virtual double feta(double nb)
The strength parameter for quartic terms in the symmetry energy.
virtual double fesym(double nb, double delta=0.0)
Calculate symmetry energy of matter in using calc_dmu_delta() .
double kprime
Skewness in .
Definition: eos_had_base.h:352
virtual double fesym_err(double nb, double delta, double &unc)
Calculate symmetry energy of matter and its uncertainty in .
bool err_nonconv
If true, call the error handler if msolve() or msolve_de() does not converge (default true) ...
Definition: eos_had_base.h:357
virtual double f_effm_prot(double nb, double delta=0.0)
Proton (reduced) effective mass.
virtual double fesym_slope(double nb, double delta=0.0)
The symmetry energy slope parameter in .
int nuc_matter_p(size_t nv, const ubvector &x, ubvector &y, double nn0, double np0)
Solve for the chemical potentials given the densities.
double calc_pr(double nn, double np)
Compute the pressure as a function of the nucleon chemical potentials.
fermion * neutron
The neutron object.
Definition: eos_had_base.h:927
mroot_hybrids def_mroot
The default solver.
Definition: eos_had_base.h:787
virtual double feoa(double nb, double delta=0.0)
Calculate the energy per baryon in using calc_e()
double msom
Effective mass (neutron)
Definition: eos_had_base.h:349
double calc_dmu_delta(double delta, double nb)
Compute the difference between neutron and proton chemical potentials as a function of the isospin as...
double eoa
Binding energy (without the rest mass) in .
Definition: eos_had_base.h:332
fermion_eff def_fet
Default fermion thermodynamics object.
virtual double fesym_curve(double nb, double delta=0.0)
The curvature of the symmetry energy in .
virtual void set_sat_deriv2(deriv_base<> &de)
Set the second deriv_base object to use to find saturation properties.
virtual void set_mroot(mroot<> &mr)
Set class mroot object for use in calculating chemical potentials from densities. ...
virtual void set_n_and_p(fermion &n, fermion &p)
Set neutron and proton.
fermion def_neutron
The defaut neutron.
Definition: eos_had_base.h:756

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