nstar_rot.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  This file is part of O2scl. It has been adapted from RNS v1.1d
5  written by N. Stergioulas and S. Morsink. The modifications made in
6  this version from the original are copyright (C) 2015-2017, Andrew
7  W. Steiner.
8 
9  O2scl is free software; you can redistribute it and/or modify
10  it under the terms of the GNU General Public License as published by
11  the Free Software Foundation; either version 3 of the License, or
12  (at your option) any later version.
13 
14  O2scl is distributed in the hope that it will be useful,
15  but WITHOUT ANY WARRANTY; without even the implied warranty of
16  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17  GNU General Public License for more details.
18 
19  You should have received a copy of the GNU General Public License
20  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
21 
22  -------------------------------------------------------------------
23 */
24 /*
25  -------------------------------------------------------------------
26  Relativistic models of rapidly rotating compact stars,
27  using tabulated or polytropic equations of state.
28 
29  Author: Nikolaos Stergioulas
30 
31  Current Address:
32 
33  Department of Physics
34  University of Wisconsin-Milwaukee
35  PO Box 413, Milwaukee, WI 53201, USA
36 
37  E-mail: niksterg@csd.uwm.edu, or
38  niksterg@pauli.phys.uwm.edu
39 
40  Version: 1.1
41 
42  Date: June, 1995
43 
44  Changes made to code by Sharon Morsink
45 
46  03-03-97: Corrected the units for polytropic stars
47  10-28-98: Added the star's quadrupole moment to the output.
48 
49  References:
50  KEH : H. Komatsu, Y. Eriguchi and I. Hachisu, Mon. Not. R. astr. Soc.
51  (1989) 237, 355-379.
52  CST : G. Cook, S. Shapiro and S. Teukolsky, Ap. J (1992) 398, 203-223.
53 
54  -------------------------------------------------------------------
55 */
56 /** \file nstar_rot.h
57  \brief File defining \ref o2scl::nstar_rot
58 */
59 #ifndef NSTAR_ROT_H
60 #define NSTAR_ROT_H
61 
62 #include <cmath>
63 #include <iostream>
64 
65 #include <boost/numeric/ublas/vector.hpp>
66 #include <boost/numeric/ublas/matrix.hpp>
67 
68 #include <o2scl/err_hnd.h>
69 #include <o2scl/search_vec.h>
70 #include <o2scl/test_mgr.h>
71 #include <o2scl/root_bkt_cern.h>
72 #include <o2scl/lib_settings.h>
73 #include <o2scl/interp.h>
74 #include <o2scl/eos_tov.h>
75 #include <o2scl/table3d.h>
76 #include <o2scl/tensor.h>
77 #include <o2scl/mroot_hybrids.h>
78 
79 namespace o2scl {
80 
81  /** \brief An EOS for \ref nstar_rot
82  */
83  class eos_nstar_rot : public eos_tov {
84 
85  public:
86 
87  /** \brief From the pressure, return the enthalpy
88  */
89  virtual double enth_from_pr(double pr)=0;
90 
91  /** \brief From the enthalpy, return the pressure
92  */
93  virtual double pr_from_enth(double enth)=0;
94 
95  /** \brief From the baryon density, return the enthalpy
96  */
97  virtual double enth_from_nb(double nb)=0;
98  };
99 
100  /** \brief Create a tabulated EOS for \ref nstar_rot using interpolation
101 
102  \todo Document what the \ref nstar_rot code requirements are
103  for the low-density part of the EOS.
104 
105  \future Replace arrays with vectors
106  */
108 
109  protected:
110 
111  /// Array search object
113 
114  /// Search in array \c x of length \c n for value \c val
115  int new_search(int n, double *x, double val);
116 
117  /** \brief number of tabulated EOS points */
118  int n_tab;
119  /** \brief rho points in tabulated EOS */
120  double log_e_tab[201];
121  /** \brief p points in tabulated EOS */
122  double log_p_tab[201];
123  /** \brief h points in EOS file */
124  double log_h_tab[201];
125  /** \brief number density in EOS file */
126  double log_n0_tab[201];
127 
128  /// \name Constants
129  //@{
130  /** \brief Speed of light in vacuum (in CGS units) */
131  double C;
132  /** \brief Gravitational constant (in CGS units) */
133  double G;
134  /** \brief Square of length scale in CGS units,
135  \f$ \kappa \equiv 10^{-15} c^2/G \f$
136  */
137  double KAPPA;
138  /** \brief The value \f$ 10^{-15}/c^2 \f$ */
139  double KSCALE;
140  //@}
141 
142  /** \brief Cache for interpolation
143  */
145 
148 
149  /** \brief Driver for the interpolation routine.
150 
151  First we find the tab. point nearest to xb, then we
152  interpolate using four points around xb.
153  */
154  double interp(double xp[], double yp[], int np, double xb);
155 
156  public:
157 
159 
160  /** \brief Set the EOS from four vectors in the native unit system
161  */
162  template<class vec1_t, class vec2_t, class vec3_t, class vec4_t>
163  void set_eos_native(vec1_t &eden, vec2_t &pres, vec3_t &enth,
164  vec4_t &nb) {
165 
166  double C=o2scl_cgs::speed_of_light;
168  double KAPPA=1.0e-15*C*C/G;
169  double KSCALE=KAPPA*G/(C*C*C*C);
170 
171  n_tab=eden.size();
172 
173  for(int i=1;i<=n_tab;i++) {
174  log_e_tab[i]=log10(eden[i-1]*C*C*KSCALE);
175  log_p_tab[i]=log10(pres[i-1]*KSCALE);
176  log_h_tab[i]=log10(enth[i-1]/(C*C));
177  log_n0_tab[i]=log10(nb[i-1]);
178  }
179 
180  return;
181  }
182 
183  /** \brief Set the EOS from energy density, pressure, and
184  baryon density stored in powers of \f$ \mathrm{fm} \f$ .
185  */
186  template<class vec1_t, class vec2_t, class vec3_t>
187  void set_eos_fm(size_t n, vec1_t &eden, vec2_t &pres, vec3_t &nb) {
188 
189  static const int n_crust=78;
190 
191  if (n>200-n_crust) {
192  O2SCL_ERR2("Too many EOS points in ",
193  "nstar_rot::set_eos().",o2scl::exc_einval);
194  }
195 
196  // Conversion factor for energy density
197  double conv1=o2scl_settings.get_convert_units().convert
198  ("1/fm^4","g/cm^3",1.0);
199  // Conversion factor for pressure
200  double conv2=o2scl_settings.get_convert_units().convert
201  ("1/fm^4","dyne/cm^2",1.0);
202 
203  n_tab=n+n_crust;
204 
205  /* Use the original RNS crust, except for the enthalpy which is
206  computed by hand below. This appears to work better than the
207  default O2scl crust, and this may have to do with the fact
208  that the default O2scl crust has decreasing mu with
209  increasing density at low densities.
210  */
211  double nst_arr[n_crust][3]={
212  {7.800e+00,1.010e+08,4.698795180722962e+24},
213  {7.860e+00,1.010e+09,4.734939759036205e+24},
214  {7.900e+00,1.010e+10,4.759036144578364e+24},
215  {8.150e+00,1.010e+11,4.909638554215315e+24},
216  {1.160e+01,1.210e+12,6.987951807098076e+24},
217  {1.640e+01,1.400e+13,9.879518070489597e+24},
218  {4.510e+01,1.700e+14,2.716867462904601e+25},
219  {2.120e+02,5.820e+15,1.277108403508764e+26},
220  {1.150e+03,1.900e+17,6.927709645088004e+26},
221  {1.044e+04,9.744e+18,6.289148562640985e+27},
222  {2.622e+04,4.968e+19,1.579513843816999e+28},
223  {6.587e+04,2.431e+20,3.968050678245718e+28},
224  {1.654e+05,1.151e+21,9.963748410271617e+28},
225  {4.156e+05,5.266e+21,2.503563031417219e+29},
226  {1.044e+06,2.318e+22,6.288917532113082e+29},
227  {2.622e+06,9.755e+22,1.579410809416864e+30},
228  {6.588e+06,3.911e+23,3.968207649843547e+30},
229  {8.293e+06,5.259e+23,4.995116726219748e+30},
230  {1.655e+07,1.435e+24,9.967984755458204e+30},
231  {3.302e+07,3.833e+24,1.988624478073943e+31},
232  {6.589e+07,1.006e+25,3.967807406359445e+31},
233  {1.315e+08,2.604e+25,7.917691186982454e+31},
234  {2.624e+08,6.676e+25,1.579648605894070e+32},
235  {3.304e+08,8.738e+25,1.988876577393412e+32},
236  {5.237e+08,1.629e+26,3.152005155076383e+32},
237  {8.301e+08,3.029e+26,4.995278531652059e+32},
238  {1.045e+09,4.129e+26,6.287859551784352e+32},
239  {1.316e+09,5.036e+26,7.917701445937253e+32},
240  {1.657e+09,6.860e+26,9.968319738044036e+32},
241  {2.626e+09,1.272e+27,1.579408507997411e+33},
242  {4.164e+09,2.356e+27,2.503766293549853e+33},
243  {6.601e+09,4.362e+27,3.967852390467774e+33},
244  {8.312e+09,5.662e+27,4.995474308724729e+33},
245  {1.046e+10,7.702e+27,6.285277578607203e+33},
246  {1.318e+10,1.048e+28,7.918132634568090e+33},
247  {1.659e+10,1.425e+28,9.964646988214994e+33},
248  {2.090e+10,1.938e+28,1.255052800774333e+34},
249  {2.631e+10,2.503e+28,1.579545673652798e+34},
250  {3.313e+10,3.404e+28,1.988488463504033e+34},
251  {4.172e+10,4.628e+28,2.503379640977065e+34},
252  {5.254e+10,5.949e+28,3.151720931652274e+34},
253  {6.617e+10,8.089e+28,3.968151735612910e+34},
254  {8.332e+10,1.100e+29,4.994995310195290e+34},
255  {1.049e+11,1.495e+29,6.286498800006776e+34},
256  {1.322e+11,2.033e+29,7.919521253825185e+34},
257  {1.664e+11,2.597e+29,9.964341016667146e+34},
258  {1.844e+11,2.892e+29,1.104024323001462e+35},
259  {2.096e+11,3.290e+29,1.254619611126682e+35},
260  {2.640e+11,4.473e+29,1.579588892045295e+35},
261  {3.325e+11,5.816e+29,1.988565738933728e+35},
262  {4.188e+11,7.538e+29,2.503561780689725e+35},
263  {4.299e+11,7.805e+29,2.569780082714395e+35},
264  {4.460e+11,7.890e+29,2.665824694449485e+35},
265  {5.228e+11,8.352e+29,3.123946525953616e+35},
266  {6.610e+11,9.098e+29,3.948222384313103e+35},
267  {7.964e+11,9.831e+29,4.755697604312120e+35},
268  {9.728e+11,1.083e+30,5.807556544067428e+35},
269  {1.196e+12,1.218e+30,7.138304213736713e+35},
270  {1.471e+12,1.399e+30,8.777653631971616e+35},
271  {1.805e+12,1.683e+30,1.076837272716171e+36},
272  {2.202e+12,1.950e+30,1.313417953138369e+36},
273  {2.930e+12,2.592e+30,1.747157788902558e+36},
274  {3.833e+12,3.506e+30,2.285004034820638e+36},
275  {4.933e+12,4.771e+30,2.939983642627298e+36},
276  {6.248e+12,6.481e+30,3.722722765704268e+36},
277  {7.801e+12,8.748e+30,4.646805278760175e+36},
278  {9.611e+12,1.170e+31,5.723413975645761e+36},
279  {1.246e+13,1.695e+31,7.417258934884369e+36},
280  {1.496e+13,2.209e+31,8.902909532230595e+36},
281  {1.778e+13,2.848e+31,1.057801059193907e+37},
282  {2.210e+13,3.931e+31,1.314278492046241e+37},
283  {2.988e+13,6.178e+31,1.775810743961577e+37},
284  {3.767e+13,8.774e+31,2.237518046976615e+37},
285  {5.081e+13,1.386e+32,3.015480061626022e+37},
286  {6.193e+13,1.882e+32,3.673108933334910e+37},
287  {7.732e+13,2.662e+32,4.582250451016437e+37},
288  {9.826e+13,3.897e+32,5.817514573447143e+37},
289  {1.262e+14,5.861e+32,7.462854442694524e+37}};
290 
291  double mu_start;
292  // Note that there is no c^2 needed in the computation of the
293  // enthalpy as the original code removes it.
294  for(size_t i=0;i<n_crust;i++) {
295  log_e_tab[i+1]=log10(nst_arr[i][0]*C*C*KSCALE);
296  log_p_tab[i+1]=log10(nst_arr[i][1]*KSCALE);
297  double mu=(nst_arr[i][0]/conv1+nst_arr[i][1]/conv2)/
298  nst_arr[i][2]*1.0e39;
299  if (i==0) {
300  mu_start=mu;
301  } else {
302  log_h_tab[i+1]=log10(log(mu/mu_start));
303  }
304  log_n0_tab[i+1]=log10(nst_arr[i][2]);
305  }
306  // This shift of 8.0 approximately reproduces the results
307  // implied by the internal RNS EOSs. It is not clear how
308  // sensitive the code is to the low-density part of the EOS
309  log_h_tab[1]=log_h_tab[2]-8.0;
310 
311  for(size_t i=0;i<n;i++) {
312  log_e_tab[i+n_crust+1]=log10(eden[i]*conv1*C*C*KSCALE);
313  log_p_tab[i+n_crust+1]=log10(pres[i]*conv2*KSCALE);
314  log_h_tab[i+n_crust+1]=log10(log((eden[i]+pres[i])/nb[i]/
315  mu_start));
316  log_n0_tab[i+n_crust+1]=log10(nb[i]*1.0e39);
317  }
318 
319  return;
320  }
321 
322  /** \brief From the pressure, return the energy density
323  */
324  virtual double ed_from_pr(double pr);
325 
326  /** \brief From the energy density, return the pressure
327  */
328  virtual double pr_from_ed(double ed);
329 
330  /** \brief From the pressure, return the baryon density
331  */
332  virtual double nb_from_pr(double pr);
333 
334  /** \brief From the baryon density, return the pressure
335  */
336  virtual double pr_from_nb(double nb);
337 
338  /** \brief From the baryon density, return the energy density
339  */
340  virtual double ed_from_nb(double nb);
341 
342  /** \brief From the energy density, return the baryon density
343  */
344  virtual double nb_from_ed(double ed);
345 
346  /** \brief From the pressure, return the enthalpy
347  */
348  virtual double enth_from_pr(double pr);
349 
350  /** \brief From the baryon density, return the enthalpy
351  */
352  virtual double enth_from_nb(double nb);
353 
354  /** \brief From the enthalpy, return the pressure
355  */
356  virtual double pr_from_enth(double enth);
357 
358  /** \brief Given the pressure, produce the energy and number densities
359 
360  The arguments \c pr and \c ed should always be in \f$
361  M_{\odot}/\mathrm{km}^3 \f$ . The argument for \c nb should be
362  in \f$ \mathrm{fm}^{-3} \f$ .
363 
364  If \ref baryon_column is false, then \c nb is unmodified.
365  */
366  virtual void ed_nb_from_pr(double pr, double &ed, double &nb);
367 
368  /** \brief Output EOS table to screen
369 
370  \comment
371  This is mostly for testing and should be replaced once
372  the arrays are properly replaced with vectors
373  \endcomment
374  */
375  void output() {
376  for(int i=n_tab;i>=1;i--) {
377  std::cout << log_e_tab[i] << " " << log_p_tab[i] << " "
378  << log_h_tab[i] << " " << log_n0_tab[i] << std::endl;
379  }
380  std::cout << std::endl;
381  return;
382  }
383 
384  };
385 
386  /** \brief Tabulated EOS for \ref nstar_rot from \ref Bethe74
387  */
389  public:
390  eos_nstar_rot_C(bool rns_constants=false);
391  };
392 
393  /** \brief Tabulated EOS for \ref nstar_rot from \ref Pandharipande75
394  */
396  public:
397  eos_nstar_rot_L(bool rns_constants=false);
398  };
399 
400  /** \brief Rotating neutron star class based on RNS v1.1d from
401  N. Stergioulas et al.
402 
403  \note This class is still experimental.
404 
405  Several changes have been made to the original code. The code
406  using Numerical Recipes has been removed and replaced with an
407  equivalent based on GSL and \o2. The overall interface has
408  been changed and some code has been updated with C++
409  equivalents.
410 
411  <b>Usage</b>
412 
413  <b>Initial guess</b>
414 
415  The original RNS code suggests that the initial guess is
416  typically a star with a smaller angular momentum.
417 
418  <b>References</b>
419 
420  The original RNS v1.1d can be obtained from
421  http://www.gravity.phys.uwm.edu/rns/ , and you may find Nick
422  Stergioulas's web page http://www.astro.auth.gr/~niksterg/ , or
423  Sharon Morsink's page http://fermi.phys.ualberta.ca/~morsink/
424  useful. See \ref Bonazzola73, \ref Bonazzola94, \ref Cook92,
425  \ref Cook94, \ref Friedman88, \ref Gourgoulhon94, \ref
426  Komatsu89, \ref Laarakkers99, \ref Nozawa98, \ref Stergioulas95,
427  and \ref Stergioulas03 .
428 
429  \todo Better documentation is needed everywhere.
430  \todo Test the resize() function
431 
432  \future Make a GSL-like set() function
433  \future Rework EOS interface and implement better
434  integration with the other \o2e EOSs.
435  \future Remove the unit-indexed arrays.
436  \future Try moving some of the storage to the heap?
437  \future Some of the arrays seem larger than necessary.
438  \future The function \ref o2scl::nstar_rot::new_search() is
439  inefficient because it has to handle the boundary conditions
440  separately. This could be improved.
441  \future Make the solvers more robust. The ang_vel() and ang_vel_alt()
442  functions appear particularly unstable.
443 
444  <b>Draft documentation</b>
445 
446  This class retains the definition of the specific enthalpy
447  from the original code, namely that
448  \f[
449  h = c^2 \log[\frac{\mu}{\mu(P=0)}] = \int \frac{c^2 dP}{\varepsilon+P}
450  \f]
451  but note that the factor of \f$ c^2 \f$ is dropped before
452  taking the log in the internal copy of the EOS table.
453  Typically, \f$ \mu(P=0) \f$ is around 931 MeV.
454 
455  For spherical stars, the isotropic radius \f$ r_{\mathrm{is}}
456  \f$ is defined by
457  \f[
458  \frac{d r}{d r_{\mathrm{is}}} =
459  \left(\frac{r}{r_{\mathrm{is}}}\right)
460  \left(1 - 2 \frac{m}{r}\right)^{1/2}
461  \f]
462 
463  <b>Quadrupole moments</b>
464 
465  Quadrupole moments computed using the method in \ref Laarakkers99.
466 
467  <b>Axisymmetric Instability</b>
468 
469  \ref Friedman88 shows that a secular axisymmetric instability
470  sets in when the mass becomes maximum along a sequence of
471  constant angular momentum. Equivalently, \ref Cook92 shows that
472  the instability occurs when the angular momentum becomes minimum
473  along a sequence of constant rest mass.
474 
475  A GR virial theorem for a stationary and axisymmetric system was
476  found in \ref Bonazzola73. A more general two-dimensional virial
477  identity was found in \ref Bonazzola94. The three-dimensional
478  virial identity found in \ref Gourgoulhon94 is a generalization
479  of the Newtonial virial theorem.
480 
481  Using the stationary and axisymmetric metric ( \f$ G = c = 1 \f$
482  )
483  \f[
484  ds^2 = - e^{\gamma+\rho} dt^2 + e^{2 \alpha} \left( dr^2 +
485  r^2 d\theta^2 \right) + e^{\gamma-\rho} r^2 \sin^2 \theta
486  ( d \phi - \omega dt) ^2
487  \f]
488  one solves for the four metric functions \f$ \rho(r,\theta) \f$,
489  \f$ \gamma(r,\theta) \f$, \f$ \alpha(r,\theta) \f$ and \f$
490  \omega(r,\theta) \f$ .
491 
492  It is assumed that matter is a perfect fluid, and the
493  stress-energy tensor is
494  \f[
495  T^{\mu \nu} = \left( \rho_0 + \rho_i + P \right) u^{\mu} u^{\nu}
496  + P g^{\mu \nu}
497  \f]
498 
499  Einstein's field equations imply four field equations for
500  a specified rotation law,
501  \f[
502  u^{t} u_{\phi} = F(\Omega)
503  \f]
504  for some function \f$ F(\omega) \f$ .
505 
506  Using Eq. (27) in \ref Cook92, one can write
507  \f[
508  \rho(s,\mu) = - e^{-\gamma/2} \sum_{n=0}^{\infty}
509  P_{2n}(\mu) \int_0^{1}~ds^{\prime} \int_0^1~d \mu
510  f_{\rho}(n,s,s^{\prime}) P_{2n}{\mu^{\prime}}
511  \tilde{S}(s^{\prime},\mu^{\prime})
512  \f]
513  where the function \f$ f_{\rho} \f$ is defined by
514  \f[
515  f_{\rho} \equiv \Theta(s^{\prime}-s)
516  \left(\frac{1-s}{s}\right)^{2 n+1} \left[\frac{s^{\prime
517  2n}}{(1-s^{\prime})^{2n+2}}\right] + \Theta(s^{\prime}-s)
518  \left(\frac{1-s}{s}\right)^{2 n+1} \left[\frac{s^{\prime
519  2n}}{(1-s^{\prime})^{2n+2}}\right]
520  \f]
521  This function is stored in \ref f_rho . Similar
522  definitions are made for \ref f_gamma and \ref f_omega .
523 
524  The Keplerial orbit at the equator is
525  \f[
526  \Omega_K = \frac{\omega^{\prime}}{2 \psi^{\prime}} ...
527  \f]
528  (eq. 31 in \ref Stergioulas03 )
529 
530  */
531  class nstar_rot {
532 
533  public:
534 
536  typedef boost::numeric::ublas::range ub_range;
537  typedef boost::numeric::ublas::vector_range
538  <boost::numeric::ublas::vector<double> > ubvector_range;
540 
541  /// The number of grid points in the \f$ \mu \f$ direction
542  int MDIV;
543  /// The number of grid points in the \f$ s \f$ direction
544  int SDIV;
545  /// The number of Legendre polynomials
546  int LMAX;
547 
548  /// Resize the grid
549  void resize(int MDIV_new, int SDIV_new, int LMAX_new,
550  int RDIV_new);
551 
552  /** \brief Default solver
553  */
555 
556  /** \brief Set new solver
557  */
559  mrootp=&m;
560  return 0;
561  }
562 
563  protected:
564 
565  /// Solver
567 
568  /// Solve for the Keplerian velocity
569  int solve_kepler(size_t nv, const ubvector &x, ubvector &y);
570 
571  /// Solve for the gravitational mass
572  int solve_grav_mass(size_t nv, const ubvector &x, ubvector &y,
573  double grav_mass);
574 
575  /// Solve for the gravitational mass
576  int solve_bar_mass(size_t nv, const ubvector &x, ubvector &y,
577  double bar_mass);
578 
579  /// Solve for the gravitational mass
580  int solve_ang_vel(size_t nv, const ubvector &x, ubvector &y,
581  double ang_vel);
582 
583  /// Solve for the gravitational mass
584  int solve_ang_mom(size_t nv, const ubvector &x, ubvector &y,
585  double ang_mom);
586 
587  /** \brief Subclass of \ref nstar_rot which specifies the function
588  to invert a polytropic EOS
589  */
591 
592  protected:
593 
594  /** \brief The polytropic index
595  */
596  double _Gamma_P;
597 
598  /** \brief The energy density
599  */
600  double _ee;
601 
602  public:
603 
604  /** \brief Create a function object with specified
605  polytropic index and ?
606  */
607  polytrope_solve(double Gamma_P, double ee) {
608  _Gamma_P=Gamma_P;
609  _ee=ee;
610  }
611 
612  /** \brief The function
613  */
614  double operator()(double rho0) {
615  return pow(rho0,_Gamma_P)/(_Gamma_P-1.0)+rho0-_ee;
616  }
617 
618  };
619 
620  /// The polytrope solver
622 
623  /// Array search object
625 
626  /** \brief The number of grid points in integration of TOV equations
627  for spherical stars
628  */
629  int RDIV;
630 
631  /** \brief Maximum value of s-coordinate (default 0.9999) */
632  double SMAX;
633  /** \brief Spacing in \f$ s \f$ direction,
634  \f$ \mathrm{SMAX}/(\mathrm{SDIV}-1) \f$
635  */
636  double DS;
637  /** \brief Spacing in \f$ \mu \f$ direction, \f$ 1/(\mathrm{MDIV}-1) \f$
638  */
639  double DM;
640 
641  /// Minimum radius for spherical stars (default \f$ 10^{-15} \f$)
642  double RMIN;
643 
644  /// \name Grid quantities set in make_grid()
645  //@{
646  /// \f$ s \f$
647  ubvector s_gp;
648  /// \f$ s (1-s) \f$
649  ubvector s_1_s;
650  /// \f$ 1-s \f$
651  ubvector one_s;
652  /// \f$ \mu \f$
653  ubvector mu;
654  /// \f$ 1-\mu^2 \f$
655  ubvector one_m2;
656  /// \f$ \theta \f$ defined by \f$ \mathrm{acos}~\mu \f$
657  ubvector theta;
658  /// \f$ \sin \theta \f$
659  ubvector sin_theta;
660  //@}
661 
662  /// \name Grid values computed in integrate() for spherical_star()
663  //@{
664  /// Isotropic radius
665  ubvector r_gp;
666  /// Radial coordinate
667  ubvector r_is_gp;
668  /// Metric function \f$ \lambda \f$
669  ubvector lambda_gp;
670  /// Metric function \f$ \nu \f$
671  ubvector nu_gp;
672  /// Enclosed gravitational mass
673  ubvector m_gp;
674  /// Energy density
675  ubvector e_d_gp;
676  //@}
677 
678  /// Desc
679  ubmatrix dgds;
680  /// Desc
681  ubmatrix dgdm;
682 
683  /// \name Metric functions
684  //@{
685  /** \brief potential \f$ \rho \f$ */
686  ubmatrix rho;
687  /** \brief potential \f$ \gamma \f$ */
688  ubmatrix gamma;
689  /** \brief potential \f$ \omega \f$ */
690  ubmatrix omega;
691  /** \brief potential \f$ \alpha \f$ */
692  ubmatrix alpha;
693  //@}
694 
695  /// \name Initial guess computed by the comp() function
696  //@{
697  /// Guess for the equatorial radius
698  double r_e_guess;
699  /** \brief Guess for \f$ \rho \f$ */
700  ubmatrix rho_guess;
701  /** \brief Guess for \f$ \gamma \f$ */
702  ubmatrix gamma_guess;
703  /** \brief Guess for \f$ \alpha \f$ */
704  ubmatrix omega_guess;
705  /** \brief Guess for \f$ \omega \f$ */
706  ubmatrix alpha_guess;
707  //@}
708 
709  /// \name EOS quantities
710  //@{
711  /** \brief Energy density \f$ \epsilon \f$ */
712  ubmatrix energy;
713  /** \brief Pressure */
714  ubmatrix pressure;
715  /** \brief Enthalpy */
716  ubmatrix enthalpy;
717  //@}
718 
719  /// \name Other quantities defined over the full two-dimensional grid
720  //@{
721  /** \brief Proper velocity squared */
722  ubmatrix velocity_sq;
723  /** \brief Derivative of \f$ \alpha \f$ with respect to \f$ \mu \f$ */
724  ubmatrix da_dm;
725  //@}
726 
727  /// \name Quantities defined for fixed values of mu
728  //@{
729  /** \brief \f$ \gamma(s) \f$ at \f$ \mu=1 \f$ */
730  ubvector gamma_mu_1;
731  /** \brief \f$ \gamma(s) \f$ at \f$ \mu=0 \f$ */
732  ubvector gamma_mu_0;
733  /** \brief \f$ \rho(s) \f$ at \f$ \mu=1 \f$ */
734  ubvector rho_mu_1;
735  /** \brief \f$ \rho(s) \f$ at \f$ \mu=0 \f$ */
736  ubvector rho_mu_0;
737  /** \brief \f$ \omega(s) \f$ at \f$ \mu=0 \f$ */
738  ubvector omega_mu_0;
739  //@}
740 
741  /** \brief The value of \f$ \hat{\gamma} \f$ at the pole */
742  double gamma_pole_h;
743  /** \brief The value of \f$ \hat{\gamma} \f$ at the center */
744  double gamma_center_h;
745  /** \brief The value of \f$ \hat{\gamma} \f$ at the equator */
747  /** \brief The value of \f$ \hat{\rho} \f$ at the pole */
748  double rho_pole_h;
749  /** \brief The value of \f$ \hat{\rho} \f$ at the center */
750  double rho_center_h;
751  /** \brief The value of \f$ \hat{\rho} \f$ at the equator */
752  double rho_equator_h;
753  /** \brief The value of \f$ \hat{\omega} \f$ at the equator */
755  /** \brief Angular velocity, \f$ \hat{\omega} \f$ */
756  double Omega_h;
757  /** \brief Central pressure */
758  double p_center;
759  /** \brief Central enthalpy */
760  double h_center;
761 
762  /// \name Desc
763  //@{
764  /** \brief \f$ f_{\rho}(s,n,s') \f$ */
766  /** \brief \f$ f_{\gamma}(s,n,s') \f$ */
768  /** \brief \f$ f_{\omega}(s,n,s') \f$ */
770  //@}
771 
772  /// \name Legendre polynomials
773  //@{
774  /** \brief Legendre polynomial \f$ P_{2n}(\mu) \f$
775  */
776  ubmatrix P_2n;
777  /** \brief Associated Legendre polynomial \f$ P^1_{2n-1}(\mu) \f$
778  */
779  ubmatrix P1_2n_1;
780  //@}
781 
782  /** \brief Integrated term over m in eqn for \f$ \rho \f$ */
783  ubmatrix D1_rho;
784  /** \brief Integrated term over m in eqn for \f$ \gamma \f$ */
785  ubmatrix D1_gamma;
786  /** \brief Integ. term over m in eqn for \f$ \omega \f$ */
787  ubmatrix D1_omega;
788  /** \brief Integrated term over s in eqn for \f$ \rho \f$ */
789  ubmatrix D2_rho;
790  /** \brief Integrated term over s in eqn for \f$ \gamma \f$ */
791  ubmatrix D2_gamma;
792  /** \brief Integ. term over s in eqn for \f$ \omega \f$ */
793  ubmatrix D2_omega;
794 
795  /** \brief source term in eqn for \f$ \gamma \f$ */
796  ubmatrix S_gamma;
797  /** \brief source term in eqn for \f$ \rho \f$ */
798  ubmatrix S_rho;
799  /** \brief source term in eqn for \f$ \omega \f$ */
800  ubmatrix S_omega;
801 
802  /** \brief The tolerance for the functions with the prefix "fix"
803  (default \f$ 10^{-4} \f$ )
804  */
805  double tol_abs;
806 
807  /// \name Thermodyanmic quantities near the surface
808  //@{
809  /// Pressure at the surface
810  double p_surface;
811  /// Energy density at the surface
812  double e_surface;
813  /** \brief Minimum specific enthalpy
814  */
815  double enthalpy_min;
816  //@}
817 
818  /// \name Polytrope parameters
819  //@{
820  /// Polytropic index
821  double n_P;
822  /// Polytropic exponent
823  double Gamma_P;
824  //@}
825 
826  /// \name Interpolation functions
827  //@{
828  /** \brief Cache for interpolation
829  */
831 
832  /// Search in array \c x of length \c n for value \c val
833  int new_search(int n, ubvector &x, double val);
834 
835  /** \brief Driver for the interpolation routine.
836 
837  First we find the tab. point nearest to xb, then we
838  interpolate using four points around xb.
839 
840  Used by \ref int_z(), \ref e_at_p(), \ref p_at_e(), \ref
841  p_at_h(), \ref h_at_p(), \ref n0_at_e(), \ref comp_omega(),
842  \ref comp_M_J(), \ref comp(), \ref spherical_star(), \ref
843  iterate().
844  */
845  double interp(ubvector &xp, ubvector &yp, int np, double xb);
846 
847  /** \brief Driver for the interpolation routine.
848 
849  Four point interpolation at a given offset the index of the
850  first point k.
851 
852  Used in \ref comp() .
853  */
854  double interp_4_k(ubvector &xp, ubvector &yp, int np, double xb, int k);
855  //@}
856 
857  /** \brief Integrate f[mu] from m-1 to m.
858 
859  This implements a 8-point closed Newton-Cotes formula.
860 
861  Used in \ref comp() .
862  */
863  double int_z(ubvector &f, int m);
864 
865  /// \name EOS functions
866  //@{
867  /** \brief Compute \f$ \varepsilon(P) \f$
868 
869  Used in \ref dm_dr_is(), \ref dp_dr_is(), \ref integrate()
870  and \ref iterate().
871  */
872  double e_at_p(double pp);
873 
874  /** \brief Compute \f$ P(\varepsilon) \f$
875 
876  Used in \ref make_center() and \ref integrate().
877  */
878  double p_at_e(double ee);
879 
880  /** \brief Pressure at fixed enthalpy
881 
882  Used in \ref iterate().
883  */
884  double p_at_h(double hh);
885 
886  /** \brief Enthalpy at fixed pressure
887 
888  Used in \ref make_center() and \ref integrate().
889  */
890  double h_at_p(double pp);
891 
892  /** \brief Baryon density at fixed energy density
893 
894  Used in \ref comp_M_J() and \ref comp() .
895  */
896  double n0_at_e(double ee);
897  //@}
898 
899  /// \name Derivatives on the grid
900  //@{
901  /** \brief Returns the derivative w.r.t. s of an array f[SDIV+1].
902  */
903  double s_deriv(ubvector &f, int s);
904 
905  /** \brief Returns the derivative w.r.t. mu of an array f[MDIV+1].
906  */
907  double m_deriv(ubvector &f, int m);
908 
909  /** \brief Returns the derivative w.r.t. s
910  */
911  double deriv_s(ubmatrix &f, int s, int m);
912 
913  /** \brief Returns the derivative w.r.t. mu
914  */
915  double deriv_m(ubmatrix &f, int s, int m);
916 
917  /** \brief Returns the derivative w.r.t. s and mu
918  */
919  double deriv_sm(ubmatrix &f, int s, int m);
920  //@}
921 
922  /// \name Initialization functions
923  //@{
924  /** \brief Returns the Legendre polynomial of degree n, evaluated at x.
925 
926  This uses the recurrence relation and is used in \ref comp_f_P()
927  which is called by the constructor.
928  */
929  double legendre(int n, double x);
930 
931  /** \brief Compute two-point functions
932 
933  This function computes the 2-point functions \f$
934  f^m_{2n}(r,r') \f$ used to integrate the potentials \f$ \rho,
935  \gamma \f$ and \f$ \omega \f$ (See \ref Komatsu89 for
936  details). Since the grid points are fixed, we can compute the
937  functions \ref f_rho, \ref f_gamma, \ref f_omega, \ref P_2n,
938  and \ref P1_2n_1 once at the beginning.
939 
940  See Eqs. 27-29 of \ref Cook92 and Eqs. 33-35 of \ref
941  Komatsu89. This function is called by the constructor.
942  */
943  void comp_f_P();
944 
945  /** \brief Create computational mesh.
946 
947  Create the computational mesh for \f$ s=r/(r+r_e) \f$
948  (where \f$ r_e \f$ is the coordinate equatorial radius)
949  and \f$ \mu = \cos \theta \f$
950  using
951  \f[
952  s[i]=\mathrm{SMAX}\left(\frac{i-1}{\mathrm{SDIV}-1}\right)
953  \f]
954  \f[
955  \mu[j]=\left(\frac{i-1}{\mathrm{MDIV}-1}\right)
956  \f]
957  When \f$ r=0 \f$, \f$ s=0 \f$, when \f$ r=r_e \f$,
958  \f$ s=1/2 \f$, and when \f$ r = \infty \f$, \f$ s=1 \f$ .
959  Inverting the relationship between \f$ r \f$ and \f$ s \f$
960  gives \f$ r = r_e s / (1-s) \f$ .
961  \comment
962  (Note that some versions of the manual have a typo,
963  giving \f$ 1-i \f$ rather than \f$ i-1 \f$ above.)
964  \endcomment
965 
966  Points in the mu-direction are stored in the array
967  <tt>mu[i]</tt>. Points in the s-direction are stored in the
968  array <tt>s_gp[j]</tt>.
969 
970  This function sets \ref s_gp, \ref s_1_s, \ref one_s,
971  \ref mu, \ref one_m2, \ref theta and \ref sin_theta .
972  All of these arrays are unit-indexed. It is called by
973  the constructor.
974  */
975  void make_grid();
976  //@}
977 
978  /** \brief Compute central pressure and enthalpy from central
979  energy density
980 
981  For polytropic EOSs, this also computes <tt>rho0_center</tt> .
982  */
983  void make_center(double e_center);
984 
985  /// \name Post-processing functions
986  //@{
987  /** \brief Compute Omega and Omega_K.
988  */
989  void comp_omega();
990 
991  /** \brief Compute rest mass and angular momentum.
992  */
993  void comp_M_J();
994 
995  /** \brief Compute various quantities.
996 
997  The main post-processing funciton
998  */
999  void comp();
1000  //@}
1001 
1002  /// \name For computing spherical stars
1003  //@{
1004  /** \brief Computes a spherically symmetric star
1005 
1006  The metric is
1007  \f[
1008  ds^2 = -e^{2\nu}dt^2 + e^{2\lambda} dr^2 + r^2 d\theta^2 +
1009  r^2 sin^2\theta d\phi^2
1010  \f]
1011  where \f$ r \f$ is an isotropic radial coordinate
1012  (corresponding to <tt>r_is</tt> in the code).
1013 
1014  This function computes \ref r_e_guess, \ref R_e,
1015  \ref Mass, and \ref Z_p .
1016  */
1017  void spherical_star();
1018 
1019  /** \brief Derivative of gravitational mass with respect to
1020  isotropic radius */
1021  double dm_dr_is(double r_is, double r, double m, double p);
1022 
1023  /** \brief Derivative of pressure with respect to isotropic radius */
1024  double dp_dr_is(double r_is, double r, double m, double p);
1025 
1026  /** \brief Derivative of radius with respect to isotropic radius */
1027  double dr_dr_is(double r_is, double r, double m);
1028 
1029  /** \brief Integrate one of the differential equations for
1030  spherical stars*/
1031  void integrate(int i_check, double &r_final, double &m_final,
1032  double &r_is_final);
1033  //@}
1034 
1035  /** \brief Main iteration function
1036  */
1037  int iterate(double r_ratio, double tol_rel);
1038 
1039  /// \name EOS member variables
1040  //@{
1041  /** \brief If true, then an EOS has been set
1042  */
1043  bool eos_set;
1044 
1045  /** \brief If true, then use a polytrope and rescale
1046  */
1048 
1049  /** \brief Pointer to the user-specified EOS
1050  */
1052  //@}
1053 
1054  /** \brief Compute masses and angular momentum
1055  */
1056  void calc_masses_J(ubmatrix &rho_0);
1057 
1058  public:
1059 
1060  nstar_rot();
1061 
1062  /** \brief Relative accuracy for the equatorial radius,
1063  \f$ r_e \f$ (default \f$ 10^{-5} \f$)
1064 
1065  Used in \ref iterate() .
1066  */
1068 
1069  /** \brief Accuracy for equatorial radius using alternate
1070  solvers (default \f$ 10^{-9} \f$)
1071  */
1072  double alt_tol_rel;
1073 
1074  /** \brief Verbosity parameter
1075  */
1076  int verbose;
1077 
1078  /** \brief Create an output table
1079  */
1080  void output_table(o2scl::table3d &t);
1081 
1082  /// \name Output
1083  //@{
1084  /** \brief Central energy density (in units of
1085  \f$ 10^{15} \mathrm{g}/\mathrm{cm}^3 \f$)
1086  */
1087  double e_center;
1088  /** \brief Ratio of polar to equatorial radius
1089  */
1090  double r_ratio;
1091  /** \brief Coordinate equatorial radius
1092  */
1093  double r_e;
1094  //@}
1095 
1096  /// \name Quantities computed by nstar_rot::comp() (in order)
1097  //@{
1098  /** \brief Radius at pole */
1099  double r_p;
1100  /** \brief The value of the s-coordinate at the pole */
1101  double s_p;
1102  /** \brief The value of the s-coordinate at the equator */
1103  double s_e;
1104  /// The velocity at the equator
1106  /** \brief Circumferential radius in cm (i.e. the radius defined
1107  such that \f$ 2 \pi R_e \f$ is the proper circumference) */
1108  double R_e;
1109  /// Proper mass (in \f$ \mathrm{g} \f$ )
1110  double Mass_p;
1111  /// Gravitational mass (in \f$ \mathrm{g} \f$ )
1112  double Mass;
1113  /// Baryonic mass (in \f$ \mathrm{g} \f$ )
1114  double Mass_0;
1115  /// Angular momentum (in \f$ \mathrm{g}~\mathrm{cm}^2 / \mathrm{s} \f$ )
1116  double J;
1117  /// Angular velocity (in \f$ \mathrm{radians}/\mathrm{s} \f$ )
1118  double Omega;
1119  /// Total rotational kinetic energy
1120  double T;
1121  /// Moment of inertia
1122  double I;
1123  /// Gravitational binding energy
1124  double W;
1125  /// Polar redshift
1126  double Z_p;
1127  /// Forward equatorial redshift
1128  double Z_f;
1129  /// Backward equatorial redshift
1130  double Z_b;
1131  /** \brief Kepler rotation frequency (in 1/s) */
1132  double Omega_K;
1133  /// The eccentricity
1135  /// Desc
1136  ubvector v_plus;
1137  /// Desc
1138  ubvector v_minus;
1139  /// Desc
1140  double vel_plus;
1141  /// Desc
1142  double vel_minus;
1143  /** \brief Height from surface of last stable co-rotating circular
1144  orbit in equatorial plane
1145 
1146  If this is zero then all orbits are stable.
1147  */
1148  double h_plus;
1149  /** \brief Height from surface of last stable counter-rotating circular
1150  orbit in equatorial plane
1151 
1152  If this is zero then all orbits are stable.
1153  */
1154  double h_minus;
1155  /// Desc
1156  double Omega_plus;
1157  /// Desc
1158  double u_phi;
1159  /// Angular velocity of a particle in a circular orbit at the equator
1160  double Omega_p;
1161  /// Desc
1162  double grv2;
1163  /// Desc
1164  double grv2_new;
1165  /// Desc
1166  double grv3;
1167  /** \brief Ratio of potential \f$ \omega \f$ to angular
1168  velocity \f$ \Omega \f$
1169  */
1170  double om_over_Om;
1171  /** \brief Mass quadrupole moment
1172  */
1174  //@}
1175 
1176  /// \name Settings
1177  //@{
1178  /// The convergence factor (default 1.0)
1179  double cf;
1180  //@}
1181 
1182  /// \name Internal constants
1183  //@{
1184  /** \brief Use the values of the constants from the original RNS
1185  code
1186  */
1187  void constants_rns();
1188  /** \brief Use the \o2 values
1189  */
1190  void constants_o2scl();
1191  /** \brief Speed of light in vacuum (in CGS units) */
1192  double C;
1193  /** \brief Gravitational constant (in CGS units) */
1194  double G;
1195  /** \brief Mass of sun (in g) */
1196  double MSUN;
1197  /** \brief Square of length scale in CGS units,
1198  \f$ \kappa \equiv 10^{-15} c^2/G \f$
1199  */
1200  double KAPPA;
1201  /** \brief The mass of one baryon (in g)
1202  */
1203  double MB;
1204  /** \brief The value \f$ \kappa G c^{-4} \f$ */
1205  double KSCALE;
1206  /// The constant \f$ \pi \f$
1207  double PI;
1208  //@}
1209 
1210  /// \name Basic Usage
1211  //@{
1212  /** \brief Set the EOS
1213  */
1214  void set_eos(eos_nstar_rot &eos) {
1215  eosp=&eos;
1216  eos_set=true;
1217  scaled_polytrope=false;
1218  return;
1219  }
1220 
1221  /** \brief Use a polytropic EOS with a specified index
1222  */
1223  void polytrope_eos(double index) {
1224  n_P=index;
1225  scaled_polytrope=true;
1226  eos_set=true;
1227  return;
1228  }
1229 
1230  /** \brief Construct a configuration with a fixed central
1231  energy density and a fixed axis ratio
1232 
1233  The central energy density should be in \f$
1234  \mathrm{g}/\mathrm{cm}^3 \f$ and the axis ratio is unitless.
1235  This is fastest of the high-level interface functions as it
1236  doesn't require an additional solver.
1237  */
1238  int fix_cent_eden_axis_rat(double cent_eden, double axis_rat,
1239  bool use_guess=false);
1240 
1241  /** \brief Construct a configuration with a fixed central
1242  energy density and a fixed gravitational mass
1243 
1244  The central energy density should be in \f$
1245  \mathrm{g}/\mathrm{cm}^3 \f$ and the gravitational
1246  mass should be in solar masses.
1247  */
1248  int fix_cent_eden_grav_mass(double cent_eden, double grav_mass);
1249 
1250  /** \brief Construct a configuration with a fixed central
1251  energy density and a fixed baryonic mass
1252 
1253  The central energy density should be in \f$
1254  \mathrm{g}/\mathrm{cm}^3 \f$ and the baryonic
1255  mass should be in solar masses.
1256  */
1257  int fix_cent_eden_bar_mass(double cent_eden, double bar_mass);
1258 
1259  /** \brief Construct a configuration with a fixed central
1260  energy density and the Keplerian rotation rate
1261 
1262  The central energy density should be in \f$
1263  \mathrm{g}/\mathrm{cm}^3 \f$ .
1264  */
1265  int fix_cent_eden_with_kepler(double cent_eden);
1266 
1267  /** \brief Experimental alternate form for
1268  \ref fix_cent_eden_with_kepler()
1269  */
1270  int fix_cent_eden_with_kepler_alt(double cent_eden,
1271  bool use_guess=false);
1272 
1273  /** \brief Experimental alternate form for
1274  \ref fix_cent_eden_grav_mass()
1275  */
1276  int fix_cent_eden_grav_mass_alt(double cent_eden, double grav_mass,
1277  bool use_guess=false);
1278 
1279  /** \brief Experimental alternate form for
1280  \ref fix_cent_eden_bar_mass()
1281  */
1282  int fix_cent_eden_bar_mass_alt(double cent_eden, double bar_mass,
1283  bool use_guess=false);
1284 
1285  /** \brief Experimental alternate form for
1286  \ref fix_cent_eden_ang_vel()
1287  */
1288  int fix_cent_eden_ang_vel_alt(double cent_eden, double ang_vel,
1289  bool use_guess=false);
1290 
1291  /** \brief Experimental alternate form for
1292  \ref fix_cent_eden_ang_mom()
1293  */
1294  int fix_cent_eden_ang_mom_alt(double cent_eden, double ang_mom,
1295  bool use_guess=false);
1296 
1297  /** \brief Construct a non-rotating configuration with a fixed central
1298  energy density
1299 
1300  The central energy density should be in \f$
1301  \mathrm{g}/\mathrm{cm}^3 \f$ .
1302  */
1303  int fix_cent_eden_non_rot(double cent_eden);
1304 
1305  /** \brief Construct a configuration with a fixed central
1306  energy density and a fixed angular velocity.
1307 
1308  The central energy density should be in \f$
1309  \mathrm{g}/\mathrm{cm}^3 \f$ and the angular
1310  velocity should be in \f$ \mathrm{rad}/\mathrm{s} \f$.
1311  The final angular velocity (possibly slightly different
1312  than <tt>ang_vel</tt> is stored in \ref Omega .
1313 
1314  \note In the original RNS code, the <tt>ang_vel</tt> argument
1315  is different because it was rescaled by a factor of \f$ 10^{4}
1316  \f$.
1317  */
1318  int fix_cent_eden_ang_vel(double cent_eden, double ang_vel);
1319 
1320  /** \brief Construct a configuration with a fixed central
1321  energy density and a fixed angular momentum.
1322 
1323  The central energy density should be in \f$
1324  \mathrm{g}/\mathrm{cm}^3 \f$. The angular momentum should be
1325  in units of \f$ G M_{\odot}^2/C \f$ .
1326  */
1327  int fix_cent_eden_ang_mom(double cent_eden, double ang_mom);
1328  //@}
1329 
1330  /** \name Testing functions
1331 
1332  All these compare with hard-coded results obtained with
1333  the RNS code.
1334  */
1335  //@{
1336  /** \brief Test determining configuration with fixed central
1337  energy density and fixed radius ratio with EOS C
1338  */
1339  void test1(o2scl::test_mgr &t);
1340 
1341  /** \brief Test configuration rotating and Keplerian frequency
1342  with a fixed central energy density and EOS C
1343  */
1344  void test2(o2scl::test_mgr &t);
1345 
1346  /** \brief Test fixed central energy density and fixed
1347  gravitational mass with EOS C
1348  */
1349  void test3(o2scl::test_mgr &t);
1350 
1351  /** \brief Test fixed central energy density and fixed baryonic
1352  mass with EOS C
1353  */
1354  void test4(o2scl::test_mgr &t);
1355 
1356  /** \brief Test fixed central energy density and fixed angular
1357  velocity with EOS C
1358  */
1359  void test5(o2scl::test_mgr &t);
1360 
1361  /** \brief Test fixed central energy density and fixed angular
1362  momentum with EOS C
1363  */
1364  void test6(o2scl::test_mgr &t);
1365 
1366  /** \brief Test a series of non-rotating stars on a energy density
1367  grid with EOS C
1368  */
1369  void test7(o2scl::test_mgr &t);
1370 
1371  /** \brief Test Keplerian frequency for a polytrope
1372  */
1373  void test8(o2scl::test_mgr &t);
1374  //@}
1375 
1376 
1377  };
1378 
1379 }
1380 
1381 #endif
ubvector rho_mu_0
at
Definition: nstar_rot.h:736
ubvector theta
defined by
Definition: nstar_rot.h:657
double rho_center_h
The value of at the center.
Definition: nstar_rot.h:750
virtual double pr_from_enth(double enth)=0
From the enthalpy, return the pressure.
Rotating neutron star class based on RNS v1.1d from N. Stergioulas et al.
Definition: nstar_rot.h:531
double h_minus
Height from surface of last stable counter-rotating circular orbit in equatorial plane.
Definition: nstar_rot.h:1154
ubvector m_gp
Enclosed gravitational mass.
Definition: nstar_rot.h:673
ubmatrix gamma
potential
Definition: nstar_rot.h:688
ubmatrix P1_2n_1
Associated Legendre polynomial .
Definition: nstar_rot.h:779
ubmatrix dgdm
Desc.
Definition: nstar_rot.h:681
double KSCALE
The value .
Definition: nstar_rot.h:1205
double MB
The mass of one baryon (in g)
Definition: nstar_rot.h:1203
double p_center
Central pressure.
Definition: nstar_rot.h:758
double Mass_p
Proper mass (in )
Definition: nstar_rot.h:1110
virtual double nb_from_ed(double ed)=0
From the energy density, return the baryon density.
double vel_plus
Desc.
Definition: nstar_rot.h:1140
ubmatrix pressure
Pressure.
Definition: nstar_rot.h:714
bool scaled_polytrope
If true, then use a polytrope and rescale.
Definition: nstar_rot.h:1047
double Z_f
Forward equatorial redshift.
Definition: nstar_rot.h:1128
double gamma_center_h
The value of at the center.
Definition: nstar_rot.h:744
o2scl::mroot * mrootp
Solver.
Definition: nstar_rot.h:566
double gamma_equator_h
The value of at the equator.
Definition: nstar_rot.h:746
ubvector s_gp
Definition: nstar_rot.h:647
virtual double nb_from_pr(double pr)=0
From the pressure, return the baryon density.
double R_e
Circumferential radius in cm (i.e. the radius defined such that is the proper circumference) ...
Definition: nstar_rot.h:1108
virtual void ed_nb_from_pr(double pr, double &ed, double &nb)=0
Given the pressure, produce the energy and number densities.
ubmatrix rho_guess
Guess for .
Definition: nstar_rot.h:700
void polytrope_eos(double index)
Use a polytropic EOS with a specified index.
Definition: nstar_rot.h:1223
double Gamma_P
Polytropic exponent.
Definition: nstar_rot.h:823
double mass_quadrupole
Mass quadrupole moment.
Definition: nstar_rot.h:1173
lib_settings_class o2scl_settings
ubvector lambda_gp
Metric function .
Definition: nstar_rot.h:669
virtual double ed_from_nb(double nb)=0
From the baryon density, return the energy density.
double tol_abs
The tolerance for the functions with the prefix "fix" (default )
Definition: nstar_rot.h:805
ubmatrix S_gamma
source term in eqn for
Definition: nstar_rot.h:796
ubvector e_d_gp
Energy density.
Definition: nstar_rot.h:675
int RDIV
The number of grid points in integration of TOV equations for spherical stars.
Definition: nstar_rot.h:629
ubvector nu_gp
Metric function .
Definition: nstar_rot.h:671
ubvector v_minus
Desc.
Definition: nstar_rot.h:1138
double cf
The convergence factor (default 1.0)
Definition: nstar_rot.h:1179
double eccentricity
The eccentricity.
Definition: nstar_rot.h:1134
double Omega_p
Angular velocity of a particle in a circular orbit at the equator.
Definition: nstar_rot.h:1160
ubmatrix D2_omega
Integ. term over s in eqn for .
Definition: nstar_rot.h:793
double Omega
Angular velocity (in )
Definition: nstar_rot.h:1118
void output()
Output EOS table to screen.
Definition: nstar_rot.h:375
o2scl::root_bkt_cern< polytrope_solve > rbc
The polytrope solver.
Definition: nstar_rot.h:621
double C
Speed of light in vacuum (in CGS units)
Definition: nstar_rot.h:131
double KAPPA
Square of length scale in CGS units, .
Definition: nstar_rot.h:137
double rho_equator_h
The value of at the equator.
Definition: nstar_rot.h:752
double J
Angular momentum (in )
Definition: nstar_rot.h:1116
double DS
Spacing in direction, .
Definition: nstar_rot.h:636
double grv2
Desc.
Definition: nstar_rot.h:1162
virtual double pr_from_ed(double ed)=0
From the energy density, return the pressure.
ubmatrix D1_rho
Integrated term over m in eqn for .
Definition: nstar_rot.h:783
int n_nearest
Cache for interpolation.
Definition: nstar_rot.h:144
double h_center
Central enthalpy.
Definition: nstar_rot.h:760
int set_solver(o2scl::mroot<> &m)
Set new solver.
Definition: nstar_rot.h:558
double T
Total rotational kinetic energy.
Definition: nstar_rot.h:1120
ubmatrix velocity_sq
Proper velocity squared.
Definition: nstar_rot.h:722
double G
Gravitational constant (in CGS units)
Definition: nstar_rot.h:133
ubvector s_1_s
Definition: nstar_rot.h:649
ubmatrix S_rho
source term in eqn for
Definition: nstar_rot.h:798
double Z_b
Backward equatorial redshift.
Definition: nstar_rot.h:1130
ubmatrix omega_guess
Guess for .
Definition: nstar_rot.h:704
ubmatrix omega
potential
Definition: nstar_rot.h:690
double I
Moment of inertia.
Definition: nstar_rot.h:1122
ubvector r_gp
Isotropic radius.
Definition: nstar_rot.h:665
void set_eos_native(vec1_t &eden, vec2_t &pres, vec3_t &enth, vec4_t &nb)
Set the EOS from four vectors in the native unit system.
Definition: nstar_rot.h:163
polytrope_solve(double Gamma_P, double ee)
Create a function object with specified polytropic index and ?
Definition: nstar_rot.h:607
double r_e_guess
Guess for the equatorial radius.
Definition: nstar_rot.h:698
virtual double enth_from_pr(double pr)=0
From the pressure, return the enthalpy.
double W
Gravitational binding energy.
Definition: nstar_rot.h:1124
double r_p
Radius at pole.
Definition: nstar_rot.h:1099
double omega_equator_h
The value of at the equator.
Definition: nstar_rot.h:754
ubvector mu
Definition: nstar_rot.h:653
Tabulated EOS for nstar_rot from Bethe74.
Definition: nstar_rot.h:388
double vel_minus
Desc.
Definition: nstar_rot.h:1142
int n_tab
number of tabulated EOS points
Definition: nstar_rot.h:118
double grv3
Desc.
Definition: nstar_rot.h:1166
double h_plus
Height from surface of last stable co-rotating circular orbit in equatorial plane.
Definition: nstar_rot.h:1148
Tabulated EOS for nstar_rot from Pandharipande75.
Definition: nstar_rot.h:395
const double speed_of_light
ubmatrix da_dm
Derivative of with respect to .
Definition: nstar_rot.h:724
double rho_pole_h
The value of at the pole.
Definition: nstar_rot.h:748
Create a tabulated EOS for nstar_rot using interpolation.
Definition: nstar_rot.h:107
#define O2SCL_ERR2(d, d2, n)
double alt_tol_rel
Accuracy for equatorial radius using alternate solvers (default )
Definition: nstar_rot.h:1072
ubvector v_plus
Desc.
Definition: nstar_rot.h:1136
ubmatrix alpha
potential
Definition: nstar_rot.h:692
double MSUN
Mass of sun (in g)
Definition: nstar_rot.h:1196
double KSCALE
The value .
Definition: nstar_rot.h:139
A EOS base class for the TOV solver.
Definition: eos_tov.h:47
int SDIV
The number of grid points in the direction.
Definition: nstar_rot.h:544
double _ee
The energy density.
Definition: nstar_rot.h:600
void set_eos(eos_nstar_rot &eos)
Set the EOS.
Definition: nstar_rot.h:1214
ubvector omega_mu_0
at
Definition: nstar_rot.h:738
ubmatrix enthalpy
Enthalpy.
Definition: nstar_rot.h:716
tensor3 f_omega
Definition: nstar_rot.h:769
double enthalpy_min
Minimum specific enthalpy.
Definition: nstar_rot.h:815
double Omega_h
Angular velocity, .
Definition: nstar_rot.h:756
ubvector gamma_mu_0
at
Definition: nstar_rot.h:732
int verbose
Verbosity parameter.
Definition: nstar_rot.h:1076
double e_center
Central energy density (in units of )
Definition: nstar_rot.h:1087
double s_p
The value of the s-coordinate at the pole.
Definition: nstar_rot.h:1101
double _Gamma_P
The polytropic index.
Definition: nstar_rot.h:596
An EOS for nstar_rot.
Definition: nstar_rot.h:83
ubvector one_s
Definition: nstar_rot.h:651
tensor3 f_rho
Definition: nstar_rot.h:765
double p_surface
Pressure at the surface.
Definition: nstar_rot.h:810
int n_nearest
Cache for interpolation.
Definition: nstar_rot.h:830
double Z_p
Polar redshift.
Definition: nstar_rot.h:1126
ubmatrix energy
Energy density .
Definition: nstar_rot.h:712
double PI
The constant .
Definition: nstar_rot.h:1207
double grv2_new
Desc.
Definition: nstar_rot.h:1164
tensor3 f_gamma
Definition: nstar_rot.h:767
virtual double enth_from_nb(double nb)=0
From the baryon density, return the enthalpy.
virtual double ed_from_pr(double pr)=0
From the pressure, return the energy density.
ubmatrix D1_omega
Integ. term over m in eqn for .
Definition: nstar_rot.h:787
Subclass of nstar_rot which specifies the function to invert a polytropic EOS.
Definition: nstar_rot.h:590
ubmatrix rho
potential
Definition: nstar_rot.h:686
o2scl::search_vec< double * > sv
Array search object.
Definition: nstar_rot.h:112
double C
Speed of light in vacuum (in CGS units)
Definition: nstar_rot.h:1192
o2scl::mroot_hybrids def_mroot
Default solver.
Definition: nstar_rot.h:554
double eq_radius_tol_rel
Relative accuracy for the equatorial radius, (default )
Definition: nstar_rot.h:1067
ubmatrix gamma_guess
Guess for .
Definition: nstar_rot.h:702
double SMAX
Maximum value of s-coordinate (default 0.9999)
Definition: nstar_rot.h:632
ubmatrix D2_rho
Integrated term over s in eqn for .
Definition: nstar_rot.h:789
bool eos_set
If true, then an EOS has been set.
Definition: nstar_rot.h:1043
void set_eos_fm(size_t n, vec1_t &eden, vec2_t &pres, vec3_t &nb)
Set the EOS from energy density, pressure, and baryon density stored in powers of ...
Definition: nstar_rot.h:187
o2scl::search_vec< ubvector_range > sv_ub
Array search object.
Definition: nstar_rot.h:624
int MDIV
The number of grid points in the direction.
Definition: nstar_rot.h:542
double G
Gravitational constant (in CGS units)
Definition: nstar_rot.h:1194
ubvector one_m2
Definition: nstar_rot.h:655
double s_e
The value of the s-coordinate at the equator.
Definition: nstar_rot.h:1103
double gamma_pole_h
The value of at the pole.
Definition: nstar_rot.h:742
double Mass
Gravitational mass (in )
Definition: nstar_rot.h:1112
ubvector sin_theta
Definition: nstar_rot.h:659
ubmatrix alpha_guess
Guess for .
Definition: nstar_rot.h:706
double operator()(double rho0)
The function.
Definition: nstar_rot.h:614
int LMAX
The number of Legendre polynomials.
Definition: nstar_rot.h:546
eos_nstar_rot * eosp
Pointer to the user-specified EOS.
Definition: nstar_rot.h:1051
double r_ratio
Ratio of polar to equatorial radius.
Definition: nstar_rot.h:1090
double Omega_plus
Desc.
Definition: nstar_rot.h:1156
ubmatrix S_omega
source term in eqn for
Definition: nstar_rot.h:800
double velocity_equator
The velocity at the equator.
Definition: nstar_rot.h:1105
double Mass_0
Baryonic mass (in )
Definition: nstar_rot.h:1114
ubmatrix D2_gamma
Integrated term over s in eqn for .
Definition: nstar_rot.h:791
double n_P
Polytropic index.
Definition: nstar_rot.h:821
double e_surface
Energy density at the surface.
Definition: nstar_rot.h:812
ubmatrix D1_gamma
Integrated term over m in eqn for .
Definition: nstar_rot.h:785
double r_e
Coordinate equatorial radius.
Definition: nstar_rot.h:1093
ubmatrix P_2n
Legendre polynomial .
Definition: nstar_rot.h:776
ubvector gamma_mu_1
at
Definition: nstar_rot.h:730
ubmatrix dgds
Desc.
Definition: nstar_rot.h:679
virtual double pr_from_nb(double nb)=0
From the baryon density, return the pressure.
double DM
Spacing in direction, .
Definition: nstar_rot.h:639
double Omega_K
Kepler rotation frequency (in 1/s)
Definition: nstar_rot.h:1132
double KAPPA
Square of length scale in CGS units, .
Definition: nstar_rot.h:1200
const double gravitational_constant
ubvector rho_mu_1
at
Definition: nstar_rot.h:734
double RMIN
Minimum radius for spherical stars (default )
Definition: nstar_rot.h:642
ubvector r_is_gp
Radial coordinate.
Definition: nstar_rot.h:667
double u_phi
Desc.
Definition: nstar_rot.h:1158
double om_over_Om
Ratio of potential to angular velocity .
Definition: nstar_rot.h:1170

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