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
242  independent of the relative scale for \f$ \bar{\omega} \f$ and
243  \f$ j \f$ . First, one rescales \f$ \bar{\omega} \f$ and
244  rewrites everything in terms of \f$ f\equiv \bar{\omega}/\Omega
245  \f$ and \f$ g \equiv r^4 j~df/dr \f$ . Then, pick a central
246  pressure, \f$ m(r=0)=g(r=0)=0 \f$, arbitrary values for \f$
247  \Phi(r=0) \f$ and \f$ f(r=0) \f$, and integrate
248  \f{eqnarray*}
249  \frac{d P}{d r} &=& - \frac{G \varepsilon m}{r^2}
250  \left( 1+\frac{P}{\varepsilon}\right)
251  \left( 1+\frac{4 \pi P r^3}{m} \right)
252  \left( 1-\frac{2 G m}{r}\right)^{-1}
253  \nonumber \\
254  \frac{d m}{d r} &=& 4 \pi r^2 \varepsilon
255  \nonumber \\
256  \frac{d \Phi}{d r} &=& - \frac{1}{\varepsilon}
257  \frac{ d P}{d r} \left(1+\frac{P}{\varepsilon}\right)^{-1}
258  \nonumber \\
259  \frac{d g}{dr} &=& -4 r^3 \frac{d j}{dr} f
260  \nonumber \\
261  \frac{d f}{dr} &=& \frac{g}{r^4 j}
262  \f}
263  Afterwards, shift \f$ \Phi \f$ by a constant to ensure
264  the correct value at \f$ r=R \f$, and multiply \f$ g \f$
265  by a constant to ensure that \f$ g=r^4 j (df/dr) \f$ holds
266  for the new potential \f$ \Phi \f$. Then, multiply
267  \f$ f \f$ by a constant to ensure that
268  \f[
269  f(r=R) + \frac{R}{3} \left(\frac{df}{dr}\right)_{r=R} = 1
270  \f]
271 
272  The functions \f$ f \f$ and \f$ g \f$ are stored in columns
273  called <tt>"omega_rat"</tt> and <tt>"rjw"</tt>, respectively.
274  One can compute the baryonic mass by integration or by adding
275  one additional differential equation, bringing the total to six.
276 
277  The moment of inertia is (see \ref tovtoc),
278  \f[
279  I = \frac{R^4}{6 G} \left.\frac{df}{dr}\right|_{r=R}
280  \f]
281  where the last fraction is stored in \ref domega_rat .
282  For an object named \c ts of type \ref tov_solve
283  after a call to \ref fixed(), \ref max(), or
284  \ref fixed_pr(), the moment of inertia can be computed
285  with, e.g.
286  \code
287  double schwarz_km=o2scl_cgs::schwarzchild_radius/1.0e5;
288  double I=ts.domega_rat*pow(ts.rad,4.0)/3.0/schwarz_km;
289  \endcode
290 
291  <b>Convergence details</b>
292 
293  By default, if the TOV solver fails to converge, the error
294  handler is called and an exception is thrown. If \ref
295  o2scl::tov_solve::err_nonconv is false, then \ref
296  o2scl::tov_solve::mvsr(), \ref o2scl::tov_solve::fixed(), and
297  \ref o2scl::tov_solve::max(), return an integer which gives some
298  information about why the solver failed to converge.
299 
300  If \ref err_nonconv is set to false, then the \ref fixed()
301  function temporarily sets the value of both \ref
302  o2scl::mroot::err_nonconv for the current solver and \ref
303  o2scl::jacobian::err_nonconv for the jacobian object, \ref
304  o2scl::mroot_hybrids::def_jac, of \ref def_solver equal to false.
305 
306  <b>Other details</b>
307 
308  The ODE solution is stored in a buffer which can be directly
309  accessed using \ref o2scl::tov_solve::get_rkx(), \ref
310  o2scl::tov_solve::get_rky(), and \ref
311  o2scl::tov_solve::get_rkdydx(). In the case that the ODE steps
312  are small enough that there isn't enough space in the buffers,
313  the ODE solution data is thinned by a factor of two to allow for
314  the remaining ODE steps to be stored. The size of the buffers
315  can be controlled in \ref buffer_size .
316 
317  If \ref o2scl::tov_solve::reformat_results is true (the
318  default), then the results are placed in a shared pointer to the
319  \ref table_units object which can be accessed using \ref
320  o2scl::tov_solve::get_results(). The maximum size of the output
321  table can be controlled with \ref max_table_size. The output
322  table may be smaller than this, as it cannot be larger than the
323  number of steps stored in the buffer.
324 
325  \note The function \ref o2scl::tov_solve::integ_star() returns
326  <tt>gsl_efailed</tt> without calling the error handler in the
327  case that the solver can recover gracefully from, for example, a
328  negative pressure.
329 
330  \future Convert to \ref o2scl::ode_iv_solve?
331  */
332  class tov_solve {
333 
334  public:
335 
338  typedef boost::numeric::ublas::matrix_column<ubmatrix> ubmatrix_column;
339  typedef boost::numeric::ublas::matrix_row<ubmatrix> ubmatrix_row;
340 
341 #ifndef DOXYGEN_INTERNAL
342 
343  protected:
344 
345  /// ODE function object
347 
348  /// Interpolation object for listed radii in \ref mvsr()
350 
351  /** \brief Set up column names and units
352 
353  When this function is used by \ref mvsr(), \c mvsr_mode is set
354  to true.
355  */
356  void column_setup(bool mvsr_mode=false);
357 
358  /** \brief If true, \ref integ_star() computes all the profile info,
359  otherwise it only computes the gravitational mass
360  */
362 
363  /// The last row index in \ref rky
364  size_t ix_last;
365 
366  /// The schwarzchild radius in km
367  double schwarz_km;
368 
369  /** \brief Target mass for integ_star()
370 
371  Has a value of zero, unless set to a non-zero value by \ref fixed().
372  */
373  double tmass;
374 
375  /** \brief Ensure \c col does not match strings in \c cnames
376 
377  Underscores are added to \c col until it matches none of
378  the strings in \c cnames.
379  */
380  void make_unique_name(std::string &col,
381  std::vector<std::string> &cnames);
382 
383  /// \name User EOS
384  //@{
385  /// The EOS
387 
388  /// True if the EOS has been set
389  bool eos_set;
390  //@}
391 
392  /// \name Units for output table
393  //@{
394  /// Units for energy density
395  std::string eunits;
396  /// Units for pressure
397  std::string punits;
398  /// Units for baryon density
399  std::string nunits;
400  /// unit conversion factor for energy density (default 1.0)
401  double efactor;
402  /// unit conversion factor for pressure (default 1.0)
403  double pfactor;
404  /// unit conversion factor for baryon density (default 1.0)
405  double nfactor;
406  //@}
407 
408  /** \brief Smallest allowed pressure for integration (default: -100)
409 
410  This quantity can't be much smaller than -100 since we need to
411  compute numbers near \f$ e^{-\mathrm{min\_log\_pres}} \f$
412 
413  \future Replace this with the proper value from std::limits?
414  */
415  double min_log_pres;
416 
417  /// \name Integration storage
418  //@{
419  /// Radial coordinate (in kilometers)
420  ubvector rkx;
421  /** \brief ODE functions
422 
423  If \c rky is viewed as a matrix, then the first column of each
424  row is the gravitational mass in solar masses, and the second
425  column is the natural logarithm of the pressure in \f$
426  \mathrm{M}_{\odot}/km^3 \f$ . When \ref calc_gpot is true, the
427  next column is the gravitational potential (which is
428  unitless), and when \ref eos_tov::baryon_column is true, the
429  next column is the baryonic mass in \f$ \mathrm{M}_{\odot}
430  \f$.
431  */
432  std::vector<ubvector> rky;
433  /// The derivatives of the ODE functions
434  std::vector<ubvector> rkdydx;
435  //@}
436 
437  /// The output table
438  std::shared_ptr<table_units<> > out_table;
439 
440  /// \name Numerical methods
441  //@{
442  /// The solver
444 
445  /// The minimizer
447 
448  /// The adaptive stepper
450  //@}
451 
452  /// The ODE step function
453  virtual int derivs(double x, size_t nv, const ubvector &y,
454  ubvector &dydx);
455 
456  /// The minimizer function to compute the maximum mass
457  virtual double max_fun(double maxx);
458 
459  /** \brief The solver function to compute the stellar profile
460  */
461  virtual int integ_star(size_t ndvar, const ubvector &ndx,
462  ubvector &ndy);
463 
464 #endif
465 
466  public:
467 
468  tov_solve();
469 
470  virtual ~tov_solve();
471 
472  /// Size of the ODE solution buffer (default \f$ 10^5 \f$)
473  size_t buffer_size;
474 
475  /// Maximum number of lines in the output table (default 400)
477 
478  /// \name Basic properties
479  //@{
480  /// Gravitational mass (in \f$ \mathrm{M}_{\odot} \f$)
481  double mass;
482  /// Radius (in km)
483  double rad;
484  /// Baryonic mass (when computed; in \f$ \mathrm{M}_{\odot} \f$)
485  double bmass;
486  /// Gravitational potential (when computed; unitless)
487  double gpot;
488  /// The value of \f$ r^4 j df / dr \f$ at the surface
489  double last_rjw;
490  /// The value of \f$ \bar{\omega} / \Omega \f$ at the surface
491  double last_f;
492  /** \brief The value of \f$ d (\bar{\omega}/\Omega)/dr \f$
493  at the surface (when \ref ang_vel is true)
494  */
495  double domega_rat;
496  /** \brief Maximum value for central pressure in
497  \f$ \mathrm{M}_{\odot}/\mathrm{km}^3 \f$ (default \f$ 10^{20} \f$ )
498 
499  This variable is set by the <tt>mvsr()</tt> and <tt>max()</tt>
500  functions and used in \ref integ_star() .
501  */
502  double pcent_max;
503  //@}
504 
505  /** \name Solution parameters
506 
507  These parameters can be changed at any time.
508  */
509  //@{
510  /** \brief If true, then fixed() and max() reformat results into
511  a \ref o2scl::table_units object
512 
513  Note that \ref mvsr() always places its results in the
514  output table independent of the value of this variable.
515  */
517  /** \brief The mass of one baryon
518 
519  The mass of one baryon in kg for the total baryon mass
520  calculation (defaults to the proton mass).
521  */
522  double baryon_mass;
523  /// If true, solve for the angular velocity (default false)
524  bool ang_vel;
525  /// Apply general relativistic corrections (default true)
526  bool gen_rel;
527  /** \brief calculate the gravitational potential (default false)
528  */
529  bool calc_gpot;
530  /// smallest allowed radial stepsize in km (default 1.0e-4)
531  double step_min;
532  /// largest allowed radial stepsize in km (default 0.05)
533  double step_max;
534  /// initial radial stepsize in km (default 4.0e-3)
535  double step_start;
536  /// control for output (default 1)
537  int verbose;
538  /// Maximum number of integration steps (default 100000)
540  /** \brief If true, call the error handler if the solution does
541  not converge (default true)
542  */
544  //@}
545 
546  /** \brief Default value of maximum pressure for maximum mass star
547  in \f$ \mathrm{M}_{\odot}/\mathrm{km}^3 \f$
548  (default \f$ 10^{20} \f$)
549  */
550  double pmax_default;
551 
552  /// \name Mass versus radius parameters
553  //@{
554  /** \brief Beginning pressure in
555  \f$ \mathrm{M}_{\odot}/\mathrm{km}^3 \f$ (default 7.0e-7)
556  */
557  double prbegin;
558  /** \brief Ending pressure in
559  \f$ \mathrm{M}_{\odot}/\mathrm{km}^3 \f$ (default 8.0e-3)
560  */
561  double prend;
562  /// Increment factor for pressure (default 1.1)
563  double princ;
564  /** \brief List of pressures at which more information should be
565  recorded
566 
567  If pressures (in the user-specified units) are added to this
568  vector, then in mvsr(), the radial location, enclosed
569  gravitational mass, and (if \ref o2scl::eos_tov::baryon_column
570  is true) enclosed baryon mass are stored in the table for each
571  central pressure. The associated columns are named
572  <tt>r0, gm0, bm0, r1, gm1, bm1,</tt> etc.
573  */
574  std::vector<double> pr_list;
575  //@}
576 
577  /// \name Fixed mass parameter
578  //@{
579  /** \brief Guess for central pressure in
580  \f$ \mathrm{M}_{\odot}/\mathrm{km}^3 \f$
581  (default \f$ 5.2 \times 10^{-5} \f$)
582 
583  This guess is used in the functions fixed().
584  */
586  //@}
587 
588  /// \name Maximum mass profile parameters
589  //@{
590  /** \brief Beginning pressure for maximum mass guess in
591  \f$ \mathrm{M}_{\odot}/\mathrm{km}^3 \f$ (default 7.0e-5)
592  */
593  double max_begin;
594  /** \brief Ending pressure for maximum mass guess in
595  \f$ \mathrm{M}_{\odot}/\mathrm{km}^3 \f$ (default 5.0e-3)
596  */
597  double max_end;
598  /** \brief Increment for pressure for maximum mass guess
599  (unitless; default 1.3)
600  */
601  double max_inc;
602  //@}
603 
604  /// \name Basic operation
605  //@{
606  /// Set the EOS to use
607  void set_eos(eos_tov &ter) {
608  te=&ter;
609  eos_set=true;
610  return;
611  }
612 
613  /** \brief Set output units for the table
614  */
615  void set_units(double s_efactor=1.0, double s_pfactor=1.0,
616  double s_nbfactor=1.0);
617 
618  /** \brief Set output units for the table
619  */
620  void set_units(std::string eunits="", std::string punits="",
621  std::string nunits="");
622 
623  /// Calculate the mass vs. radius curve
624  virtual int mvsr();
625 
626  /** \brief Calculate the profile of a star with fixed mass
627 
628  If the target mass is negative, it is interpreted as
629  subtracting from the maximum mass configuration. For a 2.15
630  solar mass neutron star, <tt>target_mass=-0.2</tt> corresponds
631  to 1.95 solar mass configuration.
632 
633  The variable \c pmax is the maximum allowable central pressure
634  in \f$ \mathrm{M}_{\odot}/\mathrm{km}^3 \f$, i.e. the central
635  pressure of the maximum mass star. This ensures that the
636  function does not unintentionally select a configuration on an
637  unstable branch. If \c pmax is greater than or equal to the
638  default value (\ref pmax_default), then the maximum mass star
639  will be explicitly computed first.
640  */
641  virtual int fixed(double target_mass, double pmax=1.0e20);
642 
643  /** \brief Calculate the profile of a star with a specified
644  central pressure
645  */
646  virtual int fixed_pr(double pcent);
647 
648  /** \brief Calculate the profile of the maximum mass star
649  */
650  virtual int max();
651 
652  /** \brief Construct a table from the results
653 
654  This function constructs a \ref table_units object from the
655  information in \ref rkx, \ref rky, and \ref rkdydx . Note that
656  the table is always constructed by default so this function
657  need not be called unless \ref tov_solve::reformat_results is
658  <tt>false</tt>.>
659  */
660  virtual void make_table();
661 
662  /** \brief Return the results data table
663  */
664  std::shared_ptr<table_units<> > get_results() {
665  return out_table;
666  }
667 
668  /** \brief Return the results data table
669 
670  This function immediately adds four constants to the table,
671  <tt>schwarz, Msun, pi</tt> and <tt>mproton</tt>.
672  */
673  void set_table(std::shared_ptr<table_units<> > t) {
674  out_table=t;
675  // Add native constants
676  out_table->add_constant("schwarz",schwarz_km);
677  out_table->add_constant("Msun",o2scl_mks::solar_mass);
678  out_table->add_constant("pi",o2scl_const::pi);
679  out_table->add_constant("mproton",o2scl_mks::mass_proton);
680  return;
681  }
682  //@}
683 
684  /// \name Convergence information flags
685  //@{
686  static const int fixed_solver_failed=128;
687  static const int fixed_integ_star_failed=256;
688  static const int fixed_wrong_mass=512;
689  static const int max_minimizer_failed=1024;
690  static const int max_integ_star_failed=2048;
691  static const int mvsr_integ_star_failed=4096;
692  static const int ang_vel_failed=8192;
693  static const int cent_press_large=16384;
694  static const int cent_press_neg=32768;
695  static const int over_max_steps=65536;
696  static const int last_step_large=131072;
697  //@}
698 
699  /// \name Get internally formatted results (in internal unit system)
700  //@{
701  /// Get the vector for the radial grid
702  const ubvector &get_rkx() const {
703  return rkx;
704  }
705  /// Get a list of vectors for the ODEs
706  const std::vector<ubvector> &get_rky() const {
707  return rky;
708  }
709  /// Get a list of vectors for the corresponding derivatives
710  const std::vector<ubvector> &get_rkdydx() const {
711  return rkdydx;
712  }
713  //@}
714 
715  /// \name Control numerical methods
716  //@{
717  /** \brief Set solver
718  */
720  mroot_ptr=&s_mrp;
721  return;
722  };
723 
724  /** \brief Set minimizer
725  */
726  void set_min(min_base<> &s_mp) {
727  min_ptr=&s_mp;
728  return;
729  };
730 
731  /** \brief Set the adaptive stepper
732  */
734  as_ptr=&sap;
735  return;
736  };
737  //@}
738 
739  /// \name Default numerical objects
740  //@{
741  /// The default minimizer
743 
744  /// The default solver
746 
747  /// The default adaptive stepper
749  //@}
750 
751  };
752 
753 #ifndef DOXYGEN_NO_O2NS
754 }
755 #endif
756 
757 #endif
758 
759 
void set_eos(eos_tov &ter)
Set the EOS to use.
Definition: tov_solve.h:607
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:706
double pcent_max
Maximum value for central pressure in (default )
Definition: tov_solve.h:502
std::string punits
Units for pressure.
Definition: tov_solve.h:397
double rad
Radius (in km)
Definition: tov_solve.h:483
const double solar_mass
double schwarz_km
The schwarzchild radius in km.
Definition: tov_solve.h:367
mroot_hybrids< mm_funct, ubvector, ubmatrix, jac_funct > def_solver
The default solver.
Definition: tov_solve.h:745
const double mass_proton
bool calc_gpot
calculate the gravitational potential (default false)
Definition: tov_solve.h:529
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:543
void set_mroot(mroot< mm_funct, ubvector, jac_funct > &s_mrp)
Set solver.
Definition: tov_solve.h:719
const double pi
std::shared_ptr< table_units<> > get_results()
Return the results data table.
Definition: tov_solve.h:664
double domega_rat
The value of at the surface (when ang_vel is true)
Definition: tov_solve.h:495
std::string nunits
Units for baryon density.
Definition: tov_solve.h:399
double max_begin
Beginning pressure for maximum mass guess in (default 7.0e-5)
Definition: tov_solve.h:593
double pmax_default
Default value of maximum pressure for maximum mass star in (default )
Definition: tov_solve.h:550
size_t max_integ_steps
Maximum number of integration steps (default 100000)
Definition: tov_solve.h:539
int verbose
control for output (default 1)
Definition: tov_solve.h:537
size_t ix_last
The last row index in rky.
Definition: tov_solve.h:364
bool gen_rel
Apply general relativistic corrections (default true)
Definition: tov_solve.h:526
double max_end
Ending pressure for maximum mass guess in (default 5.0e-3)
Definition: tov_solve.h:597
double prbegin
Beginning pressure in (default 7.0e-7)
Definition: tov_solve.h:557
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:403
std::vector< ubvector > rky
ODE functions.
Definition: tov_solve.h:432
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:415
virtual int fixed_pr(double pcent)
Calculate the profile of a star with a specified central pressure.
double prend
Ending pressure in (default 8.0e-3)
Definition: tov_solve.h:561
bool integ_star_final
If true, integ_star() computes all the profile info, otherwise it only computes the gravitational mas...
Definition: tov_solve.h:361
const std::vector< ubvector > & get_rkdydx() const
Get a list of vectors for the corresponding derivatives.
Definition: tov_solve.h:710
double bmass
Baryonic mass (when computed; in )
Definition: tov_solve.h:485
bool ang_vel
If true, solve for the angular velocity (default false)
Definition: tov_solve.h:524
bool eos_set
True if the EOS has been set.
Definition: tov_solve.h:389
virtual void make_table()
Construct a table from the results.
min_brent_gsl def_min
The default minimizer.
Definition: tov_solve.h:742
astep_gsl< ubvector, ubvector, ubvector, ode_funct > def_stepper
The default adaptive stepper.
Definition: tov_solve.h:748
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:386
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:473
bool reformat_results
If true, then fixed() and max() reformat results into a o2scl::table_units object.
Definition: tov_solve.h:516
std::string eunits
Units for energy density.
Definition: tov_solve.h:395
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:438
void set_stepper(astep_base< ubvector, ubvector, ubvector, ode_funct > &sap)
Set the adaptive stepper.
Definition: tov_solve.h:733
A EOS base class for the TOV solver.
Definition: eos_tov.h:47
double mass
Gravitational mass (in )
Definition: tov_solve.h:481
const ubvector & get_rkx() const
Get the vector for the radial grid.
Definition: tov_solve.h:702
double baryon_mass
The mass of one baryon.
Definition: tov_solve.h:522
double step_max
largest allowed radial stepsize in km (default 0.05)
Definition: tov_solve.h:533
ubvector rkx
Radial coordinate (in kilometers)
Definition: tov_solve.h:420
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:489
std::vector< ubvector > rkdydx
The derivatives of the ODE functions.
Definition: tov_solve.h:434
Class to solve the Tolman-Oppenheimer-Volkov equations.
Definition: tov_solve.h:332
interp< ubvector > iop
Interpolation object for listed radii in mvsr()
Definition: tov_solve.h:349
ode_funct ofm
ODE function object.
Definition: tov_solve.h:346
double max_inc
Increment for pressure for maximum mass guess (unitless; default 1.3)
Definition: tov_solve.h:601
double nfactor
unit conversion factor for baryon density (default 1.0)
Definition: tov_solve.h:405
double fixed_pr_guess
Guess for central pressure in (default )
Definition: tov_solve.h:585
double princ
Increment factor for pressure (default 1.1)
Definition: tov_solve.h:563
void set_table(std::shared_ptr< table_units<> > t)
Return the results data table.
Definition: tov_solve.h:673
double step_start
initial radial stepsize in km (default 4.0e-3)
Definition: tov_solve.h:535
mroot< mm_funct, ubvector, jac_funct > * mroot_ptr
The solver.
Definition: tov_solve.h:443
min_base * min_ptr
The minimizer.
Definition: tov_solve.h:446
double tmass
Target mass for integ_star()
Definition: tov_solve.h:373
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:531
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:574
size_t max_table_size
Maximum number of lines in the output table (default 400)
Definition: tov_solve.h:476
double last_f
The value of at the surface.
Definition: tov_solve.h:491
double gpot
Gravitational potential (when computed; unitless)
Definition: tov_solve.h:487
astep_base< ubvector, ubvector, ubvector, ode_funct > * as_ptr
The adaptive stepper.
Definition: tov_solve.h:449
void set_min(min_base<> &s_mp)
Set minimizer.
Definition: tov_solve.h:726
double efactor
unit conversion factor for energy density (default 1.0)
Definition: tov_solve.h:401

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