eos_tov.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2017, 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_tov.h
24  \brief File defining \ref o2scl::eos_tov
25 */
26 #ifndef O2SCL_TOV_EOS_H
27 #define O2SCL_TOV_EOS_H
28 
29 #include <cmath>
30 #include <iostream>
31 #include <fstream>
32 
33 #include <boost/numeric/ublas/vector.hpp>
34 
35 #include <o2scl/constants.h>
36 #include <o2scl/lib_settings.h>
37 #include <o2scl/interp.h>
38 #include <o2scl/table_units.h>
39 #include <o2scl/vector_derint.h>
40 
41 #ifndef DOXYGEN_NO_O2NS
42 namespace o2scl {
43 #endif
44 
45  /** \brief A EOS base class for the TOV solver
46  */
47  class eos_tov {
48 
49  protected:
50 
51  /** \brief Set to true if the baryon density is provided in the
52  EOS (default false)
53  */
55 
56  public:
57 
58  eos_tov();
59 
60  virtual ~eos_tov() {}
61 
62  /// Control for output (default 1)
63  int verbose;
64 
65  /// Return true if a baryon density is available
66  bool has_baryons() {
67  return baryon_column;
68  }
69 
70  /** \brief Check that the baryon density is consistent
71  with the \f$ P(\varepsilon) \f$
72  */
73  void check_nb(double &avg_abs_dev, double &max_abs_dev);
74 
75  /** \brief From the pressure, return the energy density
76  */
77  virtual double ed_from_pr(double pr)=0;
78 
79  /** \brief From the energy density, return the pressure
80  */
81  virtual double pr_from_ed(double ed)=0;
82 
83  /** \brief From the energy density, return the baryon density
84  */
85  virtual double nb_from_ed(double ed)=0;
86 
87  /** \brief From the pressure, return the baryon density
88  */
89  virtual double nb_from_pr(double pr)=0;
90 
91  /** \brief From the baryon density, return the energy density
92  */
93  virtual double ed_from_nb(double nb)=0;
94 
95  /** \brief From the baryon density, return the pressure
96  */
97  virtual double pr_from_nb(double nb)=0;
98 
99  /** \brief Given the pressure, produce the energy and number densities
100 
101  The arguments \c pr and \c ed should always be in \f$
102  M_{\odot}/\mathrm{km}^3 \f$ . The argument for \c nb should be
103  in \f$ \mathrm{fm}^{-3} \f$ .
104 
105  If \ref baryon_column is false, then \c nb is unmodified.
106  */
107  virtual void ed_nb_from_pr(double pr, double &ed, double &nb)=0;
108 
109  };
110 
111  /** \brief The Buchdahl EOS for the TOV solver
112 
113  Given the EOS
114  \f[
115  \varepsilon = 12 \sqrt{p_{*} P}- 5 P
116  \f]
117  the TOV equation has an analytical solution
118  \f[
119  R=(1-\beta) \sqrt{\frac{\pi}{288 p_{*} G (1-2 \beta)}}
120  \f]
121  where \f$ \beta = GM/R \f$.
122 
123  The baryon chemical potential is
124  \f[
125  \mu = \mu_1 \left[ \frac{\left(9 p_{*}-P_1\right) \left(3+t_1\right)
126  \left(3-t_2\right)}{\left(9 p_{*}-P\right)\left(3-t_1\right)
127  \left(3+t_2\right)}\right]^{1/4}
128  \f]
129  where \f$ t_1 = \sqrt{P}/\sqrt{p_{*}} \f$ and \f$ t_2 =
130  \sqrt{P_1/p_{*}} \f$ . The baryon density can then be obtained
131  directly from the thermodynamic identity. In the case that one
132  assumes \f$ \mu_1 = m_n \f$ and \f$ P_1 = 0 \f$, the baryon
133  density can be simplified to
134  \f[
135  n m_n = 12 \sqrt{P p_{*}} \left( 1-\frac{1}{3} \sqrt{P/p_{*}}
136  \right)^{3/2}
137  \f]
138  c.f. Eq. 10 in \ref Lattimer01.
139 
140  The central pressure and energy density are
141  \f[
142  P_c = 36 p_{*} \beta^2
143  \f]
144  \f[
145  {\varepsilon}_c = 72 p_{*} \beta (1-5 \beta/2)
146  \f]
147 
148  Physical solutions are obtained only for \f$ P< 25 p_{*}/144 \f$
149  (ensuring that the argument to the square root is positive)
150  and \f$ \beta<1/6 \f$ (ensuring that the EOS is not acausal).
151 
152  Based on \ref Lattimer01 .
153 
154  \todo Fix base EOS functions
155  \future Figure out what to do with the buchfun() function
156  */
157  class eos_tov_buchdahl : public eos_tov {
158 
159  protected:
160 
161  /** \brief The baryon density at \c ed1
162  */
163  double nb1;
164 
165  /** \brief The energy density for which the baryon density is known
166  */
167  double ed1;
168 
169  /** \brief The pressure at \c ed1
170  */
171  double pr1;
172 
173  public:
174 
176 
177  virtual ~eos_tov_buchdahl() {}
178 
179  /** \brief The parameter with units of pressure in units of solar
180  masses per km cubed (default value \f$ 3.2 \times 10^{-5} \f$
181  )
182  */
183  double Pstar;
184 
185  /** \brief Set the baryon density
186  */
187  void set_baryon_density(double nb, double ed);
188 
189  /** \brief From the pressure, return the energy density
190  */
191  virtual double ed_from_pr(double pr);
192 
193  /** \brief From the energy density, return the pressure
194  */
195  virtual double pr_from_ed(double ed);
196 
197  /** \brief From the energy density, return the baryon density
198  */
199  virtual double nb_from_ed(double ed);
200 
201  /** \brief From the pressure, return the baryon density
202  */
203  virtual double nb_from_pr(double pr);
204 
205  /** \brief From the baryon density, return the energy density
206  */
207  virtual double ed_from_nb(double nb);
208 
209  /** \brief From the baryon density, return the pressure
210  */
211  virtual double pr_from_nb(double nb);
212 
213  /** \brief Given the pressure, produce the energy and number densities
214 
215  If the baryon density is not specified, it should be set to
216  zero or \ref baryon_column should be set to false
217  */
218  virtual void ed_nb_from_pr(double pr, double &ed, double &nb);
219 
220  protected:
221 
222  /** \brief Solve to compute profiles
223 
224  After solving the two equations
225  \f{eqnarray*}
226  r^{\prime} &=& r \left(1-\beta+u\right)^{-1}
227  \left(1 - 2 \beta\right) \nonumber \\
228  A^2 &=& 288 \pi p_{*} G \left( 1 - 2 \beta \right)^{-1}
229  \f}
230  for \f$ u = \beta/(A r^{\prime}) \sin A r^{\prime} \f$
231  and \f$ r^{\prime} \f$,
232  one can compute the pressure and energy density
233  profiles
234  \f{eqnarray*}
235  8 \pi P &=& A^2 u^2 \left(1 - 2 \beta \right)
236  \left(1 - \beta + u \right)^{-2}
237  \nonumber \\
238  8 \pi \varepsilon &=& 2 A^2 u \left(1 - 2 \beta\right)
239  \left(1 - \beta - 3 u/2\right) \left(1 - \beta + u \right)^{-2}
240  \nonumber \\
241  \f}
242  */
243  int solve_u_rp_fun(size_t bv, const std::vector<double> &bx,
244  std::vector<double> &by);
245 
246  };
247 
248  /** \brief Standard polytropic EOS \f$ P = K \varepsilon^{1+1/n} \f$
249 
250  The quantity \f$ K \f$ must be in units of
251  \f$ \left(M_{\odot}/km^3\right)^{-1/n} \f$ .
252 
253  \comment
254  The documentation below was taken from bamr.
255  \endcomment
256  For a polytrope \f$ P = K \varepsilon^{1+1/n} \f$
257  beginning at a pressure of \f$ P_1 \f$, an energy
258  density of \f$ \varepsilon_1 \f$ and a baryon density
259  of \f$ n_{B,1} \f$, the baryon density along the polytrope
260  is
261  \f[
262  n_B = n_{B,1} \left(\frac{\varepsilon}{\varepsilon_1}\right)^{1+n}
263  \left(\frac{\varepsilon_1+P_1}{\varepsilon+P}\right)^{n} \, .
264  \f]
265  Similarly, the chemical potential is
266  \f[
267  \mu_B = \mu_{B,1} \left(1 + \frac{P_1}{\varepsilon_1}\right)^{1+n}
268  \left(1 + \frac{P}{\varepsilon}\right)^{-(1+n)} \, .
269  \f]
270  The expression for the
271  baryon density can be inverted to determine \f$ \varepsilon(n_B) \f$
272  \f[
273  \varepsilon(n_B) = \left[ \left(\frac{n_{B,1}}
274  {n_B \varepsilon_1} \right)^{1/n}
275  \left(1+\frac{P_1}{\varepsilon_1}\right)-K\right]^{-n} \, .
276  \f]
277  Sometimes the baryon susceptibility is also useful
278  \f[
279  \frac{d \mu_B}{d n_B} = \left(1+1/n\right)
280  \left( \frac{P}{\varepsilon}\right)
281  \left( \frac{\mu_B}{n_B}\right) \, .
282  \f]
283 
284  \future The simple formulation fo the expressions here more than
285  likely break down when their arguments are sufficiently extreme.
286  */
287  class eos_tov_polytrope : public eos_tov {
288 
289  protected:
290 
291  /** \brief The baryon density at \c ed1
292  */
293  double nb1;
294 
295  /** \brief The energy density for which the baryon density is known
296  */
297  double ed1;
298 
299  /** \brief The pressure at \c ed1
300  */
301  double pr1;
302 
303  /** \brief Coefficient (default 1.0)
304  */
305  double K;
306 
307  /// Index (default 3.0)
308  double n;
309 
310  public:
311 
313 
314  virtual ~eos_tov_polytrope() {}
315 
316  /** \brief Set the coefficient and polytropic index
317  */
318  void set_coeff_index(double coeff, double index);
319 
320  /** \brief Set the baryon density
321  */
322  void set_baryon_density(double nb, double ed);
323 
324  /** \brief From the pressure, return the energy density
325  */
326  virtual double ed_from_pr(double pr);
327 
328  /** \brief From the energy density, return the pressure
329  */
330  virtual double pr_from_ed(double ed);
331 
332  /** \brief From the energy density, return the baryon density
333  */
334  virtual double nb_from_ed(double ed);
335 
336  /** \brief From the pressure, return the baryon density
337  */
338  virtual double nb_from_pr(double pr);
339 
340  /** \brief From the baryon density, return the energy density
341  */
342  virtual double ed_from_nb(double nb);
343 
344  /** \brief From the baryon density, return the pressure
345  */
346  virtual double pr_from_nb(double nb);
347 
348  /** \brief Given the pressure, produce the energy and number densities
349  */
350  virtual void ed_nb_from_pr(double pr, double &ed, double &nb);
351 
352  };
353 
354  /** \brief Linear EOS \f$ P = c_s^2 (\varepsilon-\varepsilon_0) \f$
355 
356  This implements a linear EOS with a fixed speed of sound and a
357  fixed energy density at zero pressure. This will also compute
358  the baryon density, if one calls \ref set_baryon_density() to
359  set the baryon density at one fiducial energy density.
360 
361  Given a fiducial baryon density \f$ n_{B,1} \f$ at
362  some energy density \f$ \varepsilon_1 \f$ and pressure
363  \f$ P_1 \f$, the baryon density is
364  \f[
365  n_B = n_{B,1} \left[ \frac{\varepsilon(1+c_s^2) -
366  c_s^2 \varepsilon_0 } {\varepsilon_1 (1 + c_s^2) -
367  c_s^2 \varepsilon_0}\right]^{1/(1+c_s^2)} =
368  n_{B,1} \left[ \frac{ \varepsilon + P }
369  {\varepsilon_1 + P_1} \right]^{1/(1+c_s^2)}
370  \f]
371 
372  \note Experimental
373  */
374  class eos_tov_linear : public eos_tov {
375 
376  protected:
377 
378  /** \brief The baryon density at \c ed1
379  */
380  double nb1;
381 
382  /** \brief The energy density for which the baryon density is known
383  */
384  double ed1;
385 
386  /** \brief The pressure for which the baryon density is known
387  */
388  double pr1;
389 
390  /** \brief Coefficient (default 1.0)
391  */
392  double cs2;
393 
394  /// The energy density at zero pressure (default 0.0)
395  double eps0;
396 
397  public:
398 
399  eos_tov_linear();
400 
401  virtual ~eos_tov_linear() {}
402 
403  /** \brief Set the sound speed and energy density at zero pressure
404  */
405  void set_cs2_eps0(double cs2_, double eps0_);
406 
407  /** \brief Set the baryon density
408  */
409  void set_baryon_density(double nb, double ed);
410 
411  /** \brief From the pressure, return the energy density
412  */
413  virtual double ed_from_pr(double pr);
414 
415  /** \brief From the energy density, return the pressure
416  */
417  virtual double pr_from_ed(double ed);
418 
419  /** \brief From the energy density, return the baryon density
420  */
421  virtual double nb_from_ed(double ed);
422 
423  /** \brief From the pressure, return the baryon density
424  */
425  virtual double nb_from_pr(double pr);
426 
427  /** \brief From the baryon density, return the energy density
428  */
429  virtual double ed_from_nb(double nb);
430 
431  /** \brief From the baryon density, return the pressure
432  */
433  virtual double pr_from_nb(double nb);
434 
435  /** \brief Given the pressure, produce the energy and number densities
436  */
437  virtual void ed_nb_from_pr(double pr, double &ed, double &nb);
438 
439  };
440 
441  /** \brief Provide an EOS for TOV solvers based on
442  interpolation of user-supplied vectors
443  */
444  template<class vec_t> class eos_tov_vectors : public eos_tov {
445 
446  /** \brief Desc
447  */
448  void reset_interp(size_t n) {
449  pe_int.set(n,pr,ed,itp_linear);
450  ep_int.set(n,ed,pr,itp_linear);
451  return;
452  }
453 
454  /** \brief Desc
455  */
456  void reset_interp_nb(size_t n) {
457  reset_interp();
458  pn_int.set(n,pr,nb,itp_linear);
459  np_int.set(n,nb,pr,itp_linear);
460  en_int.set(n,ed,nb,itp_linear);
461  ne_int.set(n,nb,ed,itp_linear);
462  return;
463  }
464 
465  public:
466 
467  /** \brief Read the EOS from a set of equal length
468  vectors for energy density, pressure, and baryon density
469 
470  In this version, the user-specified vectors are swapped
471  with internal storage.
472  */
473  void read_vectors_swap(size_t user_n, vec_t &user_ed, vec_t &user_pr,
474  vec_t &user_nb) {
475  std::swap(user_ed,ed);
476  std::swap(user_pr,pr);
477  std::swap(user_nb,nb);
478  this->baryon_column=true;
479  reset_interp_nb(user_n);
480  return;
481  }
482 
483  /** \brief Read the EOS from a pair of equal length
484  vectors for energy density and pressure
485 
486  In this version, the user-specified vectors are swapped
487  with internal storage.
488  */
489  void read_vectors_swap(size_t user_n, vec_t &user_ed, vec_t &user_pr) {
490  std::swap(user_ed,ed);
491  std::swap(user_pr,pr);
492  this->baryon_column=false;
493  reset_interp(user_n);
494  return;
495  }
496 
497  /** \brief Read the EOS from a set of equal length
498  vectors for energy density, pressure, and baryon density
499 
500  In this version, the user-specified vectors are copied
501  to internal storage.
502  */
503  void read_vectors_copy(size_t user_n, vec_t &user_ed, vec_t &user_pr,
504  vec_t &user_nb) {
505  if (ed.size()!=user_n) ed.resize(user_n);
506  if (pr.size()!=user_n) pr.resize(user_n);
507  if (nb.size()!=user_n) nb.resize(user_n);
508  vector_copy(user_ed,ed);
509  vector_copy(user_pr,pr);
510  vector_copy(user_nb,nb);
511  this->baryon_column=true;
512  reset_interp_nb(user_n);
513  return;
514  }
515 
516  /** \brief Read the EOS from a pair of equal length
517  vectors for energy density and pressure
518 
519  In this version, the user-specified vectors are copied
520  to internal storage.
521  */
522  void read_vectors_copy(size_t user_n, vec_t &user_ed, vec_t &user_pr) {
523  if (ed.size()!=user_n) ed.resize(user_n);
524  if (pr.size()!=user_n) pr.resize(user_n);
525  vector_copy(user_ed,ed);
526  vector_copy(user_pr,pr);
527  this->baryon_column=false;
528  reset_interp(user_n);
529  return;
530  }
531 
532  /// \name Basic EOS functions
533  //@{
534  /** \brief From the pressure, return the energy density
535  */
536  virtual double ed_from_pr(double pr) {
537  return pe_int.eval(pr);
538  }
539 
540  /** \brief From the energy density, return the pressure
541  */
542  virtual double pr_from_ed(double ed) {
543  return ep_int.eval(pr);
544  }
545 
546  /** \brief From the energy density, return the baryon density
547  */
548  virtual double nb_from_ed(double ed) {
549  return en_int.eval(ed);
550  }
551 
552  /** \brief From the pressure, return the baryon density
553  */
554  virtual double nb_from_pr(double pr) {
555  return pn_int.eval(ed);
556  }
557 
558  /** \brief From the baryon density, return the energy density
559  */
560  virtual double ed_from_nb(double nb) {
561  return ne_int.eval(ed);
562  }
563 
564  /** \brief From the baryon density, return the pressure
565  */
566  virtual double pr_from_nb(double nb) {
567  return np_int.eval(ed);
568  }
569 
570  /** \brief Given the pressure, produce the energy and number densities
571 
572  The arguments \c pr and \c ed should always be in \f$
573  M_{\odot}/\mathrm{km}^3 \f$ . The argument for \c nb should be
574  in \f$ \mathrm{fm}^{-3} \f$ .
575 
576  If \ref baryon_column is false, then \c nb is unmodified.
577  */
578  virtual void ed_nb_from_pr(double pr, double &ed, double &nb) {
579  ed_from_pr(pr);
580  if (this->baryon_column) {
581  nb_from_pr(pr);
582  }
583  return;
584  }
585  //@}
586 
587  protected:
588 
589  /// \name EOS storage
590  //@{
591  /// Energy densities from full EOS
592  vec_t ed;
593  /// Pressures from full EOS
594  vec_t pr;
595  /// Baryon densities from full EOS
596  vec_t nb;
597  //@}
598 
599  /// \name Interpolators
600  //@{
601  interp<vec_t> pe_int;
602  interp<vec_t> pn_int;
603  interp<vec_t> ep_int;
604  interp<vec_t> en_int;
605  interp<vec_t> np_int;
606  interp<vec_t> ne_int;
607  //@}
608 
609  };
610 
611  /** \brief An EOS for the TOV solver using simple linear
612  interpolation and an optional crust EOS
613 
614  The simplest usage of this class is simply to use \ref
615  read_table() to read a tabulated EOS stored in a \ref
616  table_units object and optionally specify a separate crust EOS.
617 
618  Alternatively, the user can simply specify objects
619  of type <tt>std::vector<double></tt> which store the energy
620  density, pressure, and baryon density.
621 
622  There are two methods to handle the crust-core interface. The
623  default, <tt>smooth_trans</tt> uses the crust below pressure \f$
624  P_1 \f$ (equal to the value of \ref trans_pres divided by \ref
625  trans_width) and the core above pressure \f$ P_2 \f$ (the value
626  of \ref trans_pres times \ref trans_width) and then in between
627  uses
628  \f[
629  \varepsilon(P) = [1-\chi(P)] \varepsilon_{\mathrm{crust}}(P) +
630  \chi(P) \varepsilon_{\mathrm{core}}(P)
631  \f]
632  where the value \f$ \chi(P) \f$ is determined by
633  \f[
634  \chi(P) = (P-P_1)/(P_2-P_1) \, .
635  \f]
636  This method is a bit more faithful to the original EOS tables,
637  but the matching can result in pressures which decrease with
638  increasing energy density. Alternatively the <tt>match_line</tt>
639  method uses
640  \f$ \varepsilon_1=\varepsilon_{\mathrm{crust}}(P_1) \f$ and
641  \f$ \varepsilon_2=\varepsilon_{\mathrm{core}}(P_2) \f$ and
642  \f[
643  \varepsilon(P) = (\varepsilon_2 - \varepsilon_1) \chi
644  + \varepsilon_1 \, .
645  \f]
646  (using the same expression for \f$ \chi \f$ ). This method less
647  frequently results in decreasing pressures, but can deviate
648  further from the original tables.
649 
650  Internally, energy and pressure are stored in units of \f$
651  \mathrm{M}_{\odot}/\mathrm{km}^3 \f$ and baryon density is
652  stored in units of \f$ \mathrm{fm}^{-3} \f$ . The user-specified
653  EOS table is left as is, and unit conversion is performed as
654  needed in ed_nb_from_pr() and other functions from the units
655  specified in the input \ref table_units object.
656 
657  \comment
658  \todo It might be useful to exit more gracefully when non-finite
659  values are obtained in interpolation, analogous to the
660  err_nonconv mechanism elsewhere.
661  2/6/16: I'm not sure this is really necessary. It's only
662  linear interpolation, so the err_nonconv mechanism probably
663  won't be useful.
664  \endcomment
665 
666  \future Create a sanity check where core_auxp is nonzero only if
667  core_table is also nonzero. Alternatively, this complication is
668  due to the fact that this class works in two ways: one where it
669  reads a table (and adds a crust), and one where it reads in
670  vectors (with no crust). Maybe these two modes of operation
671  should be separated into two classes. (2/6/16: Now there is
672  a new eos_tov_vector class. The best way forward may be
673  to make this a child of eos_tov_vectors.)
674  */
675  class eos_tov_interp : public eos_tov {
676 
677  public:
678 
679  eos_tov_interp();
680 
681  virtual ~eos_tov_interp();
682 
683  /// \name Mode of transitioning between crust and core EOS
684  //@{
685  int transition_mode;
686  static const int smooth_trans=0;
687  static const int match_line=1;
688  //@}
689 
690  /// \name Basic EOS functions
691  //@{
692  /** \brief From the pressure, return the energy density
693  */
694  virtual double ed_from_pr(double pr);
695 
696  /** \brief From the energy density, return the pressure
697  */
698  virtual double pr_from_ed(double ed);
699 
700  /** \brief From the energy density, return the baryon density
701  */
702  virtual double nb_from_ed(double ed);
703 
704  /** \brief From the pressure, return the baryon density
705  */
706  virtual double nb_from_pr(double pr);
707 
708  /** \brief From the baryon density, return the energy density
709  */
710  virtual double ed_from_nb(double nb);
711 
712  /** \brief From the baryon density, return the pressure
713  */
714  virtual double pr_from_nb(double nb);
715  //@}
716 
717  /// \name Basic usage
718  //@{
719  /** \brief Specify the EOS through a table
720 
721  If units are specified for any of the columns, then this
722  function attempts to automatically determine the correct
723  conversion factors using the \ref o2scl::convert_units object
724  returned by \ref o2scl::o2scl_settings . If the units for any
725  of the columns are blank, then they are assumed to be the
726  native units for \ref o2scl::tov_solve .
727 
728  \note The input table must have at least 2 rows and
729  the pressure column must be strictly increasing.
730 
731  This function copies the needed information from the
732  table so if it is modified then this function
733  needs to be called again to read a new table.
734  */
735  void read_table(table_units<> &eosat, std::string s_cole,
736  std::string s_colp, std::string s_colnb="");
737  //@}
738 
739  /// \name Crust EOS functions
740  //@{
741  /// Default crust EOS from \ref Negele73 and \ref Baym71tg
742  void default_low_dens_eos();
743 
744  /// Crust EOS from \ref Shen11b
745  void sho11_low_dens_eos();
746 
747  /** \brief Crust EOS from \ref Steiner12
748 
749  This function uses the neutron star crust models from \ref
750  Steiner12 . The current acceptable values for \c model are
751  <tt>APR</tt>, <tt>Gs</tt>, <tt>Rs</tt> and <tt>SLy4</tt>.
752  */
753  void s12_low_dens_eos(std::string model="SLy4",
754  bool external=false);
755 
756  /** \brief Crust EOS from Goriely, Chamel, and Pearson
757 
758  From \ref Goriely10, \ref Pearson11, and \ref Pearson12 .
759  */
760  void gcp10_low_dens_eos(std::string model="BSk20",
761  bool external=false);
762 
763  /** \brief Crust EOS from \ref Newton13 given L in MeV
764 
765  Current acceptable values for \c model are <tt>PNM</tt>
766  and <tt>J35</tt>.
767  */
768  void ngl13_low_dens_eos(double L, std::string model="PNM",
769  bool external=false);
770 
771  /** \brief Crust EOS from \ref Newton13 given S and L in MeV
772  and a transition density
773 
774  Note that this function works only if \f$ 28 < S < 38 \f$ MeV,
775  \f$ 25 < L < 115 \f$ MeV, \f$ 0.01 < n_t < 0.15 \f$,
776  and \f$ L > 5 S-65~\mathrm{MeV} \f$
777  . If \c fname is a string of length 0 (the default),
778  then this function looks for a file named \c newton_SL.o2
779  in the \o2 data directory specified by
780  \ref o2scl::lib_settings_class::get_data_dir() .
781  */
782  void ngl13_low_dens_eos2(double S, double L, double nt,
783  std::string fname="");
784 
785  /// Compute with no crust EOS
787  use_crust=false;
788  return;
789  }
790  //@}
791 
792  /// \name Functions used by the tov_solve class
793  //@{
794  /** \brief Given the pressure, produce the energy and number densities
795 
796  The arguments \c pr and \c ed should always be in \f$
797  M_{\odot}/\mathrm{km}^3 \f$ . The argument for \c nb should be
798  in \f$ \mathrm{fm}^{-3} \f$ .
799 
800  If the baryon density is not specified, it should be set to
801  zero or \ref baryon_column should be set to false
802  */
803  virtual void ed_nb_from_pr(double pr, double &ed, double &nb);
804  //@}
805 
806  /// \name Other functions
807  //@{
808  /** \brief Get the energy density from the pressure in the
809  user-specified unit system
810  */
811  virtual void get_eden_user(double pres, double &ed, double &nb);
812 
813  /** \brief Get the transition pressure (in the user-specified
814  unit system) and width
815  */
816  void get_transition(double &ptrans, double &pwidth);
817 
818  /** \brief Set the transition pressure and "width"
819 
820  Sets the transition pressure and the width (specified as a
821  number greater than unity in \c pw) of the transition between
822  the two EOSs. The transition should be in the same units of
823  the user-specified EOS. The transition is done smoothly using
824  linear interpolation between \f$ P=\mathrm{ptrans}/pmathrm{pw}
825  \f$ and \f$ P=\mathrm{ptrans} \times pmathrm{pw} \f$.
826  */
827  void set_transition(double ptrans, double pw);
828  //@}
829 
830  /// \name User EOS
831  //@{
832  /// Energy densities from full EOS
833  std::vector<double> full_vece;
834  /// Pressures from full EOS
835  std::vector<double> full_vecp;
836  /// Baryon densities from full EOS
837  std::vector<double> full_vecnb;
838  //@}
839 
840 #ifndef DOXYGEN_INTERNAL
841 
842  protected:
843 
844  /** \brief Internal function to reinterpolate if if either the
845  core or crust tables are changed
846  */
847  void internal_read();
848 
849  /// \name Crust EOS
850  //@{
851  /// Set to true if we are using a crust EOS (default false)
852  bool use_crust;
853 
854  /// Energy densities
855  std::vector<double> crust_vece;
856  /// Pressures
857  std::vector<double> crust_vecp;
858  /// Baryon densities
859  std::vector<double> crust_vecnb;
860  //@}
861 
862  /// \name Core EOS
863  //@{
864  /// Energy densities
865  std::vector<double> core_vece;
866  /// Pressures
867  std::vector<double> core_vecp;
868  /// Baryon densities
869  std::vector<double> core_vecnb;
870  //@}
871 
872  /// \name Interpolation objects
873  //@{
876  interp<std::vector<double> > gen_int;
877  //@}
878 
879  /// \name Unit conversion factors for core EOS
880  //@{
881  /// Unit conversion factor for energy density (default 1.0)
882  double efactor;
883  /// Unit conversion factor for pressure (default 1.0)
884  double pfactor;
885  /// Unit conversion factor for baryon density (default 1.0)
886  double nfactor;
887  //@}
888 
889  /// \name Properties of transition
890  //@{
891  /** \brief Transition pressure (in \f$ M_{\odot}/\mathrm{km}^3 \f$)
892  */
893  double trans_pres;
894  /// Transition width (unitless)
895  double trans_width;
896  //@}
897 
898 #endif
899 
900  };
901 
902 #ifndef DOXYGEN_NO_O2NS
903 }
904 #endif
905 
906 #endif
907 
908 
std::vector< double > full_vecnb
Baryon densities from full EOS.
Definition: eos_tov.h:837
void read_vectors_swap(size_t user_n, vec_t &user_ed, vec_t &user_pr)
Read the EOS from a pair of equal length vectors for energy density and pressure. ...
Definition: eos_tov.h:489
std::vector< double > core_vecnb
Baryon densities.
Definition: eos_tov.h:869
double pr1
The pressure for which the baryon density is known.
Definition: eos_tov.h:388
Linear EOS .
Definition: eos_tov.h:374
vec_t ed
Energy densities from full EOS.
Definition: eos_tov.h:592
virtual double nb_from_ed(double ed)=0
From the energy density, return the baryon density.
std::vector< double > crust_vecnb
Baryon densities.
Definition: eos_tov.h:859
std::vector< double > crust_vecp
Pressures.
Definition: eos_tov.h:857
virtual double nb_from_pr(double pr)=0
From the pressure, return the baryon density.
virtual void ed_nb_from_pr(double pr, double &ed, double &nb)=0
Given the pressure, produce the energy and number densities.
virtual double ed_from_pr(double pr)
From the pressure, return the energy density.
Definition: eos_tov.h:536
std::vector< double > full_vecp
Pressures from full EOS.
Definition: eos_tov.h:835
virtual double ed_from_nb(double nb)=0
From the baryon density, return the energy density.
double eps0
The energy density at zero pressure (default 0.0)
Definition: eos_tov.h:395
virtual double ed_from_nb(double nb)
From the baryon density, return the energy density.
Definition: eos_tov.h:560
double ed1
The energy density for which the baryon density is known.
Definition: eos_tov.h:167
double pr1
The pressure at ed1.
Definition: eos_tov.h:301
double trans_pres
Transition pressure (in )
Definition: eos_tov.h:893
double K
Coefficient (default 1.0)
Definition: eos_tov.h:305
double ed1
The energy density for which the baryon density is known.
Definition: eos_tov.h:384
virtual double nb_from_ed(double ed)
From the energy density, return the baryon density.
Definition: eos_tov.h:548
double pfactor
Unit conversion factor for pressure (default 1.0)
Definition: eos_tov.h:884
std::vector< double > core_vecp
Pressures.
Definition: eos_tov.h:867
virtual double pr_from_ed(double ed)=0
From the energy density, return the pressure.
virtual void ed_nb_from_pr(double pr, double &ed, double &nb)
Given the pressure, produce the energy and number densities.
Definition: eos_tov.h:578
double Pstar
The parameter with units of pressure in units of solar masses per km cubed (default value ) ...
Definition: eos_tov.h:183
void read_vectors_copy(size_t user_n, vec_t &user_ed, vec_t &user_pr)
Read the EOS from a pair of equal length vectors for energy density and pressure. ...
Definition: eos_tov.h:522
void read_vectors_swap(size_t user_n, vec_t &user_ed, vec_t &user_pr, vec_t &user_nb)
Read the EOS from a set of equal length vectors for energy density, pressure, and baryon density...
Definition: eos_tov.h:473
std::vector< double > crust_vece
Energy densities.
Definition: eos_tov.h:855
void reset_interp_nb(size_t n)
Desc.
Definition: eos_tov.h:456
double nb1
The baryon density at ed1.
Definition: eos_tov.h:380
virtual double nb_from_pr(double pr)
From the pressure, return the baryon density.
Definition: eos_tov.h:554
void vector_copy(const vec_t &src, vec2_t &dest)
vec_t pr
Pressures from full EOS.
Definition: eos_tov.h:594
double n
Index (default 3.0)
Definition: eos_tov.h:308
double ed1
The energy density for which the baryon density is known.
Definition: eos_tov.h:297
An EOS for the TOV solver using simple linear interpolation and an optional crust EOS...
Definition: eos_tov.h:675
bool has_baryons()
Return true if a baryon density is available.
Definition: eos_tov.h:66
A EOS base class for the TOV solver.
Definition: eos_tov.h:47
virtual double pr_from_nb(double nb)
From the baryon density, return the pressure.
Definition: eos_tov.h:566
Provide an EOS for TOV solvers based on interpolation of user-supplied vectors.
Definition: eos_tov.h:444
virtual double ed_from_pr(double pr)=0
From the pressure, return the energy density.
void no_low_dens_eos()
Compute with no crust EOS.
Definition: eos_tov.h:786
Standard polytropic EOS .
Definition: eos_tov.h:287
double cs2
Coefficient (default 1.0)
Definition: eos_tov.h:392
double pr1
The pressure at ed1.
Definition: eos_tov.h:171
void check_nb(double &avg_abs_dev, double &max_abs_dev)
Check that the baryon density is consistent with the .
double nb1
The baryon density at ed1.
Definition: eos_tov.h:163
void read_vectors_copy(size_t user_n, vec_t &user_ed, vec_t &user_pr, vec_t &user_nb)
Read the EOS from a set of equal length vectors for energy density, pressure, and baryon density...
Definition: eos_tov.h:503
double efactor
Unit conversion factor for energy density (default 1.0)
Definition: eos_tov.h:882
double nfactor
Unit conversion factor for baryon density (default 1.0)
Definition: eos_tov.h:886
vec_t nb
Baryon densities from full EOS.
Definition: eos_tov.h:596
bool use_crust
Set to true if we are using a crust EOS (default false)
Definition: eos_tov.h:852
int verbose
Control for output (default 1)
Definition: eos_tov.h:63
bool baryon_column
Set to true if the baryon density is provided in the EOS (default false)
Definition: eos_tov.h:54
std::vector< double > full_vece
Energy densities from full EOS.
Definition: eos_tov.h:833
virtual double pr_from_ed(double ed)
From the energy density, return the pressure.
Definition: eos_tov.h:542
void reset_interp(size_t n)
Desc.
Definition: eos_tov.h:448
double nb1
The baryon density at ed1.
Definition: eos_tov.h:293
virtual double pr_from_nb(double nb)=0
From the baryon density, return the pressure.
double trans_width
Transition width (unitless)
Definition: eos_tov.h:895
std::vector< double > core_vece
Energy densities.
Definition: eos_tov.h:865
itp_linear
The Buchdahl EOS for the TOV solver.
Definition: eos_tov.h:157

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