eos_tov.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_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 Internal function to reset the interpolation
447  */
448  void reset_interp(size_t n) {
449  pe_int.set(n,pr_vec,ed_vec,itp_linear);
450  ep_int.set(n,ed_vec,pr_vec,itp_linear);
451  return;
452  }
453 
454  /** \brief Internal function to reset the interpolation
455  with baryon density
456  */
457  void reset_interp_nb(size_t n) {
458  reset_interp(n);
459  pn_int.set(n,pr_vec,nb_vec,itp_linear);
460  np_int.set(n,nb_vec,pr_vec,itp_linear);
461  en_int.set(n,ed_vec,nb_vec,itp_linear);
462  ne_int.set(n,nb_vec,ed_vec,itp_linear);
463  return;
464  }
465 
466  public:
467 
468  /** \brief Read the EOS from a set of equal length
469  vectors for energy density, pressure, and baryon density
470 
471  In this version, the user-specified vectors are swapped
472  with internal storage.
473  */
474  void read_vectors_swap(size_t user_n, vec_t &user_ed, vec_t &user_pr,
475  vec_t &user_nb) {
476  std::swap(user_ed,ed_vec);
477  std::swap(user_pr,pr_vec);
478  std::swap(user_nb,nb_vec);
479  this->baryon_column=true;
480  reset_interp_nb(user_n);
481  return;
482  }
483 
484  /** \brief Read the EOS from a pair of equal length
485  vectors for energy density and pressure
486 
487  In this version, the user-specified vectors are swapped
488  with internal storage.
489  */
490  void read_vectors_swap(size_t user_n, vec_t &user_ed, vec_t &user_pr) {
491  std::swap(user_ed,ed_vec);
492  std::swap(user_pr,pr_vec);
493  this->baryon_column=false;
494  reset_interp(user_n);
495  return;
496  }
497 
498  /** \brief Read the EOS from a set of equal length
499  vectors for energy density, pressure, and baryon density
500 
501  In this version, the user-specified vectors are copied
502  to internal storage.
503  */
504  void read_vectors_copy(size_t user_n, vec_t &user_ed, vec_t &user_pr,
505  vec_t &user_nb) {
506  if (ed_vec.size()!=user_n) ed_vec.resize(user_n);
507  if (pr_vec.size()!=user_n) pr_vec.resize(user_n);
508  if (nb_vec.size()!=user_n) nb_vec.resize(user_n);
509  vector_copy(user_ed,ed_vec);
510  vector_copy(user_pr,pr_vec);
511  vector_copy(user_nb,nb_vec);
512  this->baryon_column=true;
513  reset_interp_nb(user_n);
514  return;
515  }
516 
517  /** \brief Read the EOS from a pair of equal length
518  vectors for energy density and pressure
519 
520  In this version, the user-specified vectors are copied
521  to internal storage.
522  */
523  void read_vectors_copy(size_t user_n, vec_t &user_ed, vec_t &user_pr) {
524  if (ed_vec.size()!=user_n) ed_vec.resize(user_n);
525  if (pr_vec.size()!=user_n) pr_vec.resize(user_n);
526  vector_copy(user_ed,ed_vec);
527  vector_copy(user_pr,pr_vec);
528  this->baryon_column=false;
529  reset_interp(user_n);
530  return;
531  }
532 
533  /// \name Basic EOS functions
534  //@{
535  /** \brief From the pressure, return the energy density
536  */
537  virtual double ed_from_pr(double pr) {
538  return pe_int.eval(pr);
539  }
540 
541  /** \brief From the energy density, return the pressure
542  */
543  virtual double pr_from_ed(double ed) {
544  return ep_int.eval(ed);
545  }
546 
547  /** \brief From the energy density, return the baryon density
548  */
549  virtual double nb_from_ed(double ed) {
550  return en_int.eval(ed);
551  }
552 
553  /** \brief From the pressure, return the baryon density
554  */
555  virtual double nb_from_pr(double pr) {
556  return pn_int.eval(pr);
557  }
558 
559  /** \brief From the baryon density, return the energy density
560  */
561  virtual double ed_from_nb(double nb) {
562  return ne_int.eval(nb);
563  }
564 
565  /** \brief From the baryon density, return the pressure
566  */
567  virtual double pr_from_nb(double nb) {
568  return np_int.eval(nb);
569  }
570 
571  /** \brief Given the pressure, produce the energy and number densities
572 
573  The arguments \c pr and \c ed should always be in \f$
574  M_{\odot}/\mathrm{km}^3 \f$ . The argument for \c nb should be
575  in \f$ \mathrm{fm}^{-3} \f$ .
576 
577  If \ref baryon_column is false, then \c nb is unmodified.
578  */
579  virtual void ed_nb_from_pr(double pr, double &ed, double &nb) {
580  ed=ed_from_pr(pr);
581  if (this->baryon_column) {
582  nb=nb_from_pr(pr);
583  }
584  return;
585  }
586  //@}
587 
588  protected:
589 
590  /// \name EOS storage
591  //@{
592  /// Energy densities from full EOS
593  vec_t ed_vec;
594  /// Pressures from full EOS
595  vec_t pr_vec;
596  /// Baryon densities from full EOS
597  vec_t nb_vec;
598  //@}
599 
600  /// \name Interpolators
601  //@{
602  interp_vec<vec_t> pe_int;
603  interp_vec<vec_t> pn_int;
604  interp_vec<vec_t> ep_int;
605  interp_vec<vec_t> en_int;
606  interp_vec<vec_t> np_int;
607  interp_vec<vec_t> ne_int;
608  //@}
609 
610  };
611 
612  /** \brief An EOS for the TOV solver using simple linear
613  interpolation and an optional crust EOS
614 
615  The simplest usage of this class is simply to use \ref
616  read_table() to read a tabulated EOS stored in a \ref
617  table_units object and optionally specify a separate crust EOS.
618 
619  Alternatively, the user can simply specify objects
620  of type <tt>std::vector<double></tt> which store the energy
621  density, pressure, and baryon density (which should
622  include the crust if necessary).
623 
624  There are two methods to handle the crust-core interface. The
625  method labeled <tt>smooth_trans</tt> uses the crust below
626  pressure \f$ P_{\mathrm{lo}} \f$ (equal to the value of \ref
627  trans_pres divided by \ref trans_width) and the core above
628  pressure \f$ P_{\mathrm{hi}} \f$ (the value of \ref trans_pres
629  times \ref trans_width) and then in between uses
630  \f[
631  \varepsilon(P) = [1-\chi(P)] \varepsilon_{\mathrm{crust}}(P) +
632  \chi(P) \varepsilon_{\mathrm{core}}(P)
633  \f]
634  where the value \f$ \chi(P) \f$ is determined by
635  \f[
636  \chi(P) = (P-P_{\mathrm{lo}})/
637  (P_{\mathrm{hi}}-P_{\mathrm{lo}}) \, .
638  \f]
639  This method is a bit more faithful to the original EOS tables,
640  but the matching can result in pressures which decrease with
641  increasing energy density. Alternatively the <tt>match_line</tt>
642  method uses \f$
643  \varepsilon_{\mathrm{lo}}=\varepsilon_{\mathrm{crust}}
644  (P_{\mathrm{lo}}) \f$ and \f$
645  \varepsilon_{\mathrm{hi}}=\varepsilon_{\mathrm{core}}
646  (P_{\mathrm{hi}}) \f$ and
647  \f[
648  \varepsilon(P) = (\varepsilon_{\mathrm{hi}} -
649  \varepsilon_{\mathrm{lo}}) \chi
650  + \varepsilon_{\mathrm{lo}} \, .
651  \f]
652  (using the same expression for \f$ \chi \f$ ). This method less
653  frequently results in decreasing pressures, but can deviate
654  further from the original tables. For the baryon density,
655  expression similar to those used for \f$ \varepsilon(P) \f$
656  are used for \f$ n_B(P) \f$ .
657 
658  By default, no crust EOS is used. If a crust EOS is specified
659  through one of the crust EOS functions, then by default \ref
660  trans_width is 1.0 and <tt>transition_mode</tt> is set equal
661  <tt>smooth_trans</tt>. This creates a discontinuous energy
662  density between the core and crust EOS at the transition
663  pressure. A smoother transition can be chosen by increasing \ref
664  trans_width to a value larger than 1. The crust EOS can be
665  changed after the core EOS is specified.
666 
667  The value of \ref trans_pres is set either by \ref
668  set_transition(), or any of the crust EOS functions.
669 
670  Internally, energy and pressure are stored in units of \f$
671  \mathrm{M}_{\odot}/\mathrm{km}^3 \f$ and baryon density is
672  stored in units of \f$ \mathrm{fm}^{-3} \f$ . The user-specified
673  EOS table is left as is, and unit conversion is performed as
674  needed in ed_nb_from_pr() and other functions from the units
675  specified in the input \ref table_units object.
676 
677  \comment
678  \todo It might be useful to exit more gracefully when non-finite
679  values are obtained in interpolation, analogous to the
680  err_nonconv mechanism elsewhere.
681  2/6/16: I'm not sure this is really necessary. It's only
682  linear interpolation, so the err_nonconv mechanism probably
683  won't be useful.
684  \endcomment
685 
686  \future Create a sanity check where core_auxp is nonzero only if
687  core_table is also nonzero. Alternatively, this complication is
688  due to the fact that this class works in two ways: one where it
689  reads a table (and adds a crust), and one where it reads in
690  vectors (with no crust). Maybe these two modes of operation
691  should be separated into two classes. (2/6/16: Now there is
692  a new eos_tov_vector class. The best way forward may be
693  to make this a child of eos_tov_vectors.)
694  */
695  class eos_tov_interp : public eos_tov {
696 
697  public:
698 
699  eos_tov_interp();
700 
701  virtual ~eos_tov_interp();
702 
703  /// \name Mode of transitioning between crust and core EOS
704  //@{
705  int transition_mode;
706  static const int smooth_trans=0;
707  static const int match_line=1;
708  //@}
709 
710  /** \brief If true, call the error handler if the EOS
711  reports a non-finite value (default true)
712  */
714 
715  /// \name Basic EOS functions (all in the internal unit system)
716  //@{
717  /** \brief From the pressure, return the energy density
718  */
719  virtual double ed_from_pr(double pr);
720 
721  /** \brief From the energy density, return the pressure
722  */
723  virtual double pr_from_ed(double ed);
724 
725  /** \brief From the energy density, return the baryon density
726  */
727  virtual double nb_from_ed(double ed);
728 
729  /** \brief From the pressure, return the baryon density
730  */
731  virtual double nb_from_pr(double pr);
732 
733  /** \brief From the baryon density, return the energy density
734  */
735  virtual double ed_from_nb(double nb);
736 
737  /** \brief From the baryon density, return the pressure
738  */
739  virtual double pr_from_nb(double nb);
740  //@}
741 
742  /// \name Basic usage
743  //@{
744  /** \brief Specify the EOS through a table
745 
746  If units are specified for any of the columns, then this
747  function attempts to automatically determine the correct
748  conversion factors using the \ref o2scl::convert_units object
749  returned by \ref o2scl::o2scl_settings . If the units for any
750  of the columns are blank, then they are assumed to be the
751  native units for \ref o2scl::tov_solve .
752 
753  \note The input table must have at least 2 rows and
754  the pressure column must be strictly increasing.
755 
756  This function copies the needed information from the
757  table so if it is modified then this function
758  needs to be called again to read a new table.
759  */
760  void read_table(table_units<> &eosat, std::string s_cole,
761  std::string s_colp, std::string s_colnb="");
762  //@}
763 
764  /// \name Crust EOS functions
765  //@{
766  /// Standard crust EOS from \ref Negele73 and \ref Baym71tg
767  void default_low_dens_eos();
768 
769  /// Crust EOS from \ref Shen11b
770  void sho11_low_dens_eos();
771 
772  /** \brief Crust EOS from \ref Steiner12
773 
774  This function uses the neutron star crust models from \ref
775  Steiner12 . The current acceptable values for \c model are
776  <tt>APR</tt>, <tt>Gs</tt>, <tt>Rs</tt> and <tt>SLy4</tt>.
777  */
778  void s12_low_dens_eos(std::string model="SLy4",
779  bool external=false);
780 
781  /** \brief Crust EOS from Goriely, Chamel, and Pearson
782 
783  From \ref Goriely10, \ref Pearson11, and \ref Pearson12 .
784  */
785  void gcp10_low_dens_eos(std::string model="BSk20",
786  bool external=false);
787 
788  /** \brief Crust EOS from \ref Newton13 given L in MeV
789 
790  Current acceptable values for \c model are <tt>PNM</tt>
791  and <tt>J35</tt>.
792  */
793  void ngl13_low_dens_eos(double L, std::string model="PNM",
794  bool external=false);
795 
796  /** \brief Crust EOS from \ref Newton13 given S and L in MeV
797  and a transition density
798 
799  Note that this function works only if \f$ 28 < S < 38 \f$ MeV,
800  \f$ 25 < L < 115 \f$ MeV, \f$ 0.01 < n_t < 0.15 \f$,
801  and \f$ L > 5 S-65~\mathrm{MeV} \f$
802  . If \c fname is a string of length 0 (the default),
803  then this function looks for a file named \c newton_SL.o2
804  in the \o2 data directory specified by
805  \ref o2scl::lib_settings_class::get_data_dir() .
806  */
807  void ngl13_low_dens_eos2(double S, double L, double nt,
808  std::string fname="");
809 
810  /// Compute with no crust EOS (this is the default)
812  use_crust=false;
813  return;
814  }
815  //@}
816 
817  /// \name Functions used by the tov_solve class
818  //@{
819  /** \brief Given the pressure, produce the energy and number densities
820 
821  The arguments \c pr and \c ed should always be in \f$
822  M_{\odot}/\mathrm{km}^3 \f$ . The argument for \c nb should be
823  in \f$ \mathrm{fm}^{-3} \f$ .
824 
825  If the baryon density is not specified, it should be set to
826  zero or \ref baryon_column should be set to false
827  */
828  virtual void ed_nb_from_pr(double pr, double &ed, double &nb);
829  //@}
830 
831  /// \name Other functions
832  //@{
833  /** \brief Get the energy density from the pressure in the
834  user-specified unit system
835  */
836  virtual void get_eden_user(double pres, double &ed, double &nb);
837 
838  /** \brief Get the transition pressure (in the user-specified
839  unit system) and width
840  */
841  void get_transition(double &ptrans, double &pwidth);
842 
843  /** \brief Set the transition pressure and "width"
844 
845  Sets the transition pressure and the width (specified as a
846  number greater than unity in \c pw) of the transition between
847  the two EOSs. The transition should be in the same units of
848  the user-specified EOS. The transition is done smoothly using
849  linear interpolation between \f$ P=\mathrm{ptrans}/pmathrm{pw}
850  \f$ and \f$ P=\mathrm{ptrans} \times pmathrm{pw} \f$.
851  */
852  void set_transition(double ptrans, double pw);
853  //@}
854 
855  /// \name Full EOS arrays (in internal units)
856  //@{
857  /// Energy density
858  std::vector<double> full_vece;
859  /// Pressure
860  std::vector<double> full_vecp;
861  /// Baryon density
862  std::vector<double> full_vecnb;
863  //@}
864 
865 #ifndef DOXYGEN_INTERNAL
866 
867  protected:
868 
869  /** \brief Internal function to reinterpolate if if either the
870  core or crust tables are changed
871  */
872  void internal_read();
873 
874  /// \name Crust EOS
875  //@{
876  /// Set to true if we are using a crust EOS (default false)
877  bool use_crust;
878 
879  /// Energy densities
880  std::vector<double> crust_vece;
881  /// Pressures
882  std::vector<double> crust_vecp;
883  /// Baryon densities
884  std::vector<double> crust_vecnb;
885  //@}
886 
887  /// \name Core EOS
888  //@{
889  /// Energy densities
890  std::vector<double> core_vece;
891  /// Pressures
892  std::vector<double> core_vecp;
893  /// Baryon densities
894  std::vector<double> core_vecnb;
895  //@}
896 
897  /// \name Interpolation objects
898  //@{
901  interp<std::vector<double> > gen_int;
902  //@}
903 
904  /// \name Unit conversion factors for core EOS
905  //@{
906  /// Unit conversion factor for energy density (default 1.0)
907  double efactor;
908  /// Unit conversion factor for pressure (default 1.0)
909  double pfactor;
910  /// Unit conversion factor for baryon density (default 1.0)
911  double nfactor;
912  //@}
913 
914  /// \name Properties of transition
915  //@{
916  /** \brief Transition pressure (in \f$ M_{\odot}/\mathrm{km}^3 \f$)
917  */
918  double trans_pres;
919  /// Transition width (unitless)
920  double trans_width;
921  //@}
922 
923 #endif
924 
925  };
926 
927 #ifndef DOXYGEN_NO_O2NS
928 }
929 #endif
930 
931 #endif
932 
933 
bool err_nonconv
If true, call the error handler if the EOS reports a non-finite value (default true) ...
Definition: eos_tov.h:713
std::vector< double > full_vecnb
Baryon density.
Definition: eos_tov.h:862
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:490
std::vector< double > core_vecnb
Baryon densities.
Definition: eos_tov.h:894
double pr1
The pressure for which the baryon density is known.
Definition: eos_tov.h:388
Linear EOS .
Definition: eos_tov.h:374
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:884
std::vector< double > crust_vecp
Pressures.
Definition: eos_tov.h:882
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:537
std::vector< double > full_vecp
Pressure.
Definition: eos_tov.h:860
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:561
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:918
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
vec_t pr_vec
Pressures from full EOS.
Definition: eos_tov.h:595
virtual double nb_from_ed(double ed)
From the energy density, return the baryon density.
Definition: eos_tov.h:549
double pfactor
Unit conversion factor for pressure (default 1.0)
Definition: eos_tov.h:909
std::vector< double > core_vecp
Pressures.
Definition: eos_tov.h:892
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:579
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:523
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:474
std::vector< double > crust_vece
Energy densities.
Definition: eos_tov.h:880
void reset_interp_nb(size_t n)
Internal function to reset the interpolation with baryon density.
Definition: eos_tov.h:457
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:555
void vector_copy(const vec_t &src, vec2_t &dest)
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:695
bool has_baryons()
Return true if a baryon density is available.
Definition: eos_tov.h:66
vec_t ed_vec
Energy densities from full EOS.
Definition: eos_tov.h:593
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:567
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 (this is the default)
Definition: eos_tov.h:811
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:504
double efactor
Unit conversion factor for energy density (default 1.0)
Definition: eos_tov.h:907
double nfactor
Unit conversion factor for baryon density (default 1.0)
Definition: eos_tov.h:911
bool use_crust
Set to true if we are using a crust EOS (default false)
Definition: eos_tov.h:877
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 density.
Definition: eos_tov.h:858
virtual double pr_from_ed(double ed)
From the energy density, return the pressure.
Definition: eos_tov.h:543
void reset_interp(size_t n)
Internal function to reset the interpolation.
Definition: eos_tov.h:448
vec_t nb_vec
Baryon densities from full EOS.
Definition: eos_tov.h:597
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:920
std::vector< double > core_vece
Energy densities.
Definition: eos_tov.h:890
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).