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

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