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> 115 int new_search(
int n,
double *x,
double val);
120 double log_e_tab[201];
122 double log_p_tab[201];
124 double log_h_tab[201];
126 double log_n0_tab[201];
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>
168 double KAPPA=1.0e-15*C*C/G;
169 double KSCALE=KAPPA*G/(C*C*C*C);
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]);
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++) {
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;
302 log_h_tab[i+1]=log10(log(mu/mu_start));
304 log_n0_tab[i+1]=log10(nst_arr[i][2]);
309 log_h_tab[1]=log_h_tab[2]-8.0;
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]/
316 log_n0_tab[i+n_crust+1]=log10(nb[i]*1.0e39);
366 virtual void ed_nb_from_pr(
double pr,
double &ed,
double &nb);
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;
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,
569 int solve_kepler(
size_t nv,
const ubvector &x, ubvector &y);
572 int solve_grav_mass(
size_t nv,
const ubvector &x, ubvector &y,
576 int solve_bar_mass(
size_t nv,
const ubvector &x, ubvector &y,
580 int solve_ang_vel(
size_t nv,
const ubvector &x, ubvector &y,
584 int solve_ang_mom(
size_t nv,
const ubvector &x, ubvector &y,
615 return pow(rho0,_Gamma_P)/(_Gamma_P-1.0)+rho0-_ee;
833 int new_search(
int n, ubvector &x,
double val);
845 double interp(ubvector &xp, ubvector &yp,
int np,
double xb);
854 double interp_4_k(ubvector &xp, ubvector &yp,
int np,
double xb,
int k);
863 double int_z(ubvector &f,
int m);
872 double e_at_p(
double pp);
878 double p_at_e(
double ee);
884 double p_at_h(
double hh);
890 double h_at_p(
double pp);
896 double n0_at_e(
double ee);
903 double s_deriv(ubvector &f,
int s);
907 double m_deriv(ubvector &f,
int m);
911 double deriv_s(ubmatrix &f,
int s,
int m);
915 double deriv_m(ubmatrix &f,
int s,
int m);
919 double deriv_sm(ubmatrix &f,
int s,
int m);
929 double legendre(
int n,
double x);
983 void make_center(
double e_center);
1017 void spherical_star();
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);
1037 int iterate(
double r_ratio,
double tol_rel);
1056 void calc_masses_J(ubmatrix &rho_0);
1187 void constants_rns();
1190 void constants_o2scl();
1217 scaled_polytrope=
false;
1225 scaled_polytrope=
true;
1238 int fix_cent_eden_axis_rat(
double cent_eden,
double axis_rat,
1239 bool use_guess=
false);
1248 int fix_cent_eden_grav_mass(
double cent_eden,
double grav_mass);
1257 int fix_cent_eden_bar_mass(
double cent_eden,
double bar_mass);
1265 int fix_cent_eden_with_kepler(
double cent_eden);
1270 int fix_cent_eden_with_kepler_alt(
double cent_eden,
1271 bool use_guess=
false);
1276 int fix_cent_eden_grav_mass_alt(
double cent_eden,
double grav_mass,
1277 bool use_guess=
false);
1282 int fix_cent_eden_bar_mass_alt(
double cent_eden,
double bar_mass,
1283 bool use_guess=
false);
1288 int fix_cent_eden_ang_vel_alt(
double cent_eden,
double ang_vel,
1289 bool use_guess=
false);
1294 int fix_cent_eden_ang_mom_alt(
double cent_eden,
double ang_mom,
1295 bool use_guess=
false);
1303 int fix_cent_eden_non_rot(
double cent_eden);
1318 int fix_cent_eden_ang_vel(
double cent_eden,
double ang_vel);
1327 int fix_cent_eden_ang_mom(
double cent_eden,
double ang_mom);
ubvector theta
defined by
double rho_center_h
The value of at the center.
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.
double h_minus
Height from surface of last stable counter-rotating circular orbit in equatorial plane.
ubvector m_gp
Enclosed gravitational mass.
ubmatrix P1_2n_1
Associated Legendre polynomial .
double MB
The mass of one baryon (in g)
double p_center
Central pressure.
double Mass_p
Proper mass (in )
virtual double nb_from_ed(double ed)=0
From the energy density, return the baryon density.
ubmatrix pressure
Pressure.
bool scaled_polytrope
If true, then use a polytrope and rescale.
double Z_f
Forward equatorial redshift.
double gamma_center_h
The value of at the center.
o2scl::mroot * mrootp
Solver.
double gamma_equator_h
The value of at the equator.
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) ...
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 .
void polytrope_eos(double index)
Use a polytropic EOS with a specified index.
double Gamma_P
Polytropic exponent.
double mass_quadrupole
Mass quadrupole moment.
lib_settings_class o2scl_settings
ubvector lambda_gp
Metric function .
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 )
ubmatrix S_gamma
source term in eqn for
ubvector e_d_gp
Energy density.
int RDIV
The number of grid points in integration of TOV equations for spherical stars.
ubvector nu_gp
Metric function .
double cf
The convergence factor (default 1.0)
double eccentricity
The eccentricity.
double Omega_p
Angular velocity of a particle in a circular orbit at the equator.
ubmatrix D2_omega
Integ. term over s in eqn for .
double Omega
Angular velocity (in )
void output()
Output EOS table to screen.
o2scl::root_bkt_cern< polytrope_solve > rbc
The polytrope solver.
double C
Speed of light in vacuum (in CGS units)
double KAPPA
Square of length scale in CGS units, .
double rho_equator_h
The value of at the equator.
double J
Angular momentum (in )
double DS
Spacing in direction, .
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 .
int n_nearest
Cache for interpolation.
double h_center
Central enthalpy.
int set_solver(o2scl::mroot<> &m)
Set new solver.
double T
Total rotational kinetic energy.
ubmatrix velocity_sq
Proper velocity squared.
double G
Gravitational constant (in CGS units)
ubmatrix S_rho
source term in eqn for
double Z_b
Backward equatorial redshift.
ubmatrix omega_guess
Guess for .
double I
Moment of inertia.
ubvector r_gp
Isotropic radius.
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.
polytrope_solve(double Gamma_P, double ee)
Create a function object with specified polytropic index and ?
double r_e_guess
Guess for the equatorial radius.
virtual double enth_from_pr(double pr)=0
From the pressure, return the enthalpy.
double W
Gravitational binding energy.
double r_p
Radius at pole.
double omega_equator_h
The value of at the equator.
Tabulated EOS for nstar_rot from Bethe74.
int n_tab
number of tabulated EOS points
double h_plus
Height from surface of last stable co-rotating circular orbit in equatorial plane.
Tabulated EOS for nstar_rot from Pandharipande75.
const double speed_of_light
ubmatrix da_dm
Derivative of with respect to .
double rho_pole_h
The value of at the pole.
Create a tabulated EOS for nstar_rot using interpolation.
#define O2SCL_ERR2(d, d2, n)
double alt_tol_rel
Accuracy for equatorial radius using alternate solvers (default )
double MSUN
Mass of sun (in g)
A EOS base class for the TOV solver.
int SDIV
The number of grid points in the direction.
double _ee
The energy density.
void set_eos(eos_nstar_rot &eos)
Set the EOS.
ubmatrix enthalpy
Enthalpy.
double enthalpy_min
Minimum specific enthalpy.
double Omega_h
Angular velocity, .
int verbose
Verbosity parameter.
double e_center
Central energy density (in units of )
double s_p
The value of the s-coordinate at the pole.
double _Gamma_P
The polytropic index.
double p_surface
Pressure at the surface.
int n_nearest
Cache for interpolation.
double Z_p
Polar redshift.
ubmatrix energy
Energy density .
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 .
Subclass of nstar_rot which specifies the function to invert a polytropic EOS.
o2scl::search_vec< double * > sv
Array search object.
double C
Speed of light in vacuum (in CGS units)
o2scl::mroot_hybrids def_mroot
Default solver.
double eq_radius_tol_rel
Relative accuracy for the equatorial radius, (default )
ubmatrix gamma_guess
Guess for .
double SMAX
Maximum value of s-coordinate (default 0.9999)
ubmatrix D2_rho
Integrated term over s in eqn for .
bool eos_set
If true, then an EOS has been set.
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 ...
o2scl::search_vec< ubvector_range > sv_ub
Array search object.
int MDIV
The number of grid points in the direction.
double G
Gravitational constant (in CGS units)
double s_e
The value of the s-coordinate at the equator.
double gamma_pole_h
The value of at the pole.
double Mass
Gravitational mass (in )
ubmatrix alpha_guess
Guess for .
double operator()(double rho0)
The function.
int LMAX
The number of Legendre polynomials.
eos_nstar_rot * eosp
Pointer to the user-specified EOS.
double r_ratio
Ratio of polar to equatorial radius.
ubmatrix S_omega
source term in eqn for
double velocity_equator
The velocity at the equator.
double Mass_0
Baryonic mass (in )
ubmatrix D2_gamma
Integrated term over s in eqn for .
double n_P
Polytropic index.
double e_surface
Energy density at the surface.
ubmatrix D1_gamma
Integrated term over m in eqn for .
double r_e
Coordinate equatorial radius.
ubmatrix P_2n
Legendre polynomial .
virtual double pr_from_nb(double nb)=0
From the baryon density, return the pressure.
double DM
Spacing in direction, .
double Omega_K
Kepler rotation frequency (in 1/s)
double KAPPA
Square of length scale in CGS units, .
const double gravitational_constant
double RMIN
Minimum radius for spherical stars (default )
ubvector r_is_gp
Radial coordinate.
double om_over_Om
Ratio of potential to angular velocity .