Go to the documentation of this file.
65 #include <boost/numeric/ublas/vector.hpp>
66 #include <boost/numeric/ublas/matrix.hpp>
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>
154 double interp(
double xp[],
double yp[],
int np,
double xb);
162 template<
class vec1_t,
class vec2_t,
class vec3_t,
class vec4_t>
173 for(
int i=1;i<=
n_tab;i++) {
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) {
189 static const int n_crust=78;
198 (
"1/fm^4",
"g/cm^3",1.0);
201 (
"1/fm^4",
"dyne/cm^2",1.0);
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}};
294 for(
size_t i=0;i<n_crust;i++) {
297 double mu=(nst_arr[i][0]/conv1+nst_arr[i][1]/conv2)/
298 nst_arr[i][2]*1.0e39;
311 for(
size_t i=0;i<n;i++) {
314 log_h_tab[i+n_crust+1]=log10(log((eden[i]+pres[i])/nb[i]/
366 virtual void ed_nb_from_pr(
double pr,
double &ed,
double &nb);
376 for(
int i=
n_tab;i>=1;i--) {
380 std::cout << std::endl;
536 typedef boost::numeric::ublas::range ub_range;
537 typedef boost::numeric::ublas::vector_range
549 void resize(
int MDIV_new,
int SDIV_new,
int LMAX_new,
1021 double dm_dr_is(
double r_is,
double r,
double m,
double p);
1024 double dp_dr_is(
double r_is,
double r,
double m,
double p);
1027 double dr_dr_is(
double r_is,
double r,
double m);
1031 void integrate(
int i_check,
double &r_final,
double &m_final,
1032 double &r_is_final);
1239 bool use_guess=
false);
1271 bool use_guess=
false);
1277 bool use_guess=
false);
1283 bool use_guess=
false);
1289 bool use_guess=
false);
1295 bool use_guess=
false);
ubmatrix D2_rho
Integrated term over s in eqn for .
double e_surface
Energy density at the surface.
virtual double nb_from_pr(double pr)
From the pressure, return the baryon density.
double KAPPA
Square of length scale in CGS units, .
double n_P
Polytropic index.
double h_plus
Height from surface of last stable co-rotating circular orbit in equatorial plane.
int SDIV
The number of grid points in the direction.
void resize(int MDIV_new, int SDIV_new, int LMAX_new, int RDIV_new)
Resize the grid.
double interp_4_k(ubvector &xp, ubvector &yp, int np, double xb, int k)
Driver for the interpolation routine.
double KAPPA
Square of length scale in CGS units, .
double int_z(ubvector &f, int m)
Integrate f[mu] from m-1 to m.
double s_deriv(ubvector &f, int s)
Returns the derivative w.r.t. s of an array f[SDIV+1].
ubmatrix D1_omega
Integ. term over m in eqn for .
double r_e_guess
Guess for the equatorial radius.
o2scl::root_bkt_cern< polytrope_solve > rbc
The polytrope solver.
int fix_cent_eden_axis_rat(double cent_eden, double axis_rat, bool use_guess=false)
Construct a configuration with a fixed central energy density and a fixed axis ratio.
void constants_rns()
Use the values of the constants from the original RNS code.
double rho_pole_h
The value of at the pole.
double mass_quadrupole
Mass quadrupole moment.
double log_h_tab[201]
h points in EOS file
ubmatrix D1_rho
Integrated term over m in eqn for .
virtual double enth_from_pr(double pr)
From the pressure, return the enthalpy.
void comp_omega()
Compute Omega and Omega_K.
int solve_kepler(size_t nv, const ubvector &x, ubvector &y)
Solve for the Keplerian velocity.
ubvector m_gp
Enclosed gravitational mass.
virtual double enth_from_nb(double nb)=0
From the baryon density, return the enthalpy.
int solve_grav_mass(size_t nv, const ubvector &x, ubvector &y, double grav_mass)
Solve for the gravitational mass.
int set_solver(o2scl::mroot<> &m)
Set new solver.
double G
Gravitational constant (in CGS units)
ubmatrix rho_guess
Guess for .
virtual double pr_from_enth(double enth)
From the enthalpy, return the pressure.
double log_p_tab[201]
p points in tabulated EOS
double r_p
Radius at pole.
double _Gamma_P
The polytropic index.
double J
Angular momentum (in )
double C
Speed of light in vacuum (in CGS units)
#define O2SCL_ERR2(d, d2, n)
double p_surface
Pressure at the surface.
double DS
Spacing in direction, .
ubmatrix D2_gamma
Integrated term over s in eqn for .
int new_search(int n, ubvector &x, double val)
Search in array x of length n for value val.
void comp_M_J()
Compute rest mass and angular momentum.
Tabulated EOS for nstar_rot from Bethe74.
double Z_f
Forward equatorial redshift.
int solve_ang_vel(size_t nv, const ubvector &x, ubvector &y, double ang_vel)
Solve for the gravitational mass.
int RDIV
The number of grid points in integration of TOV equations for spherical stars.
double p_at_h(double hh)
Pressure at fixed enthalpy.
void integrate(int i_check, double &r_final, double &m_final, double &r_is_final)
Integrate one of the differential equations for spherical stars.
double I
Moment of inertia.
int fix_cent_eden_ang_vel(double cent_eden, double ang_vel)
Construct a configuration with a fixed central energy density and a fixed angular velocity.
Tabulated EOS for nstar_rot from Pandharipande75.
virtual double pr_from_ed(double ed)
From the energy density, return the pressure.
double velocity_equator
The velocity at the equator.
o2scl::search_vec< ubvector_range > sv_ub
Array search object.
double eq_radius_tol_rel
Relative accuracy for the equatorial radius, (default )
double Mass_0
Baryonic mass (in )
o2scl::mroot_hybrids def_mroot
Default solver.
virtual double enth_from_pr(double pr)=0
From the pressure, return the enthalpy.
bool eos_set
If true, then an EOS has been set.
double s_e
The value of the s-coordinate at the equator.
int fix_cent_eden_grav_mass(double cent_eden, double grav_mass)
Construct a configuration with a fixed central energy density and a fixed gravitational mass.
void output()
Output EOS table to screen.
double e_center
Central energy density (in units of )
double dp_dr_is(double r_is, double r, double m, double p)
Derivative of pressure with respect to isotropic radius.
bool scaled_polytrope
If true, then use a polytrope and rescale.
double r_e
Coordinate equatorial radius.
void test4(o2scl::test_mgr &t)
Test fixed central energy density and fixed baryonic mass with EOS C.
int iterate(double r_ratio, double tol_rel)
Main iteration function.
double _ee
The energy density.
virtual double enth_from_nb(double nb)
From the baryon density, return the enthalpy.
int n_tab
number of tabulated EOS points
int fix_cent_eden_non_rot(double cent_eden)
Construct a non-rotating configuration with a fixed central energy density.
double gamma_center_h
The value of at the center.
double MSUN
Mass of sun (in g)
double Mass_p
Proper mass (in )
void make_center(double e_center)
Compute central pressure and enthalpy from central energy density.
double Omega_K
Kepler rotation frequency (in 1/s)
double gamma_pole_h
The value of at the pole.
virtual double nb_from_ed(double ed)
From the energy density, return the baryon density.
void set_eos(eos_nstar_rot &eos)
Set the EOS.
ubmatrix S_rho
source term in eqn for
double log_n0_tab[201]
number density in EOS file
double m_deriv(ubvector &f, int m)
Returns the derivative w.r.t. mu of an array f[MDIV+1].
double enthalpy_min
Minimum specific enthalpy.
double Omega_h
Angular velocity, .
ubmatrix S_omega
source term in eqn for
void test1(o2scl::test_mgr &t)
Test determining configuration with fixed central energy density and fixed radius ratio with EOS C.
double h_at_p(double pp)
Enthalpy at fixed pressure.
ubmatrix P_2n
Legendre polynomial .
ubmatrix D1_gamma
Integrated term over m in eqn for .
virtual void ed_nb_from_pr(double pr, double &ed, double &nb)
Given the pressure, produce the energy and number densities.
virtual double pr_from_nb(double nb)
From the baryon density, return the pressure.
lib_settings_class o2scl_settings
double alt_tol_rel
Accuracy for equatorial radius using alternate solvers (default )
void test7(o2scl::test_mgr &t)
Test a series of non-rotating stars on a energy density grid with EOS C.
double Z_p
Polar redshift.
double interp(ubvector &xp, ubvector &yp, int np, double xb)
Driver for the interpolation routine.
ubvector r_gp
Isotropic radius.
ubmatrix alpha_guess
Guess for .
void polytrope_eos(double index)
Use a polytropic EOS with a specified index.
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.
void make_grid()
Create computational mesh.
double gamma_equator_h
The value of at the equator.
ubvector lambda_gp
Metric function .
double p_at_e(double ee)
Compute
double Gamma_P
Polytropic exponent.
double dr_dr_is(double r_is, double r, double m)
Derivative of radius with respect to isotropic radius.
double cf
The convergence factor (default 1.0)
void spherical_star()
Computes a spherically symmetric star.
int solve_bar_mass(size_t nv, const ubvector &x, ubvector &y, double bar_mass)
Solve for the gravitational mass.
double SMAX
Maximum value of s-coordinate (default 0.9999)
ubmatrix P1_2n_1
Associated Legendre polynomial .
polytrope_solve(double Gamma_P, double ee)
Create a function object with specified polytropic index and ?
ubmatrix velocity_sq
Proper velocity squared.
double operator()(double rho0)
The function.
void test8(o2scl::test_mgr &t)
Test Keplerian frequency for a polytrope.
void test5(o2scl::test_mgr &t)
Test fixed central energy density and fixed angular velocity with EOS C.
ubmatrix enthalpy
Enthalpy.
double R_e
Circumferential radius in cm (i.e. the radius defined such that is the proper circumference)
int verbose
Verbosity parameter.
double h_minus
Height from surface of last stable counter-rotating circular orbit in equatorial plane.
double C
Speed of light in vacuum (in CGS units)
virtual double pr_from_enth(double enth)=0
From the enthalpy, return the pressure.
double s_p
The value of the s-coordinate at the pole.
ubvector r_is_gp
Radial coordinate.
int new_search(int n, double *x, double val)
Search in array x of length n for value val.
double rho_equator_h
The value of at the equator.
void output_table(o2scl::table3d &t)
Create an output table.
double r_ratio
Ratio of polar to equatorial radius.
double eccentricity
The eccentricity.
double Omega_p
Angular velocity of a particle in a circular orbit at the equator.
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 .
void comp_f_P()
Compute two-point functions.
double Omega
Angular velocity (in )
int fix_cent_eden_bar_mass(double cent_eden, double bar_mass)
Construct a configuration with a fixed central energy density and a fixed baryonic mass.
void test3(o2scl::test_mgr &t)
Test fixed central energy density and fixed gravitational mass with EOS C.
double deriv_sm(ubmatrix &f, int s, int m)
Returns the derivative w.r.t. s and mu.
int fix_cent_eden_with_kepler_alt(double cent_eden, bool use_guess=false)
Experimental alternate form for fix_cent_eden_with_kepler()
ubvector e_d_gp
Energy density.
void test6(o2scl::test_mgr &t)
Test fixed central energy density and fixed angular momentum with EOS C.
double G
Gravitational constant (in CGS units)
void calc_masses_J(ubmatrix &rho_0)
Compute masses and angular momentum.
double deriv_s(ubmatrix &f, int s, int m)
Returns the derivative w.r.t. s
double Z_b
Backward equatorial redshift.
double DM
Spacing in direction, .
double om_over_Om
Ratio of potential to angular velocity .
ubmatrix S_gamma
source term in eqn for
int fix_cent_eden_grav_mass_alt(double cent_eden, double grav_mass, bool use_guess=false)
Experimental alternate form for fix_cent_eden_grav_mass()
double rho_center_h
The value of at the center.
int LMAX
The number of Legendre polynomials.
ubvector nu_gp
Metric function .
virtual double ed_from_nb(double nb)
From the baryon density, return the energy density.
int solve_ang_mom(size_t nv, const ubvector &x, ubvector &y, double ang_mom)
Solve for the gravitational mass.
double interp(double xp[], double yp[], int np, double xb)
Driver for the interpolation routine.
int fix_cent_eden_ang_mom(double cent_eden, double ang_mom)
Construct a configuration with a fixed central energy density and a fixed angular momentum.
ubmatrix pressure
Pressure.
double RMIN
Minimum radius for spherical stars (default )
ubmatrix gamma_guess
Guess for .
int fix_cent_eden_with_kepler(double cent_eden)
Construct a configuration with a fixed central energy density and the Keplerian rotation rate.
double omega_equator_h
The value of at the equator.
Rotating neutron star class based on RNS v1.1d from N. Stergioulas et al.
double log_e_tab[201]
rho points in tabulated EOS
A EOS base class for the TOV solver.
int MDIV
The number of grid points in the direction.
double p_center
Central pressure.
void comp()
Compute various quantities.
int n_nearest
Cache for interpolation.
Subclass of nstar_rot which specifies the function to invert a polytropic EOS.
void constants_o2scl()
Use the <a href='../../html/index.html'>O<span style='position: relative; top: 0.3em; font-size: 0....
int fix_cent_eden_ang_vel_alt(double cent_eden, double ang_vel, bool use_guess=false)
Experimental alternate form for fix_cent_eden_ang_vel()
double T
Total rotational kinetic energy.
o2scl::search_vec< double * > sv
Array search object.
double dm_dr_is(double r_is, double r, double m, double p)
Derivative of gravitational mass with respect to isotropic radius.
double Mass
Gravitational mass (in )
double deriv_m(ubmatrix &f, int s, int m)
Returns the derivative w.r.t. mu.
const double gravitational_constant
ubmatrix omega_guess
Guess for .
const double speed_of_light
double n0_at_e(double ee)
Baryon density at fixed energy density.
int n_nearest
Cache for interpolation.
ubmatrix energy
Energy density .
virtual double ed_from_pr(double pr)
From the pressure, return the energy density.
void test2(o2scl::test_mgr &t)
Test configuration rotating and Keplerian frequency with a fixed central energy density and EOS C.
double legendre(int n, double x)
Returns the Legendre polynomial of degree n, evaluated at x.
ubmatrix D2_omega
Integ. term over s in eqn for .
Create a tabulated EOS for nstar_rot using interpolation.
double MB
The mass of one baryon (in g)
double tol_abs
The tolerance for the functions with the prefix "fix" (default )
double W
Gravitational binding energy.
eos_nstar_rot * eosp
Pointer to the user-specified EOS.
o2scl::mroot * mrootp
Solver.
double h_center
Central enthalpy.
ubmatrix da_dm
Derivative of with respect to .
double e_at_p(double pp)
Compute
int fix_cent_eden_bar_mass_alt(double cent_eden, double bar_mass, bool use_guess=false)
Experimental alternate form for fix_cent_eden_bar_mass()
int fix_cent_eden_ang_mom_alt(double cent_eden, double ang_mom, bool use_guess=false)
Experimental alternate form for fix_cent_eden_ang_mom()
Documentation generated with Doxygen. Provided under the
GNU Free Documentation License (see License Information).