tov_solve.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 tov_solve.h
24  \brief File defining \ref o2scl::tov_solve
25 */
26 #ifndef O2SCL_TOV_SOLVE_H
27 #define O2SCL_TOV_SOLVE_H
28 
29 #include <o2scl/eos_tov.h>
30 #include <o2scl/interp.h>
31 #include <o2scl/table_units.h>
32 #include <o2scl/ode_iv_solve.h>
33 #include <o2scl/mroot_hybrids.h>
34 #include <o2scl/min_brent_gsl.h>
35 
36 #ifndef DOXYGEN_NO_O2NS
37 namespace o2scl {
38 #endif
39 
40  /** \brief Class to solve the Tolman-Oppenheimer-Volkov equations
41 
42  See the \ref tovtoc section of the User's Guide for
43  the mathematical background.
44 
45  This class uses adaptive integration to compute the
46  gravitational mass, the radius, the baryonic mass (if the EOS
47  supplies the baryon density), and the gravitational potential
48  (if requested). The equation of state may be changed at any
49  time, by specifying the appropriate \ref eos_tov object
50 
51  <b>Basic Usage</b>
52 
53  After specifying the EOS through \ref tov_solve::set_eos(), one
54  calls either \ref tov_solve::mvsr(), \ref tov_solve::fixed(),
55  \ref tov_solve::max() or \ref tov_solve::fixed_pr() to compute
56  the mass versus radius curve, the profile of a star of a given
57  fixed mass, the profile of the maximum mass star, or the profile
58  of a star with a fixed central pressure. These functions all
59  generate results in the form of a \ref table_units object, which
60  can be obtained from \ref tov_solve::get_results().
61 
62  Screen output:
63  - \c verbose=0 - Nothing
64  - \c verbose=1 - Basic information
65  - \c verbose=2 - For each profile computation, report solution
66  information at every kilometer.
67  - \c verbose=3 - Report profile information at every 20 grid
68  points and output the final interpolation to zero pressure. A
69  keypress is required after each profile.
70 
71  <b>Mass versus radius curve</b>
72 
73  The neutron star mass versus radius curve is constructed by
74  computing profiles of all neutron stars with a range of central
75  pressures. This range is from \ref prbegin to \ref prend, and
76  the ratio between successive central pressures is specified in
77  \ref princ (making a logarithmic central pressure grid).
78 
79  \note If \ref prend is too close to (or smaller than) the
80  central pressure of the maximum mass star, then the mass-radius
81  curve generated by mvsr() might not include the maximum mass
82  star.
83 
84  \note The table generated by mvsr() may include
85  unstable configurations (including those with central
86  pressures larger than that of the maximum mass star).
87 
88  <b>Profiles of fixed mass</b>
89 
90  Profiles for a fixed gravitational mass are computed by solving
91  for the correct central pressure. In order to ensure that the
92  solver does not accidentally select a central pressure beyond
93  the maximum mass neutron star, the profile with the maximum mass
94  is computed beforehand automatically (using \ref max()). The
95  intial guess to the solver is always the value of \ref
96  fixed_pr_guess, which defaults to \f$ 5.2 \times
97  10^{-5}~\mathrm{Msun}/\mathrm{km}^3 \f$ . Alternatively, the
98  user can specify the central pressure of the maximum mass star
99  so that it does not have to be computed.
100 
101  In order to handle multiply-branched mass-radius relations, the
102  value of \ref fixed_pr_guess must be specified in order to
103  ensure that the correct profile is generated. Even if \ref
104  fixed_pr_guess is specified, there is no guarantee that the
105  Newton-Raphson will select the correct profile (though it is
106  very likely if the guess for the central pressure is
107  sufficiently accurate).
108 
109  The \ref fixed() function does not support computing profiles
110  with central pressures with gravitational masses larger than the
111  maximum value. For this, use \ref fixed_pr() .
112 
113  <b>Profile for the maximum mass star</b>
114 
115  This is provided by the function \ref max() . Internally, this
116  uses the minimizer specified by \ref set_min() or the default
117  minimizer, \ref def_min, to minimize the function \ref max_fun()
118  . In order to generate a good initial guess, the \ref max()
119  function begins by looping over central pressures from \ref
120  max_begin to \ref max_end and choosing the best guess from that
121  set of configurations.
122 
123  <b>Profile for a fixed central pressure</b>
124 
125  This is provided by the function \ref fixed_pr() and
126  is relatively fast because it does not require any solving
127  or minimizing. This function allows central pressures
128  larger than that of the maximum mass star.
129 
130  <b>Output tables</b>
131 
132  The functions \ref tov_solve::fixed(), \ref tov_solve::fixed_pr(),
133  and \ref tov_solve::max()
134  produce output tables which represent the profile of the neutron
135  star of the requested mass. The columns are
136  - \c gm, the enclosed gravitational mass in \f$ \mathrm{M}_{\odot} \f$
137  - \c r, the radial coordinate in \f$ \mathrm{km} \f$
138  - \c gp, the gravitational potential (unitless) when
139  \ref calc_gpot is true
140  - \c rjw, the value of \f$ g \f$ (see definition below; present
141  if \ref ang_vel is true)
142  - \c omega_rat, the value of \f$ f \f$ (see definition below;
143  present if \ref ang_vel is true)
144  - \c bm, the baryonic mass in \f$ \mathrm{M}_{\odot} \f$ (when
145  \ref eos_tov::baryon_column is true).
146  - \c pr, the pressure in user-specified units
147  - \c ed, the energy density in user-specified units
148  - \c nb, the baryon density in user-specified units
149  (if \ref eos_tov::baryon_column is true)
150  - \c sg, the local surface gravity
151  (in \f$ \mathrm{g}/\mathrm{cm}^{2} \f$ )
152  - \c rs, the local redshift (unitless),
153  - \c dmdr, the derivative of the enclosed gravitational mass
154  in \f$ \mathrm{M}_{\odot}/\mathrm{km} \f$
155  - \c dlogpdr, the derivative of the natural logarithm of the
156  pressure
157  - \c dgpdr, the derivative of the gravitational potential
158  in \f$ 1/\mathrm{km} \f$ (if \ref calc_gpot is true)
159  - \c dbmdr, the derivative of the enclosed baryonic mass
160  (if \ref eos_tov::baryon_column is true).
161 
162  The function \ref tov_solve::mvsr() produces a different kind of
163  output table corresponding to the mass versus radius curve. Some
164  points on the curve may correspond to unstable configurations.
165 
166  - \c gm, the total gravitational mass in \f$ \mathrm{M}_{\odot} \f$
167  - \c r, the radius in \f$ \mathrm{km} \f$
168  - \c gp, the gravitational potential in the center (unitless) when
169  \ref calc_gpot is true
170  - \c rjw, the value of \f$ g \f$ at the surface
171  (see definition below; present if \ref ang_vel is true)
172  - \c omega_rat, the value of \f$ f \f$ at the surface
173  (see definition below; present if \ref ang_vel is true)
174  - \c bm, total the baryonic mass in \f$ \mathrm{M}_{\odot} \f$ (when
175  \ref eos_tov::baryon_column is true).
176  - \c pr, the central pressure in user-specified units
177  - \c ed, the central energy density in user-specified units
178  - \c nb, the central baryon density in user-specified units
179  (if \ref eos_tov::baryon_column is true)
180  - \c sg, the surface gravity
181  (in \f$ \mathrm{g}/\mathrm{cm}^{2} \f$ )
182  - \c rs, the redshift at the surface,
183  - \c dmdr, the derivative of the gravitational mass
184  - \c dlogpdr, the derivative of the natural logarithm of the
185  pressure
186  - \c dgpdr, the derivative of the gravitational potential
187  in \f$ 1/\mathrm{km} \f$ (if \ref calc_gpot is true)
188  - \c dbmdr, the derivative of the enclosed baryonic mass
189  (if \ref eos_tov::baryon_column is true).
190 
191  <b>Unit systems</b>
192 
193  By default, this class operates with energy density and
194  pressure in units of \f$ \mathrm{M}_{\odot}/\mathrm{km}^3 \f$
195  and baryon density in units of \f$ \mathrm{fm}^{-3} \f$.
196 
197  The function \ref set_units(std::string,std::string,std::string)
198  allows one to use different unit systems for energy density,
199  pressure, and baryon density. The following list of units of
200  energy density and pressure are hard-coded into the library and
201  always work:
202  - <tt>"g/cm^3"</tt>,
203  - <tt>"erg/cm^3"</tt>,
204  - <tt>"dyne/cm^2"</tt>,
205  - <tt>"MeV/fm^3"</tt>,
206  - <tt>"1/fm^4"</tt>, and
207  - <tt>"Msun/km^3"</tt> (i.e. \f$ \mathrm{M}_{\odot}/\mathrm{km}^3 \f$ )
208 
209  The list of hard-coded units for baryon density are:
210  - <tt>"1/m^3"</tt>,
211  - <tt>"1/cm^3"</tt>, and
212  - <tt>"1/fm^3"</tt>.
213 
214  Other units choices will work if the conversion is either
215  already added to the global unit conversion object (from
216  <tt>o2scl_settings.get_convert_units()</tt> ) or the global unit
217  conversion object is able to compute them by a <tt>system()</tt>
218  call to GNU <tt>units</tt> (see documentation in \ref
219  convert_units). Note that the choice of what units the tables
220  are produced in is independent of the unit system specified in
221  the associated \ref eos_tov object, i.e. the input EOS and
222  output EOS units need not be the same.
223 
224  Alternatively, using \ref set_units(double,double,double)
225  allows one to specify the conversion factors directly without
226  using the global unit conversion object.
227 
228  <b>Accuracy</b>
229 
230  The present code, as demonstrated in the tests, gives the
231  correct central pressure and energy density of the analytical
232  solution by Buchdahl to within less than 1 \part in \f$ 10^8 \f$.
233 
234  <b>Rotation</b>
235 
236  Rotation is considered if \ref tov_solve::ang_vel is set to
237  <tt>true</tt>. This adds two more differential equations to
238  solve for each central pressure.
239 
240  The differential equation for \f$ \bar{\omega} \f$ (see the
241  section in the User's Guide called \ref tovtoc ) is independent
242  of the relative scale for \f$ \bar{\omega} \f$ and \f$ j \f$ .
243  (Note that \f$ j \f$ is a metric function not simply related to
244  the angular momentum, \f$ J \f$ .) First, one rescales \f$
245  \bar{\omega} \f$ and rewrites everything in terms of \f$ f\equiv
246  \bar{\omega}/\Omega \f$ and \f$ g \equiv r^4 j~df/dr \f$ . The
247  quantity \f$ f \f$ is unitless and \f$ g \f$ has units of \f$
248  \mathrm{km}^3 \f$ . Second, pick a central pressure, \f$
249  m(r=0)=g(r=0)=0 \f$, arbitrary values for \f$ \Phi(r=0) \f$ and
250  \f$ f(r=0) \f$, and integrate
251  \f{eqnarray*}
252  \frac{d P}{d r} &=& - \frac{G \varepsilon m}{r^2}
253  \left( 1+\frac{P}{\varepsilon}\right)
254  \left( 1+\frac{4 \pi P r^3}{m} \right)
255  \left( 1-\frac{2 G m}{r}\right)^{-1}
256  \nonumber \\
257  \frac{d m}{d r} &=& 4 \pi r^2 \varepsilon
258  \nonumber \\
259  \frac{d \Phi}{d r} &=& - \frac{1}{\varepsilon}
260  \frac{ d P}{d r} \left(1+\frac{P}{\varepsilon}\right)^{-1}
261  \nonumber \\
262  \frac{d g}{dr} &=& -4 r^3 \frac{d j}{dr} f
263  \nonumber \\
264  \frac{d f}{dr} &=& \frac{g}{r^4 j}
265  \f}
266  Afterwards, shift \f$ \Phi \f$ by a constant to ensure
267  the correct value at \f$ r=R \f$, and multiply \f$ g \f$
268  by a constant to ensure that \f$ g=r^4 j (df/dr) \f$ holds
269  for the new potential \f$ \Phi \f$. Then, multiply
270  \f$ f \f$ by a constant to ensure that
271  \f[
272  f(r=R) + \frac{R}{3} \left(\frac{df}{dr}\right)_{r=R} = 1
273  \f]
274 
275  The functions \f$ f \f$ and \f$ g \f$ are stored in columns
276  called <tt>"omega_rat"</tt> and <tt>"rjw"</tt>, respectively.
277  One can compute the baryonic mass by integration or by adding
278  one additional differential equation, bringing the total to six.
279 
280  The moment of inertia is (see \ref tovtoc),
281  \f[
282  I = \frac{R^4}{6 G} \left.\frac{df}{dr}\right|_{r=R}
283  \f]
284  where the last fraction is stored in \ref domega_rat .
285  For an object named \c ts of type \ref tov_solve
286  after a call to \ref fixed(), \ref max(), or
287  \ref fixed_pr(), the moment of inertia can be computed
288  with, e.g.
289  \code
290  tov_solve ts;
291  ts.max();
292  double schwarz_km=o2scl_cgs::schwarzchild_radius/1.0e5;
293  double I=ts.domega_rat*pow(ts.rad,4.0)/3.0/schwarz_km;
294  \endcode
295 
296  After a call to \ref mvsr(), the values of \f$ f(r=R) \f$
297  and \f$ g(r=R) \f$ are stored in columns labeled
298  <tt>"omega_rat"</tt> and <tt>"rjw"</tt>. The moment
299  of inertia of a 1.4 solar mass neutron star
300  can be computed with, e.g.
301  \code
302  tov_solve ts;
303  ts.mvsr();
304  std::shared_ptr<table_units<> > tab=ts.get_results();
305  double schwarz_km=o2scl_cgs::schwarzchild_radius/1.0e5;
306  double I_14=tab->interp("gm",1.4,"rjw")/3.0/schwarz_km;
307  \endcode
308 
309  <b>Convergence details</b>
310 
311  By default, if the TOV solver fails to converge, the error
312  handler is called and an exception is thrown. If \ref
313  o2scl::tov_solve::err_nonconv is false, then \ref
314  o2scl::tov_solve::mvsr(), \ref o2scl::tov_solve::fixed(), and
315  \ref o2scl::tov_solve::max(), return an integer which gives some
316  information about why the solver failed to converge.
317 
318  If \ref err_nonconv is set to false, then the \ref fixed()
319  function temporarily sets the value of both \ref
320  o2scl::mroot::err_nonconv for the current solver and \ref
321  o2scl::jacobian::err_nonconv for the jacobian object, \ref
322  o2scl::mroot_hybrids::def_jac, of \ref def_solver equal to false.
323 
324  <b>Other details</b>
325 
326  The ODE solution is stored in a buffer which can be directly
327  accessed using \ref o2scl::tov_solve::get_rkx(), \ref
328  o2scl::tov_solve::get_rky(), and \ref
329  o2scl::tov_solve::get_rkdydx(). In the case that the ODE steps
330  are small enough that there isn't enough space in the buffers,
331  the ODE solution data is thinned by a factor of two to allow for
332  the remaining ODE steps to be stored. The size of the buffers
333  can be controlled in \ref buffer_size .
334 
335  If \ref o2scl::tov_solve::reformat_results is true (the
336  default), then the results are placed in a shared pointer to the
337  \ref table_units object which can be accessed using \ref
338  o2scl::tov_solve::get_results(). The maximum size of the output
339  table can be controlled with \ref max_table_size. The output
340  table may be smaller than this, as it cannot be larger than the
341  number of steps stored in the buffer.
342 
343  \note The function \ref o2scl::tov_solve::integ_star() returns
344  <tt>gsl_efailed</tt> without calling the error handler in the
345  case that the solver can recover gracefully from, for example, a
346  negative pressure.
347 
348  \future Convert to \ref o2scl::ode_iv_solve?
349  */
350  class tov_solve {
351 
352  public:
353 
356  typedef boost::numeric::ublas::matrix_column<ubmatrix> ubmatrix_column;
357  typedef boost::numeric::ublas::matrix_row<ubmatrix> ubmatrix_row;
358 
359 #ifndef DOXYGEN_INTERNAL
360 
361  protected:
362 
363  /// ODE function object
365 
366  /// Interpolation object for listed radii in \ref mvsr()
368 
369  /** \brief Set up column names and units
370 
371  When this function is used by \ref mvsr(), \c mvsr_mode is set
372  to true.
373  */
374  void column_setup(bool mvsr_mode=false);
375 
376  /** \brief If true, \ref integ_star() computes all the profile info,
377  otherwise it only computes the gravitational mass
378  */
380 
381  /// The last row index in \ref rky
382  size_t ix_last;
383 
384  /// The schwarzchild radius in km
385  double schwarz_km;
386 
387  /** \brief Target mass for integ_star()
388 
389  Has a value of zero, unless set to a non-zero value
390  by \ref fixed().
391  */
392  double tmass;
393 
394  /** \brief Ensure \c col does not match strings in \c cnames
395 
396  Underscores are added to \c col until it matches none of
397  the strings in \c cnames.
398  */
399  void make_unique_name(std::string &col,
400  std::vector<std::string> &cnames);
401 
402  /// \name User EOS
403  //@{
404  /// The EOS
406 
407  /// True if the EOS has been set
408  bool eos_set;
409  //@}
410 
411  /// \name Units for output table
412  //@{
413  /// Units for energy density
414  std::string eunits;
415  /// Units for pressure
416  std::string punits;
417  /// Units for baryon density
418  std::string nunits;
419  /// unit conversion factor for energy density (default 1.0)
420  double efactor;
421  /// unit conversion factor for pressure (default 1.0)
422  double pfactor;
423  /// unit conversion factor for baryon density (default 1.0)
424  double nfactor;
425  //@}
426 
427  /** \brief Smallest allowed pressure for integration (default: -100)
428 
429  This quantity can't be much smaller than -100 since we need to
430  compute numbers near \f$ e^{-\mathrm{min\_log\_pres}} \f$
431 
432  \future Replace this with the proper value from std::limits?
433  */
434  double min_log_pres;
435 
436  /// \name Integration storage
437  //@{
438  /// Radial coordinate (in kilometers)
439  ubvector rkx;
440  /** \brief ODE functions
441 
442  If \c rky is viewed as a matrix, then the first column of each
443  row is the gravitational mass in solar masses, and the second
444  column is the natural logarithm of the pressure in \f$
445  \mathrm{M}_{\odot}/km^3 \f$ . When \ref calc_gpot is true, the
446  next column is the gravitational potential (which is
447  unitless), and when \ref eos_tov::baryon_column is true, the
448  next column is the baryonic mass in \f$ \mathrm{M}_{\odot}
449  \f$.
450  */
451  std::vector<ubvector> rky;
452  /// The derivatives of the ODE functions
453  std::vector<ubvector> rkdydx;
454  //@}
455 
456  /// The output table
457  std::shared_ptr<table_units<> > out_table;
458 
459  /// \name Numerical methods
460  //@{
461  /// The solver for fixed gravitational masses
463 
464  /// The minimizer for maximum masses
466 
467  /// The adaptive stepper
469  //@}
470 
471  /// The ODE step function
472  virtual int derivs(double x, size_t nv, const ubvector &y,
473  ubvector &dydx);
474 
475  /// The minimizer function to compute the maximum mass
476  virtual double max_fun(double maxx);
477 
478  /** \brief The solver function to compute the stellar profile
479  */
480  virtual int integ_star(size_t ndvar, const ubvector &ndx,
481  ubvector &ndy);
482 
483 #endif
484 
485  public:
486 
487  tov_solve();
488 
489  virtual ~tov_solve();
490 
491  /// Size of the ODE solution buffer (default \f$ 10^5 \f$)
492  size_t buffer_size;
493 
494  /// Maximum number of lines in the output table (default 400)
496 
497  /// \name Basic properties
498  //@{
499  /// Gravitational mass (in \f$ \mathrm{M}_{\odot} \f$)
500  double mass;
501  /// Radius (in km)
502  double rad;
503  /// Baryonic mass (when computed; in \f$ \mathrm{M}_{\odot} \f$)
504  double bmass;
505  /// Gravitational potential (when computed; unitless)
506  double gpot;
507  /// The value of \f$ r^4 j df / dr \f$ at the surface
508  double last_rjw;
509  /// The value of \f$ \bar{\omega} / \Omega \f$ at the surface
510  double last_f;
511  /** \brief The value of \f$ d (\bar{\omega}/\Omega)/dr \f$
512  at the surface (when \ref ang_vel is true)
513  */
514  double domega_rat;
515 
516  /** \brief Maximum value for central pressure in
517  \f$ \mathrm{M}_{\odot}/\mathrm{km}^3 \f$ (default \f$ 10^{20} \f$ )
518 
519  This variable is set by the <tt>mvsr()</tt>, <tt>max()</tt>,
520  <tt>fixed()</tt> and <tt>fixed_pr()</tt>
521  functions and used in \ref integ_star() .
522  */
523  double pcent_max;
524  //@}
525 
526  /** \name Solution parameters
527 
528  These parameters can be changed at any time.
529  */
530  //@{
531  /** \brief If true, then fixed() and max() reformat results into
532  a \ref o2scl::table_units object
533 
534  Note that \ref mvsr() always places its results in the
535  output table independent of the value of this variable.
536  */
538  /** \brief The mass of one baryon
539 
540  The mass of one baryon in kg for the total baryon mass
541  calculation (defaults to the proton mass).
542  */
543  double baryon_mass;
544  /// If true, solve for the angular velocity (default false)
545  bool ang_vel;
546  /// Apply general relativistic corrections (default true)
547  bool gen_rel;
548  /** \brief calculate the gravitational potential (default false)
549  */
550  bool calc_gpot;
551  /// smallest allowed radial stepsize in km (default 1.0e-4)
552  double step_min;
553  /// largest allowed radial stepsize in km (default 0.05)
554  double step_max;
555  /// initial radial stepsize in km (default 4.0e-3)
556  double step_start;
557  /// control for output (default 1)
558  int verbose;
559  /// Maximum number of integration steps (default 100000)
561  /** \brief If true, call the error handler if the solution does
562  not converge (default true)
563  */
565  //@}
566 
567  /** \brief Default value of maximum pressure for maximum mass star
568  in \f$ \mathrm{M}_{\odot}/\mathrm{km}^3 \f$
569  (default \f$ 10^{20} \f$)
570  */
571  double pmax_default;
572 
573  /// \name Mass versus radius parameters
574  //@{
575  /** \brief Beginning pressure in
576  \f$ \mathrm{M}_{\odot}/\mathrm{km}^3 \f$ (default 7.0e-7)
577  */
578  double prbegin;
579  /** \brief Ending pressure in
580  \f$ \mathrm{M}_{\odot}/\mathrm{km}^3 \f$ (default 8.0e-3)
581  */
582  double prend;
583  /// Increment factor for pressure (default 1.1)
584  double princ;
585  /** \brief List of pressures at which more information should be
586  recorded
587 
588  If pressures (in the user-specified units) are added to this
589  vector, then in mvsr(), the radial location, enclosed
590  gravitational mass, and (if \ref o2scl::eos_tov::baryon_column
591  is true) enclosed baryon mass are stored in the table for each
592  central pressure. The associated columns are named
593  <tt>r0, gm0, bm0, r1, gm1, bm1,</tt> etc.
594  */
595  std::vector<double> pr_list;
596  //@}
597 
598  /// \name Fixed mass parameter
599  //@{
600  /** \brief Guess for central pressure in
601  \f$ \mathrm{M}_{\odot}/\mathrm{km}^3 \f$
602  (default \f$ 5.2 \times 10^{-5} \f$)
603 
604  This guess is used in the functions fixed().
605  */
607  //@}
608 
609  /// \name Maximum mass profile parameters
610  //@{
611  /** \brief Beginning pressure for maximum mass guess in
612  \f$ \mathrm{M}_{\odot}/\mathrm{km}^3 \f$ (default 7.0e-5)
613  */
614  double max_begin;
615  /** \brief Ending pressure for maximum mass guess in
616  \f$ \mathrm{M}_{\odot}/\mathrm{km}^3 \f$ (default 5.0e-3)
617  */
618  double max_end;
619  /** \brief Increment for pressure for maximum mass guess
620  (unitless; default 1.3)
621  */
622  double max_inc;
623  //@}
624 
625  /// \name Basic operation
626  //@{
627  /// Set the EOS to use
628  void set_eos(eos_tov &ter) {
629  te=&ter;
630  eos_set=true;
631  return;
632  }
633 
634  /** \brief Set output units for the table
635  */
636  void set_units(double s_efactor=1.0, double s_pfactor=1.0,
637  double s_nbfactor=1.0);
638 
639  /** \brief Set output units for the table
640  */
641  void set_units(std::string eunits="", std::string punits="",
642  std::string nunits="");
643 
644  /** \brief Get output units for the table
645  */
646  void get_units(std::string &eunits, std::string &punits,
647  std::string &nunits);
648 
649  /** \brief Calculate the mass vs. radius curve
650  */
651  virtual int mvsr();
652 
653  /** \brief Calculate the profile of a star with fixed mass
654 
655  This function computes the profile for a star with a fixed
656  mass. If the target mass is negative, it is interpreted as
657  subtracting from the maximum mass configuration. For a 2.15
658  solar mass neutron star, <tt>target_mass=-0.2</tt> corresponds
659  to 1.95 solar mass configuration.
660 
661  The variable \c pmax is the maximum allowable central pressure
662  in \f$ \mathrm{M}_{\odot}/\mathrm{km}^3 \f$ (This often, but
663  not always, equal to the central pressure of the maximum mass
664  star.) This ensures that the function does not unintentionally
665  select an unstable configuration. If \c pmax is greater than
666  or equal to the default value (\ref pmax_default), then the
667  maximum mass star will be computed with \ref max() first
668  in order to determine the maximum allowable central pressure.
669 
670  Note that this function will likely fail when the mass-radius
671  curve has two central pressures with the same gravitational
672  mass.
673  */
674  virtual int fixed(double target_mass, double pmax=1.0e20);
675 
676  /** \brief Calculate the profile of a star with a specified
677  central pressure
678 
679  This function computes the profile of a star with a fixed
680  central pressure. The central pressure, \c pcent, should be in
681  the unit system specified by the user which defaults to solar
682  masses per cubic kilometer "Msun/km^3" but can be changed with
683  a call to one of the <tt>set_units()</tt> functions.
684 
685  The variable \c pmax is the maximum allowable central pressure
686  in \f$ \mathrm{M}_{\odot}/\mathrm{km}^3 \f$, and must
687  be larger than the value of \c pcent converted to to
688  \f$ \mathrm{M}_{\odot}/\mathrm{km}^3 \f$ .
689  */
690  virtual int fixed_pr(double pcent, double pmax=1.0e20);
691 
692  /** \brief Calculate the profile of the maximum mass star
693 
694  Note that this function maximizes the gravitational mass. If
695  the M-R curve has two stable branches, then this function does
696  not necessarily give the configuration with the largest
697  central pressure.
698 
699  This function may also depend on the accuracy of the initial
700  interval determined by \ref max_begin and \ref max_end.
701  */
702  virtual int max();
703 
704  /** \brief Construct a table from the results
705 
706  This function constructs a \ref table_units object from the
707  information in \ref rkx, \ref rky, and \ref rkdydx . Note that
708  the table is always constructed by default so this function
709  need not be called unless \ref tov_solve::reformat_results is
710  <tt>false</tt>.>
711  */
712  virtual void make_table();
713 
714  /** \brief Return the results data table
715  */
716  std::shared_ptr<table_units<> > get_results() {
717  return out_table;
718  }
719 
720  /** \brief Return the results data table
721 
722  This function immediately adds four constants to the table,
723  <tt>schwarz, Msun, pi</tt> and <tt>mproton</tt>.
724  */
725  void set_table(std::shared_ptr<table_units<> > t) {
726  out_table=t;
727  // Add native constants
728  out_table->add_constant("schwarz",schwarz_km);
729  out_table->add_constant("Msun",o2scl_mks::solar_mass);
730  out_table->add_constant("pi",o2scl_const::pi);
731  out_table->add_constant("mproton",o2scl_mks::mass_proton);
732  return;
733  }
734  //@}
735 
736  /// \name Convergence information flags
737  //@{
738  static const int fixed_solver_failed=128;
739  static const int fixed_integ_star_failed=256;
740  static const int fixed_wrong_mass=512;
741  static const int max_minimizer_failed=1024;
742  static const int max_integ_star_failed=2048;
743  static const int mvsr_integ_star_failed=4096;
744  static const int ang_vel_failed=8192;
745  static const int cent_press_large=16384;
746  static const int cent_press_neg=32768;
747  static const int over_max_steps=65536;
748  static const int last_step_large=131072;
749  //@}
750 
751  /// \name Get internally formatted results (in internal unit system)
752  //@{
753  /// Get the vector for the radial grid
754  const ubvector &get_rkx() const {
755  return rkx;
756  }
757  /// Get a list of vectors for the ODEs
758  const std::vector<ubvector> &get_rky() const {
759  return rky;
760  }
761  /// Get a list of vectors for the corresponding derivatives
762  const std::vector<ubvector> &get_rkdydx() const {
763  return rkdydx;
764  }
765  //@}
766 
767  /// \name Control numerical methods
768  //@{
769  /** \brief Set solver
770  */
772  mroot_ptr=&s_mrp;
773  return;
774  };
775 
776  /** \brief Set minimizer
777  */
778  void set_min(min_base<> &s_mp) {
779  min_ptr=&s_mp;
780  return;
781  };
782 
783  /** \brief Set the adaptive stepper
784  */
786  as_ptr=&sap;
787  return;
788  };
789  //@}
790 
791  /// \name Default numerical objects
792  //@{
793  /// The default minimizer
795 
796  /// The default solver
798 
799  /// The default adaptive stepper
801  //@}
802 
803  };
804 
805 #ifndef DOXYGEN_NO_O2NS
806 }
807 #endif
808 
809 #endif
810 
811 
void set_eos(eos_tov &ter)
Set the EOS to use.
Definition: tov_solve.h:628
virtual int mvsr()
Calculate the mass vs. radius curve.
const std::vector< ubvector > & get_rky() const
Get a list of vectors for the ODEs.
Definition: tov_solve.h:758
double pcent_max
Maximum value for central pressure in (default )
Definition: tov_solve.h:523
std::string punits
Units for pressure.
Definition: tov_solve.h:416
double rad
Radius (in km)
Definition: tov_solve.h:502
const double solar_mass
double schwarz_km
The schwarzchild radius in km.
Definition: tov_solve.h:385
mroot_hybrids< mm_funct, ubvector, ubmatrix, jac_funct > def_solver
The default solver.
Definition: tov_solve.h:797
const double mass_proton
bool calc_gpot
calculate the gravitational potential (default false)
Definition: tov_solve.h:550
virtual int fixed(double target_mass, double pmax=1.0e20)
Calculate the profile of a star with fixed mass.
bool err_nonconv
If true, call the error handler if the solution does not converge (default true)
Definition: tov_solve.h:564
void set_mroot(mroot< mm_funct, ubvector, jac_funct > &s_mrp)
Set solver.
Definition: tov_solve.h:771
const double pi
std::shared_ptr< table_units<> > get_results()
Return the results data table.
Definition: tov_solve.h:716
double domega_rat
The value of at the surface (when ang_vel is true)
Definition: tov_solve.h:514
std::string nunits
Units for baryon density.
Definition: tov_solve.h:418
double max_begin
Beginning pressure for maximum mass guess in (default 7.0e-5)
Definition: tov_solve.h:614
double pmax_default
Default value of maximum pressure for maximum mass star in (default )
Definition: tov_solve.h:571
size_t max_integ_steps
Maximum number of integration steps (default 100000)
Definition: tov_solve.h:560
int verbose
control for output (default 1)
Definition: tov_solve.h:558
size_t ix_last
The last row index in rky.
Definition: tov_solve.h:382
bool gen_rel
Apply general relativistic corrections (default true)
Definition: tov_solve.h:547
double max_end
Ending pressure for maximum mass guess in (default 5.0e-3)
Definition: tov_solve.h:618
double prbegin
Beginning pressure in (default 7.0e-7)
Definition: tov_solve.h:578
void set_units(double s_efactor=1.0, double s_pfactor=1.0, double s_nbfactor=1.0)
Set output units for the table.
double pfactor
unit conversion factor for pressure (default 1.0)
Definition: tov_solve.h:422
std::vector< ubvector > rky
ODE functions.
Definition: tov_solve.h:451
virtual int integ_star(size_t ndvar, const ubvector &ndx, ubvector &ndy)
The solver function to compute the stellar profile.
double min_log_pres
Smallest allowed pressure for integration (default: -100)
Definition: tov_solve.h:434
void get_units(std::string &eunits, std::string &punits, std::string &nunits)
Get output units for the table.
double prend
Ending pressure in (default 8.0e-3)
Definition: tov_solve.h:582
bool integ_star_final
If true, integ_star() computes all the profile info, otherwise it only computes the gravitational mas...
Definition: tov_solve.h:379
const std::vector< ubvector > & get_rkdydx() const
Get a list of vectors for the corresponding derivatives.
Definition: tov_solve.h:762
double bmass
Baryonic mass (when computed; in )
Definition: tov_solve.h:504
bool ang_vel
If true, solve for the angular velocity (default false)
Definition: tov_solve.h:545
bool eos_set
True if the EOS has been set.
Definition: tov_solve.h:408
virtual void make_table()
Construct a table from the results.
min_brent_gsl def_min
The default minimizer.
Definition: tov_solve.h:794
astep_gsl< ubvector, ubvector, ubvector, ode_funct > def_stepper
The default adaptive stepper.
Definition: tov_solve.h:800
virtual int derivs(double x, size_t nv, const ubvector &y, ubvector &dydx)
The ODE step function.
eos_tov * te
The EOS.
Definition: tov_solve.h:405
virtual double max_fun(double maxx)
The minimizer function to compute the maximum mass.
size_t buffer_size
Size of the ODE solution buffer (default )
Definition: tov_solve.h:492
bool reformat_results
If true, then fixed() and max() reformat results into a o2scl::table_units object.
Definition: tov_solve.h:537
std::string eunits
Units for energy density.
Definition: tov_solve.h:414
std::function< int(double, size_t, const boost::numeric::ublas::vector< double > &, boost::numeric::ublas::vector< double > &)> ode_funct
std::shared_ptr< table_units<> > out_table
The output table.
Definition: tov_solve.h:457
void set_stepper(astep_base< ubvector, ubvector, ubvector, ode_funct > &sap)
Set the adaptive stepper.
Definition: tov_solve.h:785
A EOS base class for the TOV solver.
Definition: eos_tov.h:47
double mass
Gravitational mass (in )
Definition: tov_solve.h:500
const ubvector & get_rkx() const
Get the vector for the radial grid.
Definition: tov_solve.h:754
double baryon_mass
The mass of one baryon.
Definition: tov_solve.h:543
virtual int fixed_pr(double pcent, double pmax=1.0e20)
Calculate the profile of a star with a specified central pressure.
double step_max
largest allowed radial stepsize in km (default 0.05)
Definition: tov_solve.h:554
ubvector rkx
Radial coordinate (in kilometers)
Definition: tov_solve.h:439
void column_setup(bool mvsr_mode=false)
Set up column names and units.
double last_rjw
The value of at the surface.
Definition: tov_solve.h:508
std::vector< ubvector > rkdydx
The derivatives of the ODE functions.
Definition: tov_solve.h:453
Class to solve the Tolman-Oppenheimer-Volkov equations.
Definition: tov_solve.h:350
interp< ubvector > iop
Interpolation object for listed radii in mvsr()
Definition: tov_solve.h:367
ode_funct ofm
ODE function object.
Definition: tov_solve.h:364
double max_inc
Increment for pressure for maximum mass guess (unitless; default 1.3)
Definition: tov_solve.h:622
double nfactor
unit conversion factor for baryon density (default 1.0)
Definition: tov_solve.h:424
double fixed_pr_guess
Guess for central pressure in (default )
Definition: tov_solve.h:606
double princ
Increment factor for pressure (default 1.1)
Definition: tov_solve.h:584
void set_table(std::shared_ptr< table_units<> > t)
Return the results data table.
Definition: tov_solve.h:725
double step_start
initial radial stepsize in km (default 4.0e-3)
Definition: tov_solve.h:556
mroot< mm_funct, ubvector, jac_funct > * mroot_ptr
The solver for fixed gravitational masses.
Definition: tov_solve.h:462
min_base * min_ptr
The minimizer for maximum masses.
Definition: tov_solve.h:465
double tmass
Target mass for integ_star()
Definition: tov_solve.h:392
void make_unique_name(std::string &col, std::vector< std::string > &cnames)
Ensure col does not match strings in cnames.
double step_min
smallest allowed radial stepsize in km (default 1.0e-4)
Definition: tov_solve.h:552
virtual int max()
Calculate the profile of the maximum mass star.
std::vector< double > pr_list
List of pressures at which more information should be recorded.
Definition: tov_solve.h:595
size_t max_table_size
Maximum number of lines in the output table (default 400)
Definition: tov_solve.h:495
double last_f
The value of at the surface.
Definition: tov_solve.h:510
double gpot
Gravitational potential (when computed; unitless)
Definition: tov_solve.h:506
astep_base< ubvector, ubvector, ubvector, ode_funct > * as_ptr
The adaptive stepper.
Definition: tov_solve.h:468
void set_min(min_base<> &s_mp)
Set minimizer.
Definition: tov_solve.h:778
double efactor
unit conversion factor for energy density (default 1.0)
Definition: tov_solve.h:420

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