part_deriv.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2020, 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 #ifndef O2SCL_PART_DERIV_H
24 #define O2SCL_PART_DERIV_H
25 
26 /** \file part_deriv.h
27  \brief File defining \ref o2scl::part_deriv_tl
28 */
29 
30 #include <string>
31 #include <iostream>
32 #include <fstream>
33 #include <cmath>
34 #include <o2scl/part.h>
35 #include <o2scl/fermion.h>
36 #include <o2scl/fermion_rel.h>
37 
38 #ifndef DOXYGEN_NO_O2NS
39 namespace o2scl {
40 #endif
41 
42  /** \brief Object to store second derivatives of
43  \f$ P(\mu_n,\mu_p,T) \f$
44  */
45  template<class fp_t=double>
47 
48  public:
49 
50  /// The quantity \f$ (\partial^2 P)/(\partial T^2) \f$
51  fp_t dsdT;
52  /// The quantity \f$ (\partial^2 P)/(\partial T \partial \mu_n) \f$
53  fp_t dnndT;
54  /// The quantity \f$ (\partial^2 P)/(\partial T \partial \mu_p) \f$
55  fp_t dnpdT;
56  /// The quantity \f$ (\partial^2 P)/(\partial \mu_n^2) \f$
57  fp_t dnndmun;
58  /// The quantity \f$ (\partial^2 P)/(\partial \mu_n \partial \mu_p) \f$
60  /// The quantity \f$ (\partial^2 P)/(\partial \mu_p^2) \f$
61  fp_t dnpdmup;
62 
63  };
64 
66 
67  /** \brief Object to store second derivatives of
68  \f$ f(n_n,n_p,T) \f$
69  */
70  template<class fp_t=double>
72 
73  public:
74 
75  /// The quantity \f$ (\partial^2 P)/(\partial T^2) \f$
76  fp_t dsdT;
77  /// The quantity \f$ (\partial^2 P)/(\partial T \partial n_n) \f$
78  fp_t dmundT;
79  /// The quantity \f$ (\partial^2 P)/(\partial T \partial n_p) \f$
80  fp_t dmupdT;
81  /// The quantity \f$ (\partial^2 P)/(\partial n_n^2) \f$
82  fp_t dmundnn;
83  /// The quantity \f$ (\partial^2 P)/(\partial n_n \partial n_p) \f$
85  /// The quantity \f$ (\partial^2 P)/(\partial n_p^2) \f$
86  fp_t dmupdnp;
87  };
88 
90 
91  /** \brief Particle derivatives in the pressure representation
92 
93  This class adds the derivatives \ref dndmu, \ref dndT, and
94  \ref dsdT, which correspond to
95  \f[
96  \left(\frac{d n}{d \mu}\right)_{T,V}, \quad
97  \left(\frac{d n}{d T}\right)_{\mu,V}, \quad \mathrm{and} \quad
98  \left(\frac{d s}{d T}\right)_{\mu,V}
99  \f]
100  respectively. All other first-order thermodynamic derivatives
101  can be expressed in terms of the first three derivatives.
102 
103  \comment
104 
105  (This is no longer required)
106 
107  In the
108  case that the particle is interacting (i.e. \ref
109  part::non_interacting is \c false), then the derivatives which
110  are computed are
111  \f[
112  \left(\frac{d n}{d \nu}\right)_{T,V}, \quad
113  \left(\frac{d n}{d T}\right)_{\nu,V}, \quad
114  \left(\frac{d s}{d T}\right)_{\nu,V}, \quad \mathrm{and} \quad
115  \left(\frac{d n}{d m^{*}}\right)_{T,\nu,V},
116  \f]
117  If the particles are interacting, no derivative with respect to
118  the bare mass is given, since classes cannot know how to relate
119  the effective mass to the bare mass.
120 
121  \endcomment
122  */
123  template<class fp_t=double>
125 
126  public:
127 
128  /// Derivative of number density with respect to chemical potential
129  fp_t dndmu;
130 
131  /// Derivative of number density with respect to temperature
132  fp_t dndT;
133 
134  /// Derivative of entropy density with respect to temperature
135  fp_t dsdT;
136 
138  }
139 
140  /// Copy constructor
142  dndmu=p.dndmu;
143  dndT=p.dndT;
144  dsdT=p.dsdT;
145  }
146 
147  /// Copy construction with operator=()
149  if (this!=&p) {
150  dndmu=p.dndmu;
151  dndT=p.dndT;
152  dsdT=p.dsdT;
153  }
154  return *this;
155  }
156 
157  /** \brief Compute derivatives in the Helmholtz free energy
158  representation from the derivatives in the pressure
159  representation
160  */
161  void deriv_f(fp_t &dmudn, fp_t &dmudT, fp_t &dsdT_n) {
162  dmudn=1.0/dndmu;
163  dmudT=-dndT/dndmu;
164  dsdT_n=dsdT-dndT*dndT/dndmu;
165  return;
166  }
167  };
168 
169  /** \brief Double-precision version of \ref o2scl::part_deriv_press_tl
170  */
172 
173  /** \brief A fermion with derivative information
174  */
175  template<class fp_t=double>
176  class fermion_deriv_tl : public fermion_tl<fp_t>,
177  public part_deriv_press_tl<fp_t> {
178 
179  public:
180 
181  /// Make a particle of mass \c mass and degeneracy \c dof.
182  fermion_deriv_tl(fp_t mass=0.0, fp_t dof=0.0) : fermion(mass,dof) {
183  }
184 
185  /// Copy constructor
187  this->g=p.g;
188  this->m=p.m;
189  this->ms=p.ms;
190  this->n=p.n;
191  this->ed=p.ed;
192  this->pr=p.pr;
193  this->mu=p.mu;
194  this->en=p.en;
195  this->nu=p.nu;
196  this->dndmu=p.dndmu;
197  this->dndT=p.dndT;
198  this->dsdT=p.dsdT;
201  }
202 
203  /// Copy constructor
205  this->g=p.g;
206  this->m=p.m;
207  this->ms=p.ms;
208  this->n=p.n;
209  this->ed=p.ed;
210  this->pr=p.pr;
211  this->mu=p.mu;
212  this->en=p.en;
213  this->nu=p.nu;
214  this->dndmu=0.0;
215  this->dndT=0.0;
216  this->dsdT=0.0;
219  }
220 
221  /// Copy construction with operator=()
223  if (this!=&p) {
224  this->g=p.g;
225  this->m=p.m;
226  this->ms=p.ms;
227  this->n=p.n;
228  this->ed=p.ed;
229  this->pr=p.pr;
230  this->mu=p.mu;
231  this->en=p.en;
232  this->nu=p.nu;
233  this->dndmu=p.dndmu;
234  this->dndT=p.dndT;
235  this->dsdT=p.dsdT;
238  }
239  return *this;
240  }
241 
242  /// Copy construction with operator=()
244  if (this!=&p) {
245  this->g=p.g;
246  this->m=p.m;
247  this->ms=p.ms;
248  this->n=p.n;
249  this->ed=p.ed;
250  this->pr=p.pr;
251  this->mu=p.mu;
252  this->en=p.en;
253  this->nu=p.nu;
254  this->dndmu=0.0;
255  this->dndT=0.0;
256  this->dsdT=0.0;
259  }
260  return *this;
261  }
262 
263  };
264 
265  typedef fermion_deriv_tl<double> fermion_deriv;
266 
267  /** \brief A part with derivative information
268  */
269  template<class fp_t=double>
270  class part_deriv_tl : public part_tl<fp_t>,
271  public part_deriv_press_tl<fp_t> {
272 
273  public:
274 
275  /// Make a particle of mass \c mass and degeneracy \c dof.
276  part_deriv_tl(fp_t mass=0.0, fp_t dof=0.0) : part(mass,dof) {
277  }
278 
279  /// Copy constructor
281  this->g=p.g;
282  this->m=p.m;
283  this->ms=p.ms;
284  this->n=p.n;
285  this->ed=p.ed;
286  this->pr=p.pr;
287  this->mu=p.mu;
288  this->en=p.en;
289  this->nu=p.nu;
290  this->dndmu=p.dndmu;
291  this->dndT=p.dndT;
292  this->dsdT=p.dsdT;
295  }
296 
297  /// Copy construction with operator=()
299  if (this!=&p) {
300  this->g=p.g;
301  this->m=p.m;
302  this->ms=p.ms;
303  this->n=p.n;
304  this->ed=p.ed;
305  this->pr=p.pr;
306  this->mu=p.mu;
307  this->en=p.en;
308  this->nu=p.nu;
309  this->dndmu=p.dndmu;
310  this->dndT=p.dndT;
311  this->dsdT=p.dsdT;
314  }
315  return *this;
316  }
317 
318  };
319 
320  typedef part_deriv_tl<double> part_deriv;
321 
322  /** \brief Base quantities for thermodynamic derivatives
323 
324  The quantities \f$ c_P \f$ computed by
325  \ref heat_cap_ppart_const_press(), \f$ c_V \f$
326  computed by \ref heat_cap_ppart_const_vol(),
327  \f$ \beta_T \f$ computed by \ref compress_const_tptr(),
328  \f$ \beta_S \f$ computed by \ref compress_adiabatic(),
329  and \f$ \alpha_V \f$ computed by \ref coeff_thermal_exp
330  are related by
331  \f[
332  c_P - c_V = \frac{T \alpha_V^2}{n \beta_T}
333  \f]
334  and
335  \f[
336  \beta_T - \beta_S = \frac{T \alpha_V^2}{n c_P}
337  \f]
338 
339  For the derivatives below, the following
340  Jacobian is useful
341  \f{eqnarray*}
342  \frac{\partial (P,S,N)}{\partial (V,\mu,T)}
343  &=& -n \left[
344  \left(\frac{\partial S}{\partial V}\right)_{\mu,T}
345  \left(\frac{\partial N}{\partial T}\right)_{\mu,V}
346  - \left(\frac{\partial S}{\partial T}\right)_{\mu,V}
347  \left(\frac{\partial N}{\partial V}\right)_{\mu,T}
348  \right] +
349  s \left[
350  \left(\frac{\partial S}{\partial V}\right)_{\mu,T}
351  \left(\frac{\partial N}{\partial \mu}\right)_{V,T}
352  - \left(\frac{\partial S}{\partial \mu}\right)_{V,T}
353  \left(\frac{\partial N}{\partial V}\right)_{\mu,T}
354  \right]
355  \nonumber \\
356  &=& - V n \left[
357  s
358  \left(\frac{\partial n}{\partial T}\right)_{\mu}
359  - n \left(\frac{\partial s}{\partial T}\right)_{\mu}
360  \right] + V s \left[ s
361  \left(\frac{\partial n}{\partial \mu}\right)_{T}
362  - n \left(\frac{\partial n}{\partial T}\right)_{T}
363  \right]
364  = V n^2 \left(\frac{\partial s}{\partial T}\right)_{\mu}
365  - 2 V n s \left(\frac{\partial n}{\partial T}\right)_{\mu}
366  + V s^2 \left(\frac{\partial n}{\partial \mu}\right)_{T}
367  \f}
368  For convenience, we define the quantity
369  \f[
370  X \equiv \frac{1}{V}
371  \left[ \frac{\partial (P,S,N)}{\partial (V,\mu,T)} \right]
372  \f]
373  Another common combination of derivatives is
374  \f[
375  Y \equiv \left(\frac{\partial n}{\partial T}\right)_{\mu}^2 -
376  \left(\frac{\partial s}{\partial T}\right)_{\mu}
377  \left(\frac{\partial n}{\partial \mu}\right)_{T}
378  \f]
379 
380  */
381  template<class fp_t=double> class deriv_thermo_base_tl {
382 
383  public:
384 
385  /** \brief The heat capacity per particle at
386  constant volume (unitless)
387 
388  This function returns
389  \f[
390  c_V = \frac{T}{N}
391  \left(\frac{\partial S}{\partial T}\right)_{V,N} =
392  \frac{T}{n} \left(\frac{\partial s}{\partial T}\right)_{V,n} =
393  \frac{1}{N} \left(\frac{\partial E}{\partial T}\right)_{V,N}
394  \f]
395 
396  To write this in terms of the three derivatives in
397  \ref o2scl::part_deriv_press_tl,
398  \f[
399  \frac{T}{n} \left(\frac{\partial s}{\partial T}\right)_{V,n}
400  = \frac{T}{n} \frac{\partial(s,n,V)}{\partial(T,n,V)} =
401  \frac{T}{n} \left[\frac{\partial(s,n,V)}{\partial(T,\mu,V)}\right]
402  \left[\frac{\partial(T,n,V)}{\partial(T,\mu,V)}\right]^{-1}
403  \f]
404  \f[
405  = \frac{T}{n}
406  \left[
407  \left(\frac{\partial s}{\partial T}\right)_{\mu} -
408  \left(\frac{\partial n}{\partial T}\right)_{\mu}^2
409  \left(\frac{\partial n}{\partial \mu}\right)_{T}^{-1}
410  \right]
411  \f]
412 
413  This is \f$ 3/2 \f$ for an ideal gas.
414  */
415  template<class part_deriv_t>
416  fp_t heat_cap_ppart_const_vol(part_deriv_t &p, fp_t temper) {
417  return (p.dsdT-p.dndT*p.dndT/p.dndmu)*temper/p.n;
418  }
419 
420  /** \brief The heat capacity per particle
421  at constant pressure (unitless)
422 
423  This function returns
424  \f[
425  c_P = \frac{T}{N}
426  \left(\frac{\partial S}{\partial T}\right)_{P,N} =
427  \frac{1}{N} \left(\frac{\partial H}{\partial T}\right)_{P,N}
428  \f]
429 
430  To write this in terms of the three derivatives in
431  \ref o2scl::part_deriv_press_tl,
432  \f[
433  \frac{T}{N} \left(\frac{\partial S}{\partial T}\right)_{P,N}
434  = \frac{T}{N} \frac{\partial(S,N,P)}{\partial(T,N,P)} =
435  \frac{T}{N} \left[\frac{\partial(S,N,P)}{\partial(T,\mu,V)}\right]
436  \left[\frac{\partial(T,N,P)}{\partial(T,\mu,V)}\right]^{-1}
437  \f]
438  The first Jacobian was computed above since
439  \f[
440  \frac{\partial(S,N,P)}{\partial(T,\mu,V)} = -
441  \frac{\partial(P,S,N)}{\partial(V,\mu,T)}
442  \f]
443  The second is
444  \f[
445  \frac{\partial(T,N,P)}{\partial(T,\mu,V)}
446  =
447  \left[
448  \left(\frac{\partial N}{\partial \mu}\right)_{T,V}
449  \left(\frac{\partial P}{\partial V}\right)_{\mu,T}
450  - \left(\frac{\partial N}{\partial V}\right)_{\mu,T}
451  \left(\frac{\partial P}{\partial \mu}\right)_{T,V}
452  \right] = - n^2
453  \f]
454  The final result is
455  \f[
456  c_P = \frac{T X}{n^3} =
457  \frac{T}{n} \left(\frac{\partial s}{\partial T}\right)_{\mu}
458  + \frac{s^2 T}{n^3} \left(\frac{\partial n}{\partial \mu}\right)_{T}
459  - \frac{2 s T}{n^2} \left(\frac{\partial n}{\partial T}\right)_{\mu}
460  \f]
461 
462  This is \f$ 5/2 \f$ for an ideal gas.
463  */
464  template<class part_deriv_t>
465  fp_t heat_cap_ppart_const_press(part_deriv_t &p, fp_t temper) {
466  return temper/p.n*p.dsdT+p.en*p.en*temper/p.n/p.n/p.n*p.dndmu-
467  2.0*p.en*temper/p.n/p.n*p.dndT;
468  }
469 
470  /** \brief The adiabatic compressibility
471 
472  This function computes
473  \f[
474  \beta_S \equiv - \frac{1}{V}
475  \left(\frac{\partial V}{\partial P}\right)_{S,N}
476  \f]
477  (sometimes referred to as \f$ \kappa_S \f$ or
478  \f$ \chi_S \f$)
479 
480  To write this in terms of the three derivatives in
481  \ref o2scl::part_deriv_press_tl,
482  \f[
483  \left(\frac{\partial V}{\partial P}\right)_{S,N} =
484  \frac{\partial (V,S,N)}{\partial (P,S,N)} =
485  \frac{\partial (V,S,N)}{\partial (V,\mu,T)}
486  \left[ \frac{\partial (P,S,N)}{\partial (V,\mu,T)}\right]^{-1}
487  \f]
488  The first Jacobian
489  \f[
490  \frac{\partial (V,S,N)}{\partial (V,\mu,T)} = V^2
491  \left[
492  \left(\frac{\partial s}{\partial T}\right)_{\mu,V}
493  \left(\frac{\partial n}{\partial \mu}\right)_{T,V}
494  - \left(\frac{\partial n}{\partial T}\right)_{\mu,V}^2
495  \right]
496  \f]
497  and the second Jacobian was computed above.
498  The result is
499  \f[
500  \beta_S = Y/X = \left[
501  \left(\frac{\partial n}{\partial T}\right)_{\mu}^2 -
502  \left(\frac{\partial s}{\partial T}\right)_{\mu}
503  \left(\frac{\partial n}{\partial \mu}\right)_{T}
504  \right]
505  \left[
506  n^2 \left(\frac{\partial s}{\partial T}\right)_{\mu,V}
507  - 2 n s \left(\frac{\partial n}{\partial T}\right)_{\mu,V}
508  + s^2 \left(\frac{\partial n}{\partial \mu}\right)_{T,V}
509  \right]^{-1}
510  \f]
511  */
512  template<class part_deriv_t>
513  fp_t compress_adiabatic(part_deriv_t &p, fp_t temper) {
514  return (p.dndT*p.dndT-p.dndmu*p.dsdT)/
515  (p.n*p.n*p.dsdT-2.0*p.n*p.en*p.dndT+p.en*p.en*p.dndmu);
516  }
517 
518  /** \brief The isothermal compressibility
519 
520  This function computes
521  \f[
522  \beta_T \equiv - \frac{1}{V}
523  \left(\frac{\partial V}{\partial P}\right)_{T,N}
524  \f]
525  (sometimes referred to as \f$ \kappa_T \f$ or
526  \f$ \chi_T \f$) in units of inverse length to the fourth
527  power.
528 
529  To write this in terms of the three derivatives in
530  \ref o2scl::part_deriv_press_tl,
531  \f{eqnarray*}
532  - \frac{1}{V} \left(\frac{\partial V}{\partial P}\right)_{T,N} &=&
533  \frac{\partial (V,T,N)}{\partial (P,T,N)} =
534  \frac{1}{V}
535  \frac{\partial (V,T,N)}{\partial (V,T,\mu)}
536  \left[\frac{\partial (N,P,T)}{\partial (V,\mu,T)}\right]^{-1}
537  \nonumber \\
538  &=& \left(\frac{\partial n}{\partial \mu}\right)_{T,V}
539  \left[
540  \left(\frac{\partial N}{\partial V}\right)_{\mu,T}
541  \left(\frac{\partial P}{\partial \mu}\right)_{V,T}
542  - \left(\frac{\partial P}{\partial V}\right)_{\mu,T}
543  \left(\frac{\partial N}{\partial \mu}\right)_{V,T}
544  \right]^{-1} =
545  \frac{1}{n^2} \left(\frac{\partial n}{\partial \mu}\right)_{T}
546  \f}
547  */
548  template<class part_deriv_t>
549  fp_t compress_const_tptr(part_deriv_t &p, fp_t temper) {
550  return p.dndmu/p.n/p.n;
551  }
552 
553  /** \brief The coefficient of thermal expansion
554 
555  This function computes
556  \f[
557  \alpha_V =
558  \frac{1}{V} \left(\frac{\partial V}{\partial T}\right)_{P,N}
559  \f]
560  in units of length.
561 
562  To write this in terms of the three derivatives in
563  \ref o2scl::part_deriv_press_tl,
564  \f{eqnarray*}
565  \left(\frac{\partial V}{\partial T}\right)_{P,N} &=&
566  \frac{\partial (V,P,N)}{\partial (T,P,N)} =
567  -\frac{\partial (V,P,N)}{\partial (V,T,\mu)}
568  \left[ \frac{\partial (T,P,N)}{\partial (T,V,\mu)} \right]^{-1}
569  \nonumber \\
570  & = &
571  - \left[
572  \left(\frac{\partial P}{\partial T}\right)_{\mu,V}
573  \left(\frac{\partial N}{\partial \mu}\right)_{T,V} -
574  \left(\frac{\partial N}{\partial T}\right)_{\mu,V}
575  \left(\frac{\partial P}{\partial \mu}\right)_{T,V}
576  \right]
577  \left[
578  \left(\frac{\partial P}{\partial V}\right)_{\mu,T}
579  \left(\frac{\partial N}{\partial \mu}\right)_{V,T} -
580  \left(\frac{\partial P}{\partial \mu}\right)_{V,T}
581  \left(\frac{\partial N}{\partial V}\right)_{\mu,T}
582  \right]^{-1}
583  \nonumber \\
584  &=& \frac{s}{n^2}
585  \left(\frac{\partial n}{\partial \mu}\right)_{T} -
586  \frac{1}{n} \left(\frac{\partial n}{\partial T}\right)_{\mu}
587  \f}
588  */
589  template<class part_deriv_t>
590  fp_t coeff_thermal_exp(part_deriv_t &p, fp_t temper) {
591  return p.en/p.n/p.n*p.dndmu-p.dndT/p.n;
592  }
593 
594  /** \brief The squared sound speed (unitless)
595 
596  This function computes the squared sound speed
597  (including relativistic effects)
598  \f[
599  c_s^2 = \left(\frac{\partial P}
600  {\partial \varepsilon}\right)_{S,N}
601  \f]
602  The result is unitless. To get the units of a squared velocity,
603  one must multiply by \f$ c^2 \f$ .
604 
605  The
606  nonrelativistic squared sound speed
607  is
608  \f[
609  c_{s,\mathrm{NR}}^2 = \left[\frac{\partial P}
610  {\partial (N/V)}\right]_{S,N} =
611  - \frac{V^2}{N} \left(\frac{\partial P}
612  {\partial V}\right)_{S,N} = \frac{1}{n \beta_S}
613  \f]
614  where \f$ \beta_S \f$
615  is computed in \ref compress_adiabatic() .
616 
617  To write \f$ c_s^2 \f$ in terms of the three derivatives in
618  \ref o2scl::part_deriv_press_tl,
619  \f[
620  \left(\frac{\partial P}
621  {\partial \varepsilon}\right)_{S,N} =
622  \frac{\partial (P,S,N)}{\partial (\varepsilon,S,N)} =
623  \frac{\partial (P,S,N)}{\partial (V,T,\mu)}
624  \left[ \frac{\partial (\varepsilon,S,N)}
625  {\partial (V,T,\mu)} \right]^{-1}
626  \f]
627  The first Jacobian was computed above (up to a sign).
628  The second is the determinant of
629  \f[
630  \left(
631  \begin{array}{ccc}
632  0
633  & \frac{\partial \varepsilon}{\partial T}
634  & \frac{\partial \varepsilon}{\partial \mu} \\
635  s & V \frac{\partial s}{\partial T}
636  & V \frac{\partial n}{\partial T} \\
637  n & V \frac{\partial n}{\partial T}
638  & V \frac{\partial n}{\partial \mu}
639  \end{array}
640  \right)
641  \f]
642  with
643  \f[
644  \frac{\partial \varepsilon}{\partial T} =
645  -T \frac{\partial s}{\partial T}
646  + \mu \frac{\partial n}{\partial T}
647  \quad \mathrm{and} \quad
648  \frac{\partial \varepsilon}{\partial \mu} =
649  T \frac{\partial n}{\partial T}
650  + \mu \frac{\partial n}{\partial \mu}
651  \f]
652  giving
653  \f[
654  \frac{\partial (\varepsilon,S,N)}
655  {\partial (V,T,\mu)} = V
656  (P + \varepsilon)
657  \left[ \left(\frac{\partial n}{\partial T}\right)^2
658  - \left(\frac{\partial n}{\partial \mu}\right)
659  \left(\frac{\partial s}{\partial T}\right)
660  \right] = V Y \left(P+\varepsilon\right)
661  \f]
662  The final result is
663  \f[
664  c_s^2 =
665  - \frac{X}{(P+\varepsilon)Y}
666  =
667  \frac{
668  n^2 \left(\frac{\partial s}{\partial T}\right)
669  - 2 n s \left(\frac{\partial n}{\partial T}\right)
670  + s^2 \left(\frac{\partial n}{\partial \mu}\right)
671  }{
672  \left(P + \varepsilon\right)
673  \left[
674  \left(\frac{\partial n}{\partial \mu}\right)
675  \left(\frac{\partial s}{\partial T}\right) -
676  \left(\frac{\partial n}{\partial T}\right)^2
677  \right]
678  }
679  \f]
680 
681  */
682  template<class part_deriv_t>
683  fp_t squared_sound_speed(part_deriv_t &p, fp_t temper) {
684  fp_t edt;
685  if (p.inc_rest_mass) {
686  edt=p.ed;
687  } else {
688  edt=p.ed+p.n*p.m;
689  }
690  return (p.n*p.n*p.dsdT-2.0*p.n*p.en*p.dndT+p.en*p.en*p.dndmu)/
691  (edt+p.pr)/(p.dndmu*p.dsdT-p.dndT*p.dndT);
692  }
693 
694  };
695 
696  typedef deriv_thermo_base_tl<double> deriv_thermo_base;
697 
698  /** \brief Compute properties of a fermion including derivatives
699  [abstract base]
700 
701  \future Include explicit zero-temperature calculation, maybe
702  by making this a child of fermion_zerot or by making a
703  new fermion_deriv_zerot?
704  \comment
705  dn/dmu is just g*mu*kf/2/pi^2
706  \endcomment
707  \future There is also a closed form for the derivatives
708  of massless fermions with pairs at finite temperature
709  in Constantinou et al. 2014 which could be implemented here.
710  */
711  template<class fp_t=double>
713 
714  protected:
715 
716  /** \brief A fermion_thermo object
717 
718  This is for access to fermion_thermo::ndeg_terms().
719  */
721 
722  /// Desc
723  fp_t pi;
724 
725  /// Desc
726  fp_t pi2;
727 
728  public:
729 
731  pi=boost::math::constants::pi<fp_t>();
732  pi2=boost::math::constants::pi_sqr<fp_t>();
733  }
734 
735  virtual ~fermion_deriv_thermo_tl() {
736  }
737 
738  /** \brief Calculate properties as function of chemical potential
739  */
740  virtual int calc_mu(fermion_deriv &f, fp_t temper)=0;
741 
742  /** \brief Calculate properties as function of density
743  */
744  virtual int calc_density(fermion_deriv &f, fp_t temper)=0;
745 
746  /** \brief Calculate properties with antiparticles as function of
747  chemical potential
748  */
749  virtual int pair_mu(fermion_deriv &f, fp_t temper)=0;
750 
751  /** \brief Calculate properties with antiparticles as function of
752  density
753  */
754  virtual int pair_density(fermion_deriv &f, fp_t temper)=0;
755 
756  /// Calculate effective chemical potential from density
757  virtual int nu_from_n(fermion_deriv &f, fp_t temper)=0;
758 
759  /** \brief Calculate properties as a function of chemical
760  potential using a degenerate expansion
761 
762  \future There is some repetition of the code
763  for this function and the function
764  \ref o2scl::fermion_thermo_tl::calc_mu_deg() .
765  which could be avoided.
766  */
767  virtual bool calc_mu_deg(fermion_deriv &f, fp_t temper,
768  fp_t prec) {
769 
770  if (fr.calc_mu_deg_tlate<fermion_deriv>(f,temper,prec)==false) {
771  return false;
772  }
773 
774  // Compute psi and tt
775  fp_t psi;
776  if (f.inc_rest_mass) psi=(f.nu-f.ms)/temper;
777  else psi=(f.nu+f.m-f.ms)/temper;
778  fp_t tt=temper/f.ms;
779 
780  // Prefactor 'd' in Johns96
781  fp_t prefac=f.g/2.0/o2scl_const::pi2*pow(f.ms,4.0);
782 
783  // Define x = psi * t = (mu/m - 1) and related values
784  fp_t x=psi*tt;
785  fp_t sx=sqrt(x);
786  fp_t s2x=sqrt(2.0+x);
787  fp_t x2=x*x;
788  fp_t x3=x2*x;
789  fp_t x4=x2*x2;
790 
791  // First order density term (first order entropy term is zero)
792  fp_t dndmu_term1=sx*s2x*(1.0+x)/f.ms/f.ms;
793 
794  // Second order terms
795  fp_t dndmu_term2=tt*tt*o2scl_const::pi2/6.0*(1.0+x)*
796  (-1.0+2.0*x*(2.0+x))/
797  f.ms/f.ms/sx/s2x/x/(2.0+x);
798  fp_t dndT_term2=tt*o2scl_const::pi2/3.0*(1.0+2.0*x*(2.0+x))/
799  f.ms/f.ms/sx/s2x;
800  fp_t dsdT_term2=o2scl_const::pi2/3.0*(1.0+x)*sx*s2x/
801  f.ms/f.ms;
802 
803  // Third order terms
804  fp_t dndmu_term3=-7.0*pow(o2scl_const::pi*tt,4.0)/24.0*
805  (1.0+x)/sx/s2x/x3/(x+2.0)/(x+2.0)/(x+2.0)/f.ms/f.ms;
806  fp_t dndT_term3=7.0*pow(o2scl_const::pi*tt,4.0)/tt/30.0/
807  f.ms/f.ms/pow(x*(2.0+x),2.5);
808  fp_t dsdT_term3=7.0*pow(o2scl_const::pi*tt,2.0)*
809  o2scl_const::pi2/30.0/f.ms/f.ms*(1.0+x)*(-1.0+2.0*x*(2.0+x))/
810  x/(2.0+x)/sx/s2x;
811 
812  // Fourth order terms for density and entropy
813  fp_t dndmu_term4=-31.0*pow(o2scl_const::pi*tt,6.0)/48.0*(1.0+x)*
814  (3.0+2.0*x*(2.0+x))/f.ms/f.ms/pow(x*(2.0+x),5.5);
815  fp_t dndT_term4=31.0*pow(o2scl_const::pi*tt,6.0)/tt/168.0*
816  (7.0+6.0*x*(2.0+x))/f.ms/f.ms/pow(x*(2.0+x),4.5);
817  fp_t dsdT_term4=-155.0*pow(o2scl_const::pi*tt,4.0)*
818  o2scl_const::pi2/168.0*
819  (1.0+x)/f.ms/f.ms/pow(x*(2.0+x),3.5);
820 
821  // Add up all the terms
822  f.dndmu=prefac*(dndmu_term1+dndmu_term2+dndmu_term3+dndmu_term4);
823  f.dndT=prefac*(dndT_term2+dndT_term3+dndT_term4);
824  f.dsdT=prefac*(dsdT_term2+dsdT_term3+dsdT_term4);
825 
826  return true;
827  }
828 
829  /** \brief Calculate properties as a function of chemical
830  potential using a nondegenerate expansion
831 
832  \future There is some repetition of the code
833  for this function and the function
834  \ref o2scl::fermion_thermo_tl::calc_mu_ndeg() .
835  which could be avoided.
836  */
837  virtual bool calc_mu_ndeg(fermion_deriv &f, fp_t temper,
838  fp_t prec, bool inc_antip=false) {
839 
840  if (fr.calc_mu_ndeg_tlate<fermion_deriv>(f,temper,prec,
841  inc_antip)==false) {
842  return false;
843  }
844 
845  // Compute psi and tt
846  fp_t psi, psi_num;
847  if (f.inc_rest_mass) {
848  psi_num=f.nu-f.ms;
849  } else {
850  psi_num=f.nu+f.m-f.ms;
851  }
852  psi=psi_num/temper;
853  fp_t tt=temper/f.ms;
854  fp_t xx=psi*tt;
855 
856  // Prefactor 'd' in Johns96
857  fp_t prefac=f.g/2.0/o2scl_const::pi2*pow(f.ms,4.0);
858 
859  // One term is always used, so only values of max_term greater than
860  // 0 are useful.
861  static const size_t max_term=200;
862 
863  fp_t first_dndT=0.0;
864  fp_t first_dsdT=0.0;
865  fp_t first_dndmu=0.0;
866 
867  fp_t nu2=f.nu;
868  if (f.inc_rest_mass==false) nu2+=f.m;
869 
870  f.dndmu=0.0;
871  f.dndT=0.0;
872  f.dsdT=0.0;
873 
874  for(size_t j=1;j<=max_term;j++) {
875 
876  fp_t dj=((fp_t)j);
877  fp_t jot=dj/tt;
878 
879  // Here, we are only computing the derivatives, but we need to
880  // compute the terms in the pressure, density, and entropy because
881  // they are used in computing the terms for the derivatives.
882  fp_t pterm, nterm, enterm;
883  fp_t dndmu_term, dndT_term, dsdT_term;
884 
885  fr.ndeg_terms(j,tt,psi*tt,f.ms,f.inc_rest_mass,inc_antip,
886  pterm,nterm,enterm);
887 
888  if (inc_antip==false) {
889  dndmu_term=nterm*jot;
890  dndT_term=jot*enterm-nterm/tt;
891  dsdT_term=(3.0*tt-2.0*dj*xx-2.0*dj)/tt/tt*enterm+
892  (5.0*dj*tt-2.0*dj*dj*xx+5.0*dj*tt*xx-dj*dj*xx*xx)/
893  dj/tt/tt/tt*nterm;
894  } else {
895  dndmu_term=nterm*jot;
896  dndT_term=jot*enterm*tanh(jot*(xx+1.0))-
897  (tt+2.0*dj*(1.0+xx))/sinh(jot*(xx+1.0))*nterm/tt/tt;
898  dsdT_term=(2.0*dj*(1.0+xx)*tanh(jot*(xx+1.0))-3.0*tt)*enterm/tt/tt+
899  (2.0*pow(dj*1.0+xx,2.0)*tanh(jot*(xx+1.0))-
900  dj*dj*(2.0+2.0*xx+xx*xx)*cosh(jot*(xx+1.0))-
901  5.0*dj*(1.0+xx)*tt)*nterm/dj/tt/tt/tt;
902  }
903  dndmu_term/=f.ms;
904  dndT_term/=f.ms;
905  dsdT_term/=f.ms;
906 
907  if (j==1) {
908  first_dndT=dndT_term;
909  first_dsdT=dsdT_term;
910  first_dndmu=dndmu_term;
911  }
912  f.dndmu+=dndmu_term;
913  f.dndT+=dndT_term;
914  f.dsdT+=dsdT_term;
915 
916  /*
917  cout << j << " " << dj << " " << tt << " " << xx << " "
918  << f.ms << " " << pterm << " " << nterm << " "
919  << enterm << " "
920  << dndmu_term << " " << dndT_term << endl;
921  */
922 
923  // If the first terms are zero, then the rest of the terms
924  // will be zero so just return early
925  if (first_dndT==0.0 && first_dndmu==0.0 && first_dsdT==0.0) {
926  f.dndmu=0.0;
927  f.dndT=0.0;
928  f.dsdT=0.0;
929  return true;
930  }
931 
932  // Stop if the last term is sufficiently small compared to
933  // the first term
934  if (j>1 &&
935  fabs(dndT_term)<prec*fabs(first_dndT) &&
936  fabs(dndmu_term)<prec*fabs(first_dndmu) &&
937  fabs(dsdT_term)<prec*fabs(first_dsdT)) {
938  f.dndT*=prefac;
939  f.dndmu*=prefac;
940  f.dsdT*=prefac;
941  return true;
942  }
943 
944  // End of 'for(size_t j=1;j<=max_term;j++)'
945  }
946 
947  // We failed to add enough terms, so return false
948  return false;
949  }
950 
951  };
952 
953  typedef fermion_deriv_thermo_tl<double> fermion_deriv_thermo;
954 
955  /** \brief Object to organize calibration of derivative quantities
956  in particle classes to results stored in a table
957  */
959 
960  public:
961 
962  /** \brief Calibrate a particle thermodynamics class with
963  derivatives with results stored in a table
964 
965  This compares the approximation to the exact results using
966  calc_density(), calc_mu(), pair_density() and pair_mu(). It
967  tries each function twelve times. It tries three different
968  temperatures, setting both <tt>inc_rest_mass</tt> and
969  <tt>non_interacting</tt> equal to <tt>true</tt> and
970  <tt>false</tt>.
971 
972  The <tt>verbose</tt> parameter controls the amount of output.
973  */
974  template<class part_t, class thermo_t>
975  double part_deriv_calibrate(part_t &p, thermo_t &th, bool test_pair,
976  std::string file, bool nr_mode=false,
977  int verbose=0, bool external=false) {
978 
979  double ret=0;
980 
981  // ----------------------------------------------------------------
982  // Will return to these original values afterwards
983 
984  part_t orig=p;
985 
986  // ----------------------------------------------------------------
987  // Read data file
988 
989  std::string fname;
990  if (external==false) {
991  fname=o2scl_settings.get_data_dir()+file;
992  } else {
993  fname=file;
994  }
995 
996  if (verbose>1) {
997  std::cout << "In part_calibrate(), loading file named\n\t'"
998  << fname << "'.\n" << std::endl;
999  }
1000  o2scl::table<> tab;
1002  hf.open(fname);
1003  std::string name;
1004 #ifndef O2SCL_NO_HDF_INPUT
1005  hdf_input(hf,tab,name);
1006 #endif
1007  hf.close();
1008 
1009  if (tab.get_nlines()==0) {
1010  std::string str="Failed to load data from file '"+fname+
1011  "' in part_calibrate(). Bad filename?";
1012  O2SCL_ERR(str.c_str(),exc_efilenotfound);
1013  }
1014 
1015  if (!tab.is_column("ed")) {
1016  tab.function_column("ed_mot*mot","ed");
1017  if (test_pair) {
1018  tab.function_column("pair_ed_mot*mot","pair_ed");
1019  }
1020  }
1021 
1022  // ----------------------------------------------------------------
1023 
1024  p.g=2.0;
1025 
1026  size_t cnt=0;
1027  part_t bad, dev, exact;
1028  double m_bad=0.0, mu_bad=0.0, T_bad=0.0, mot_bad=0.0, psi_bad=0.0;
1029  p.non_interacting=true;
1030 
1031  // ----------------------------------------------------------------
1032  // First pass, test calc_mu()
1033 
1034  // k=0,2 are with rest mass, k=1,3 are without
1035  // k=0,1 are non-interacting, k=2,3 are interacting
1036  for(size_t k=0;k<4;k++) {
1037 
1038  double ret_local=0.0;
1039 
1040  // Initialize storage
1041  dev.n=0.0; dev.dndmu=0.0; dev.dndT=0.0; dev.dsdT=0.0;
1042  bad.n=0.0; bad.dndmu=0.0; bad.dndT=0.0; bad.dsdT=0.0;
1043 
1044  // Temperature loop
1045  for(double T=1.0e-2;T<=1.001e2;T*=1.0e2) {
1046 
1047  // Loop over each point in the data file
1048  for(size_t i=0;i<tab.get_nlines();i++) {
1049 
1050  double mot=tab.get("mot",i);
1051  double psi=tab.get("psi",i);
1052  exact.n=tab.get("n",i);
1053  exact.dndmu=tab.get("dndmu",i);
1054  exact.dndT=tab.get("dndT",i);
1055  exact.dsdT=tab.get("dsdT",i);
1056 
1057  set_mass_flags(p,mot,T,k);
1058  set_chem_pot(p,psi,T,k,nr_mode);
1059 
1060  th.calc_mu(p,T);
1061 
1062  exact.n*=pow(T,3.0);
1063  exact.dndmu*=pow(T,2.0);
1064  exact.dndT*=pow(T,2.0);
1065  exact.dsdT*=pow(T,2.0);
1066 
1067  dev.dndmu+=fabs((p.dndmu-exact.dndmu)/exact.dndmu);
1068  dev.dndT+=fabs((p.dndT-exact.dndT)/exact.dndT);
1069  dev.dsdT+=fabs((p.dsdT-exact.dsdT)/exact.dsdT);
1070 
1071  cnt++;
1072 
1073  check_derivs<part_t>(p,exact,bad,k,T,mot,psi,mu_bad,m_bad,T_bad,
1074  mot_bad,psi_bad,ret_local);
1075 
1076  if (verbose>1) {
1077  std::cout.precision(5);
1078  if (k>=2) {
1079  std::cout << "T,ms,nu,psi,mot: " << T << " "
1080  << p.ms << " " << p.nu
1081  << " " << psi << " " << mot << std::endl;
1082  } else {
1083  std::cout << "T,m,mu,psi,mot: " << T << " "
1084  << p.m << " " << p.mu
1085  << " " << psi << " " << mot << std::endl;
1086  }
1087  std::cout.precision(5);
1088  std::cout << "n,dndmu,dndT,dsdT: " << std::endl;
1089  std::cout << "approx: " << p.n << " " << p.dndmu << " "
1090  << p.dndT << " "
1091  << p.dsdT << std::endl;
1092  std::cout << "exact : " << exact.n << " " << exact.dndmu << " "
1093  << exact.dndT << " " << exact.dsdT << std::endl;
1094  std::cout << "bad : " << bad.n << " " << bad.dndmu << " "
1095  << bad.dndT << " " << bad.dsdT << std::endl;
1096  std::cout << std::endl;
1097  if (verbose>2) {
1098  char ch;
1099  std::cin >> ch;
1100  }
1101  }
1102 
1103  if (ret_local>ret) {
1104  ret=ret_local;
1105  }
1106 
1107  // End of loop over points in data file
1108  }
1109  // End of temperature loop
1110  }
1111 
1112  dev.n/=cnt;
1113  dev.dndmu/=cnt;
1114  dev.dndT/=cnt;
1115  dev.dsdT/=cnt;
1116 
1117  if (verbose>0) {
1118  if (k==0) {
1119  std::cout << "Function calc_mu(), include rest mass:" << std::endl;
1120  } else if (k==1) {
1121  std::cout << "Function calc_mu(), without rest mass:" << std::endl;
1122  } else if (k==2) {
1123  std::cout << "Function calc_mu(), include rest mass, "
1124  << "interacting:" << std::endl;
1125  } else {
1126  std::cout << "Function calc_mu(), without rest mass, "
1127  << "interacting:" << std::endl;
1128  }
1129 
1130  std::cout << "Average performance: " << std::endl;
1131  std::cout << "dndmu: " << dev.dndmu << " dndT: "
1132  << dev.dndT << " dsdT: " << dev.dsdT << std::endl;
1133  std::cout << "Worst case: " << std::endl;
1134  std::cout << "dndmu: " << bad.dndmu << " dndT: "
1135  << bad.dndT << " dsdT: " << bad.dsdT << std::endl;
1136  std::cout << "mu: " << mu_bad << " m: " << m_bad << " T: " << T_bad
1137  << " mot: " << mot_bad
1138  << "\n\tpsi: " << psi_bad << std::endl;
1139  std::cout << std::endl;
1140  if (verbose>2) {
1141  char ch;
1142  std::cin >> ch;
1143  }
1144  }
1145 
1146  // Reset p.non_interacting
1147  p.non_interacting=true;
1148 
1149  // End of k loop
1150  }
1151 
1152  // ----------------------------------------------------------------
1153  // Second pass, test calc_density()
1154 
1155  // k=0,2 are with rest mass, k=1,3 are without
1156  // k=0,1 are non-interacting, k=2,3 are interacting
1157  for(size_t k=0;k<4;k++) {
1158 
1159  double ret_local=0.0;
1160 
1161  // Initialize storage
1162  dev.mu=0.0; dev.dndmu=0.0; dev.dndT=0.0; dev.dsdT=0.0;
1163  bad.mu=0.0; bad.dndmu=0.0; bad.dndT=0.0; bad.dsdT=0.0;
1164 
1165  // Temperature loop
1166  for(double T=1.0e-2;T<=1.001e2;T*=1.0e2) {
1167 
1168  // Loop over each point in the data file
1169  for(size_t i=0;i<tab.get_nlines();i++) {
1170 
1171  double mot=tab.get("mot",i);
1172  double psi=tab.get("psi",i);
1173  p.n=tab.get("n",i);
1174  exact.dndmu=tab.get("dndmu",i);
1175  exact.dndT=tab.get("dndT",i);
1176  exact.dsdT=tab.get("dsdT",i);
1177 
1178  set_mass_flags(p,mot,T,k);
1179  set_chem_pot(exact,psi,T,k,nr_mode);
1180 
1181  p.n*=pow(T,3.0);
1182  exact.dndmu*=pow(T,2.0);
1183  exact.dndT*=pow(T,2.0);
1184  exact.dsdT*=pow(T,2.0);
1185 
1186  // Give it a guess for the chemical potential
1187  if (k>=2) {
1188  p.nu=p.m;
1189  } else {
1190  p.mu=p.m;
1191  }
1192 
1193  th.calc_density(p,T);
1194 
1195  dev.dndmu+=fabs((p.dndmu-exact.dndmu)/exact.dndmu);
1196  dev.dndT+=fabs((p.dndT-exact.dndT)/exact.dndT);
1197  dev.dsdT+=fabs((p.dsdT-exact.dsdT)/exact.dsdT);
1198 
1199  cnt++;
1200 
1201  check_derivs<part_t>(p,exact,bad,k,T,mot,psi,mu_bad,m_bad,
1202  T_bad,mot_bad,psi_bad,ret_local);
1203 
1204  if (verbose>1) {
1205  std::cout.precision(5);
1206  if (k>=2) {
1207  std::cout << "T,ms,n,psi,mot: " << T << " "
1208  << p.ms << " " << p.n << " "
1209  << psi << " " << mot << std::endl;
1210  } else {
1211  std::cout << "T,m,n,psi,mot: " << T << " " << p.m << " "
1212  << p.n << " " << psi << " " << mot << std::endl;
1213  }
1214  std::cout.precision(6);
1215  std::cout << "mu,dndmu,dndT,dsdT: " << std::endl;
1216  std::cout << "approx: " << p.mu << " "
1217  << p.dndmu << " " << p.dndT << " "
1218  << p.dsdT << std::endl;
1219  std::cout << "exact : " << exact.mu << " "
1220  << exact.dndmu << " " << exact.dndT << " "
1221  << exact.dsdT << std::endl;
1222  std::cout << "bad : " << bad.mu << " " << bad.dndmu << " "
1223  << bad.dndT << " " << bad.dsdT << std::endl;
1224  std::cout << std::endl;
1225  if (verbose>2) {
1226  char ch;
1227  std::cin >> ch;
1228  }
1229  }
1230 
1231  if (ret_local>ret) {
1232  ret=ret_local;
1233  }
1234 
1235  // End of loop over points in data file
1236  }
1237  // End of temperature loop
1238  }
1239 
1240  dev.mu/=cnt;
1241  dev.dndmu/=cnt;
1242  dev.dndT/=cnt;
1243  dev.dsdT/=cnt;
1244 
1245  if (verbose>0) {
1246  if (k==0) {
1247  std::cout << "Function calc_density(), include rest mass:"
1248  << std::endl;
1249  } else if (k==1) {
1250  std::cout << "Function calc_density(), without rest mass:"
1251  << std::endl;
1252  } else if (k==2) {
1253  std::cout << "Function calc_density(), include "
1254  << "rest mass, interacting:"
1255  << std::endl;
1256  } else {
1257  std::cout << "Function calc_density(), without rest mass, "
1258  << "interacting:"
1259  << std::endl;
1260  }
1261 
1262  std::cout << "Average performance: " << std::endl;
1263  std::cout << "dndmu: " << dev.dndmu << " dndT: "
1264  << dev.dndT << " dsdT: " << dev.dsdT << std::endl;
1265  std::cout << "Worst case: " << std::endl;
1266  std::cout << "dndmu: " << bad.dndmu << " dndT: "
1267  << bad.dndT << " dsdT: " << bad.dsdT << std::endl;
1268  std::cout << "mu: " << mu_bad << " m: " << m_bad
1269  << " T: " << T_bad << " mot: " << mot_bad
1270  << "\n\tpsi: " << psi_bad << std::endl;
1271  std::cout << std::endl;
1272  if (verbose>2) {
1273  char ch;
1274  std::cin >> ch;
1275  }
1276  }
1277 
1278  // End of k loop
1279  }
1280 
1281  if (test_pair) {
1282 
1283  // ----------------------------------------------------------------
1284  // Third pass, test pair_mu()
1285 
1286  // k=0,2 are with rest mass, k=1,3 are without
1287  // k=0,1 are non-interacting, k=2,3 are interacting
1288  for(size_t k=0;k<4;k++) {
1289 
1290  double ret_local=0.0;
1291 
1292  // Initialize storage
1293  dev.n=0.0; dev.dndmu=0.0; dev.dndT=0.0; dev.dsdT=0.0;
1294  bad.n=0.0; bad.dndmu=0.0; bad.dndT=0.0; bad.dsdT=0.0;
1295 
1296  // Temperature loop
1297  for(double T=1.0e-2;T<=1.001e2;T*=1.0e2) {
1298 
1299  // Loop over each point in the data file
1300  for(size_t i=0;i<tab.get_nlines();i++) {
1301 
1302  double mot=tab.get("mot",i);
1303  double psi=tab.get("psi",i);
1304  exact.n=tab.get("pair_n",i);
1305  exact.dndmu=tab.get("pair_dndmu",i);
1306  exact.dndT=tab.get("pair_dndT",i);
1307  exact.dsdT=tab.get("pair_dsdT",i);
1308 
1309  set_mass_flags(p,mot,T,k);
1310  set_chem_pot(p,psi,T,k,nr_mode);
1311 
1312  th.pair_mu(p,T);
1313 
1314  exact.n*=pow(T,3.0);
1315  exact.dndmu*=pow(T,2.0);
1316  exact.dndT*=pow(T,2.0);
1317  exact.dsdT*=pow(T,2.0);
1318 
1319  dev.dndmu+=fabs((p.dndmu-exact.dndmu)/exact.dndmu);
1320  dev.dndT+=fabs((p.dndT-exact.dndT)/exact.dndT);
1321  dev.dsdT+=fabs((p.dsdT-exact.dsdT)/exact.dsdT);
1322 
1323  cnt++;
1324 
1325  check_derivs<part_t>(p,exact,bad,k,T,mot,psi,mu_bad,m_bad,
1326  T_bad,mot_bad,psi_bad,ret_local);
1327 
1328  if (verbose>1) {
1329  std::cout.precision(5);
1330  std::cout << "T,m,mu,psi,mot: " << T << " " << p.m
1331  << " " << p.mu
1332  << " " << psi << " " << mot << std::endl;
1333  std::cout.precision(6);
1334  std::cout << "n,dndmu,dndT,dsdT: " << std::endl;
1335  std::cout << "approx: " << p.n << " " << p.dndmu << " "
1336  << p.dndT << " "
1337  << p.dsdT << std::endl;
1338  std::cout << "exact : " << exact.n << " "
1339  << exact.dndmu << " " << exact.dndT << " "
1340  << exact.dsdT << std::endl;
1341  std::cout << "bad : " << bad.n << " " << bad.dndmu << " "
1342  << bad.dndT << " " << bad.dsdT << std::endl;
1343  std::cout << std::endl;
1344  if (verbose>2) {
1345  char ch;
1346  std::cin >> ch;
1347  }
1348  }
1349 
1350  if (ret_local>ret) {
1351  ret=ret_local;
1352  }
1353 
1354  // End of loop over points in data file
1355  }
1356  // End of temperature loop
1357  }
1358 
1359  dev.n/=cnt;
1360  dev.dndmu/=cnt;
1361  dev.dndT/=cnt;
1362  dev.dsdT/=cnt;
1363 
1364  if (verbose>0) {
1365  if (k==0) {
1366  std::cout << "Function pair_mu(), include rest mass"
1367  << std::endl;
1368  } else if (k==1) {
1369  std::cout << "Function pair_mu(), without rest mass"
1370  << std::endl;
1371  } else if (k==2) {
1372  std::cout << "Function pair_mu(), include rest mass, "
1373  << "interacting" << std::endl;
1374  } else {
1375  std::cout << "Function pair_mu(), without rest mass, "
1376  << "interacting" << std::endl;
1377  }
1378 
1379  std::cout << "Average performance: " << std::endl;
1380  std::cout << "dndmu: " << dev.dndmu << " dndT: "
1381  << dev.dndT << " dsdT: " << dev.dsdT << std::endl;
1382  std::cout << "Worst case: " << std::endl;
1383  std::cout << "dndmu: " << bad.dndmu << " dndT: "
1384  << bad.dndT << " dsdT: " << bad.dsdT << std::endl;
1385  std::cout << "mu: " << mu_bad << " m: " << m_bad
1386  << " T: " << T_bad << " mot: " << mot_bad
1387  << "\n\tpsi: " << psi_bad << std::endl;
1388  std::cout << std::endl;
1389  if (verbose>2) {
1390  char ch;
1391  std::cin >> ch;
1392  }
1393  }
1394 
1395  // End of k loop
1396  }
1397 
1398  // ----------------------------------------------------------------
1399  // Fourth pass, test pair_density()
1400 
1401  // k=0,2 are with rest mass, k=1,3 are without
1402  // k=0,1 are non-interacting, k=2,3 are interacting
1403  for(size_t k=0;k<4;k++) {
1404 
1405  double ret_local=0.0;
1406 
1407  // Initialize storage
1408  dev.mu=0.0; dev.dndmu=0.0; dev.dndT=0.0; dev.dsdT=0.0;
1409  bad.mu=0.0; bad.dndmu=0.0; bad.dndT=0.0; bad.dsdT=0.0;
1410 
1411  // Temperature loop
1412  for(double T=1.0e-2;T<=1.001e2;T*=1.0e2) {
1413 
1414  // Loop over each point in the data file
1415  for(size_t i=0;i<tab.get_nlines();i++) {
1416 
1417  double mot=tab.get("mot",i);
1418  double psi=tab.get("psi",i);
1419  p.n=tab.get("pair_n",i);
1420  exact.dndmu=tab.get("pair_dndmu",i);
1421  exact.dndT=tab.get("pair_dndT",i);
1422  exact.dsdT=tab.get("pair_dsdT",i);
1423 
1424  set_mass_flags(p,mot,T,k);
1425  set_chem_pot(exact,psi,T,k,nr_mode);
1426 
1427  p.n*=pow(T,3.0);
1428  exact.dndmu*=pow(T,2.0);
1429  exact.dndT*=pow(T,2.0);
1430  exact.dsdT*=pow(T,2.0);
1431 
1432  // Give it a guess for the chemical potential
1433  p.mu=p.m;
1434 
1435  th.pair_density(p,T);
1436 
1437  dev.dndmu+=fabs((p.dndmu-exact.dndmu)/exact.dndmu);
1438  dev.dndT+=fabs((p.dndT-exact.dndT)/exact.dndT);
1439  dev.dsdT+=fabs((p.dsdT-exact.dsdT)/exact.dsdT);
1440 
1441  cnt++;
1442 
1443  check_derivs<part_t>(p,exact,bad,k,T,mot,psi,mu_bad,m_bad,
1444  T_bad,mot_bad,psi_bad,ret_local);
1445 
1446  if (verbose>1) {
1447  std::cout.precision(5);
1448  std::cout << "T,m,n,psi,mot: " << T << " " << p.m << " "
1449  << p.n << " " << psi << " " << mot << std::endl;
1450  std::cout.precision(6);
1451  std::cout << "mu,dndmu,dndT,dsdT: " << std::endl;
1452  std::cout << "approx: " << p.mu << " " << p.dndmu << " "
1453  << p.dndT << " " << p.dsdT << std::endl;
1454  std::cout << "exact : " << exact.mu << " "
1455  << exact.dndmu << " " << exact.dndT << " "
1456  << exact.dsdT << std::endl;
1457  std::cout << "bad : " << bad.mu << " " << bad.dndmu << " "
1458  << bad.dndT << " " << bad.dsdT << std::endl;
1459  std::cout << std::endl;
1460  if (verbose>2) {
1461  char ch;
1462  std::cin >> ch;
1463  }
1464  }
1465 
1466  if (ret_local>ret) {
1467  ret=ret_local;
1468  }
1469 
1470  // End of loop over points in data file
1471  }
1472  // End of temperature loop
1473  }
1474 
1475  dev.mu/=cnt;
1476  dev.dndmu/=cnt;
1477  dev.dndT/=cnt;
1478  dev.dsdT/=cnt;
1479 
1480  if (verbose>0) {
1481  if (k==0) {
1482  std::cout << "Function pair_density(), include rest mass"
1483  << std::endl;
1484  } else if (k==1) {
1485  std::cout << "Function pair_density(), without rest mass"
1486  << std::endl;
1487  } else if (k==2) {
1488  std::cout << "Function pair_density(), include rest mass, "
1489  << "interacting" << std::endl;
1490  } else {
1491  std::cout << "Function pair_density(), without rest mass, "
1492  << "interacting" << std::endl;
1493  }
1494 
1495  std::cout << "Average performance: " << std::endl;
1496  std::cout << "dndmu: "
1497  << dev.dndmu << " dndT: "
1498  << dev.dndT << " dsdT: " << dev.dsdT << std::endl;
1499  std::cout << "Worst case: " << std::endl;
1500  std::cout << "dndmu: " << bad.dndmu
1501  << " dndT: " << bad.dndT
1502  << " dsdT: " << bad.dsdT << std::endl;
1503  std::cout << "mu: " << mu_bad << " m: " << m_bad
1504  << " T: " << T_bad << " mot: " << mot_bad
1505  << "\n\tpsi: " << psi_bad << std::endl;
1506  std::cout << std::endl;
1507  if (verbose>2) {
1508  char ch;
1509  std::cin >> ch;
1510  }
1511  }
1512 
1513  // End of k loop
1514  }
1515 
1516  // End of 'if (test_pair)'
1517  }
1518 
1519  // ----------------------------------------------------------------
1520  // Return to the original values
1521 
1522  p=orig;
1523 
1524  return ret;
1525  }
1526 
1527  };
1528 
1529 #ifndef DOXYGEN_NO_O2NS
1530 }
1531 #endif
1532 
1533 #endif
o2scl::fermion_deriv_thermo_tl::nu_from_n
virtual int nu_from_n(fermion_deriv &f, fp_t temper)=0
Calculate effective chemical potential from density.
o2scl::part_tl
Particle base class.
Definition: part.h:103
o2scl::part_deriv_press_tl::dsdT
fp_t dsdT
Derivative of entropy density with respect to temperature.
Definition: part_deriv.h:135
o2scl::table::function_column
void function_column(std::string function, std::string scol)
o2scl::fermion_thermo_tl::calc_mu_deg_tlate
bool calc_mu_deg_tlate(fermion_t &f, fp_t temper, fp_t prec)
Calculate thermodynamic properties from the chemical potential using a degenerate expansion.
Definition: fermion.h:446
o2scl::thermo_np_deriv_helm_tl::dmudn_mixed
fp_t dmudn_mixed
The quantity .
Definition: part_deriv.h:84
o2scl::table
o2scl::fermion_deriv_thermo_tl
Compute properties of a fermion including derivatives [abstract base].
Definition: part_deriv.h:712
o2scl::table::get_nlines
size_t get_nlines() const
o2scl::thermo_np_deriv_helm_tl::dmundnn
fp_t dmundnn
The quantity .
Definition: part_deriv.h:82
o2scl::thermo_np_deriv_press_tl
Object to store second derivatives of .
Definition: part_deriv.h:46
o2scl::deriv_thermo_base_tl::squared_sound_speed
fp_t squared_sound_speed(part_deriv_t &p, fp_t temper)
The squared sound speed (unitless)
Definition: part_deriv.h:683
o2scl::fermion_deriv_thermo_tl::pair_density
virtual int pair_density(fermion_deriv &f, fp_t temper)=0
Calculate properties with antiparticles as function of density.
o2scl::part_deriv_press_tl::dndT
fp_t dndT
Derivative of number density with respect to temperature.
Definition: part_deriv.h:132
o2scl_const::pi2
const double pi2
o2scl::part_tl< double >::mu
double mu
Chemical potential.
Definition: part.h:118
o2scl::part_tl< double >::ms
double ms
Effective mass (Dirac unless otherwise specified)
Definition: part.h:122
o2scl::deriv_thermo_base_tl::heat_cap_ppart_const_press
fp_t heat_cap_ppart_const_press(part_deriv_t &p, fp_t temper)
The heat capacity per particle at constant pressure (unitless)
Definition: part_deriv.h:465
o2scl::thermo_np_deriv_helm_tl
Object to store second derivatives of .
Definition: part_deriv.h:71
o2scl::part_tl< double >::n
double n
Number density.
Definition: part.h:112
o2scl::part_deriv_tl::part_deriv_tl
part_deriv_tl(fp_t mass=0.0, fp_t dof=0.0)
Make a particle of mass mass and degeneracy dof.
Definition: part_deriv.h:276
o2scl::part_deriv_calibrate_class
Object to organize calibration of derivative quantities in particle classes to results stored in a ta...
Definition: part_deriv.h:958
o2scl::fermion_deriv_thermo_tl::calc_mu_ndeg
virtual bool calc_mu_ndeg(fermion_deriv &f, fp_t temper, fp_t prec, bool inc_antip=false)
Calculate properties as a function of chemical potential using a nondegenerate expansion.
Definition: part_deriv.h:837
o2scl_const::pi
const double pi
o2scl::fermion_deriv_thermo_tl::fr
fermion_rel fr
A fermion_thermo object.
Definition: part_deriv.h:720
o2scl::fermion_deriv_tl::operator=
fermion_deriv_tl & operator=(const fermion &p)
Copy construction with operator=()
Definition: part_deriv.h:243
o2scl::part_tl< double >::non_interacting
bool non_interacting
True if the particle is non-interacting (default true)
Definition: part.h:130
o2scl::fermion_thermo_tl::calc_mu_ndeg_tlate
bool calc_mu_ndeg_tlate(fermion_t &f, fp_t temper, fp_t prec, bool inc_antip)
Calculate thermodynamic properties from the chemical potential using a nondegenerate expansion.
Definition: fermion.h:324
o2scl::fermion_deriv_tl::fermion_deriv_tl
fermion_deriv_tl(fp_t mass=0.0, fp_t dof=0.0)
Make a particle of mass mass and degeneracy dof.
Definition: part_deriv.h:182
o2scl::thermo_np_deriv_helm_tl::dmupdT
fp_t dmupdT
The quantity .
Definition: part_deriv.h:80
o2scl::part_deriv_tl::operator=
part_deriv_tl & operator=(const part_deriv_tl &p)
Copy construction with operator=()
Definition: part_deriv.h:298
o2scl::part_tl< double >::en
double en
Entropy density.
Definition: part.h:120
o2scl::thermo_np_deriv_press_tl::dnndT
fp_t dnndT
The quantity .
Definition: part_deriv.h:53
o2scl::part_calibrate_class::set_mass_flags
void set_mass_flags(part_t &p, double mot, double T, size_t k)
Set mass and flags from mot, T, and the index k.
Definition: part.h:248
o2scl::deriv_thermo_base_tl
Base quantities for thermodynamic derivatives.
Definition: part_deriv.h:381
o2scl_hdf::hdf_file::close
void close()
hdf_input
void hdf_input(hdf_file &hf, o2scl::table< vec_t > &t, std::string name)
o2scl::fermion_deriv_thermo_tl::pair_mu
virtual int pair_mu(fermion_deriv &f, fp_t temper)=0
Calculate properties with antiparticles as function of chemical potential.
o2scl::fermion_deriv_thermo_tl::pi
fp_t pi
Desc.
Definition: part_deriv.h:723
o2scl::part_tl< double >::g
double g
Degeneracy (e.g. spin and color if applicable)
Definition: part.h:108
o2scl::part_tl< double >::ed
double ed
Energy density.
Definition: part.h:114
o2scl::thermo_np_deriv_press_tl::dnpdT
fp_t dnpdT
The quantity .
Definition: part_deriv.h:55
o2scl_hdf::hdf_file::open
int open(std::string fname, bool write_access=false, bool err_on_fail=true)
o2scl::part_deriv_tl::part_deriv_tl
part_deriv_tl(const part_deriv_tl &p)
Copy constructor.
Definition: part_deriv.h:280
o2scl::fermion_rel_tl< double >
o2scl_settings
lib_settings_class o2scl_settings
o2scl::part_deriv_press_tl::part_deriv_press_tl
part_deriv_press_tl(const part_deriv_press_tl &p)
Copy constructor.
Definition: part_deriv.h:141
o2scl::deriv_thermo_base_tl::compress_const_tptr
fp_t compress_const_tptr(part_deriv_t &p, fp_t temper)
The isothermal compressibility.
Definition: part_deriv.h:549
o2scl::fermion_tl
Fermion class.
Definition: fermion.h:55
o2scl::thermo_np_deriv_press_tl::dsdT
fp_t dsdT
The quantity .
Definition: part_deriv.h:51
o2scl::fermion_deriv_tl
A fermion with derivative information.
Definition: part_deriv.h:176
o2scl::part_deriv_press_tl::operator=
part_deriv_press_tl & operator=(const part_deriv_press_tl &p)
Copy construction with operator=()
Definition: part_deriv.h:148
o2scl::part_tl< double >::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_deriv_tl::fermion_deriv_tl
fermion_deriv_tl(const fermion_deriv_tl &p)
Copy constructor.
Definition: part_deriv.h:186
o2scl::part_deriv_press_tl::dndmu
fp_t dndmu
Derivative of number density with respect to chemical potential.
Definition: part_deriv.h:129
o2scl::fermion_deriv_thermo_tl::pi2
fp_t pi2
Desc.
Definition: part_deriv.h:726
o2scl::fermion_deriv_tl::operator=
fermion_deriv_tl & operator=(const fermion_deriv_tl &p)
Copy construction with operator=()
Definition: part_deriv.h:222
o2scl::thermo_np_deriv_press_tl::dnndmun
fp_t dnndmun
The quantity .
Definition: part_deriv.h:57
o2scl::part_tl< double >::m
double m
Mass.
Definition: part.h:110
o2scl::part_tl< double >::pr
double pr
Pressure.
Definition: part.h:116
o2scl::thermo_np_deriv_helm_tl::dmundT
fp_t dmundT
The quantity .
Definition: part_deriv.h:78
o2scl::fermion_deriv_thermo_tl::calc_mu_deg
virtual bool calc_mu_deg(fermion_deriv &f, fp_t temper, fp_t prec)
Calculate properties as a function of chemical potential using a degenerate expansion.
Definition: part_deriv.h:767
exc_efilenotfound
exc_efilenotfound
o2scl::thermo_np_deriv_press_tl::dndmu_mixed
fp_t dndmu_mixed
The quantity .
Definition: part_deriv.h:59
o2scl::part_deriv_tl
A part with derivative information.
Definition: part_deriv.h:270
o2scl_hdf::hdf_file
o2scl::deriv_thermo_base_tl::heat_cap_ppart_const_vol
fp_t heat_cap_ppart_const_vol(part_deriv_t &p, fp_t temper)
The heat capacity per particle at constant volume (unitless)
Definition: part_deriv.h:416
o2scl::part_deriv_calibrate_class::part_deriv_calibrate
double part_deriv_calibrate(part_t &p, thermo_t &th, bool test_pair, std::string file, bool nr_mode=false, int verbose=0, bool external=false)
Calibrate a particle thermodynamics class with derivatives with results stored in a table.
Definition: part_deriv.h:975
o2scl::fermion_deriv_thermo_tl::calc_mu
virtual int calc_mu(fermion_deriv &f, fp_t temper)=0
Calculate properties as function of chemical potential.
O2SCL_ERR
#define O2SCL_ERR(d, n)
o2scl::table::is_column
bool is_column(std::string scol) const
o2scl::part_deriv_press_tl::deriv_f
void deriv_f(fp_t &dmudn, fp_t &dmudT, fp_t &dsdT_n)
Compute derivatives in the Helmholtz free energy representation from the derivatives in the pressure ...
Definition: part_deriv.h:161
o2scl::fermion_thermo_tl::ndeg_terms
void ndeg_terms(size_t j, fp_t tt, fp_t xx, fp_t m, bool inc_rest_mass, bool inc_antip, fp_t &pterm, fp_t &nterm, fp_t &enterm)
Compute a term in the nondegenerate expansion.
Definition: fermion.h:838
o2scl::part_calibrate_class::set_chem_pot
void set_chem_pot(part_t &p, double psi, double T, size_t k, bool nr_mode)
Set chemical potential from psi, T, the index k, and the flag nr_mode.
Definition: part.h:270
o2scl::thermo_np_deriv_press_tl::dnpdmup
fp_t dnpdmup
The quantity .
Definition: part_deriv.h:61
o2scl::fermion_deriv_tl::fermion_deriv_tl
fermion_deriv_tl(const fermion &p)
Copy constructor.
Definition: part_deriv.h:204
o2scl::deriv_thermo_base_tl::compress_adiabatic
fp_t compress_adiabatic(part_deriv_t &p, fp_t temper)
The adiabatic compressibility.
Definition: part_deriv.h:513
o2scl::table::get
double get(std::string scol, size_t row) const
o2scl::part_tl< double >::nu
double nu
Effective chemical potential.
Definition: part.h:124
o2scl::part_calibrate_class
Object to organize calibration of particle classes to results stored in a table.
Definition: part.h:241
psi
const double psi
o2scl::fermion_deriv_thermo_tl::calc_density
virtual int calc_density(fermion_deriv &f, fp_t temper)=0
Calculate properties as function of density.
o2scl::deriv_thermo_base_tl::coeff_thermal_exp
fp_t coeff_thermal_exp(part_deriv_t &p, fp_t temper)
The coefficient of thermal expansion.
Definition: part_deriv.h:590
o2scl::thermo_np_deriv_helm_tl::dmupdnp
fp_t dmupdnp
The quantity .
Definition: part_deriv.h:86
o2scl::thermo_np_deriv_helm_tl::dsdT
fp_t dsdT
The quantity .
Definition: part_deriv.h:76
o2scl::part_deriv_press_tl
Particle derivatives in the pressure representation.
Definition: part_deriv.h:124

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