27#include <gsl/gsl_errno.h>
28#include <gsl/gsl_odeiv.h>
30#include <lal/LALStdlib.h>
31#include <lal/LALConstants.h>
32#include <lal/LALSimNeutronStar.h>
41 num = (8.0 / 5.0) * pow(1 - 2 *
c, 2.0) * pow(
c, 5)
42 * (2 *
c * (
y - 1) -
y + 2);
43 den = 2 *
c * (4 * (
y + 1) * pow(
c, 4) + (6 *
y - 4) * pow(
c, 3)
44 + (26 - 22 *
y) *
c *
c + 3 * (5 *
y - 8) *
c - 3 *
y + 6);
45 den -= 3 * pow(1 - 2 *
c, 2) * (2 *
c * (
y - 1) -
y + 2)
46 * log(1.0 / (1 - 2 *
c));
60#define TOV_ODE_VARS_DIM (sizeof(struct tov_ode_vars)/sizeof(double))
97 double A = 1.0 / (1.0 - 2.0 *
m /
r);
99 double C1 = 2.0 /
r + A * (2.0 *
m / (
r *
r) + 4.0 *
LAL_PI *
r * (
p -
e));
102 A * (-(2) * (2 + 1) / (
r *
r) + 4.0 *
LAL_PI * (
e +
p) * dedp +
103 4.0 *
LAL_PI * (5.0 *
e + 9.0 *
p)) - pow(2.0 * (
m +
105 double dr = -
r * (
r - 2.0 *
m) / (
m + 4.0 *
LAL_PI *
r *
r *
r *
p);
106 double dm = 4.0 *
LAL_PI *
r *
r *
e * dr;
108 double db = -(
C0 *
H +
C1 *
b) * dr;
134 double *love_number_k2,
double central_pressure_si,
138 const double epsabs = 0.0;
143 gsl_odeiv_step *step =
145 gsl_odeiv_control *ctrl = gsl_odeiv_control_y_new(epsabs, epsrel);
158 double dhdp_c = 1.0 / (ec +
pc);
159 double dedh_c = dedp_c / dhdp_c;
160 double dh = -1
e-12 * hc;
162 double h1 = 0.0 - dh;
163 double r0 = sqrt(-3.0 * dh / (2.0 *
LAL_PI * (ec + 3.0 *
pc)));
164 double m0 = 4.0 *
LAL_PI * r0 * r0 * r0 * ec / 3.0;
166 double b0 = 2.0 * r0;
176 r0 *= 1.0 + 0.25 * dh * (ec - 3.0 *
pc - 0.6 * dedh_c) / (ec + 3.0 *
pc);
178 m0 *= 1.0 + 0.6 * dh * dedh_c / ec;
188 gsl_odeiv_evolve_apply(evolv, ctrl, step, &sys, &h, h1, &dh,
y);
189 if (
s != GSL_SUCCESS)
191 "Error encountered in GSL's ODE integrator\n");
197 y[
i] += dy[
i] * (0.0 - h1);
200 c = vars->
m / vars->
r;
201 yy = vars->
r * vars->
b / vars->
H;
209 gsl_odeiv_evolve_free(evolv);
210 gsl_odeiv_control_free(ctrl);
211 gsl_odeiv_step_free(step);
228#define TOV_VIRIAL_ODE_VARS_DIM (sizeof(struct tov_virial_ode_vars)/sizeof(double))
266 double A = 1.0 / (1.0 - 2.0 *
m /
r);
268 double C1 = 2.0 /
r + A * (2.0 *
m / (
r *
r) + 4.0 *
LAL_PI *
r * (
p -
e));
271 A * (-(2) * (2 + 1) / (
r *
r) + 4.0 *
LAL_PI * (
e +
p) * dedp +
272 4.0 *
LAL_PI * (5.0 *
e + 9.0 *
p)) - pow(2.0 * (
m +
274 double dr = -
r * (
r - 2.0 *
m) / (
m + 4.0 *
LAL_PI *
r *
r *
r *
p);
275 double dm = 4.0 *
LAL_PI *
r *
r *
e * dr;
277 double db = -(
C0 *
H +
C1 *
b) * dr;
279 double alpha = 1.0 - 2.0 *
m /
r;
282 double dI1 = 8.0 *
LAL_PI *
r * pow(
alpha, (-1.0/2.0)) *
p * dr;
283 double dI2 =
r * pow(
alpha, (-1.5)) * pow(
beta, (2.0)) * dr;
284 double dJ1 = 4.0 *
LAL_PI *
r *
r * pow(
alpha, (-0.5)) * 3.0 *
p * dr;
285 double dJ2 = pow(
alpha, (-0.5)) * (pow(
alpha, (-1.0)) * pow((
beta *
r), (2.0)) - 0.5 * pow((sqrt(
alpha) - 1.0), (2.0))) * dr;
322 double *int1,
double *int2,
double *int3,
double *int4,
double *int5,
double *int6,
323 double *love_number_k2,
double central_pressure_si,
327 const double epsabs = 0.0;
332 gsl_odeiv_step *step =
334 gsl_odeiv_control *ctrl = gsl_odeiv_control_y_new(epsabs, epsrel);
347 double dhdp_c = 1.0 / (ec +
pc);
348 double dedh_c = dedp_c / dhdp_c;
349 double dh = -1
e-12 * hc;
351 double h1 = 0.0 - dh;
352 double r0 = sqrt(-3.0 * dh / (2.0 *
LAL_PI * (ec + 3.0 *
pc)));
353 double m0 = 4.0 *
LAL_PI * r0 * r0 * r0 * ec / 3.0;
355 double b0 = 2.0 * r0;
365 r0 *= 1.0 + 0.25 * dh * (ec - 3.0 *
pc - 0.6 * dedh_c) / (ec + 3.0 *
pc);
367 m0 *= 1.0 + 0.6 * dh * dedh_c / ec;
380 double I1_0 = - 8.0 *
LAL_PI * r0 * r0 *
pc * dh - dh * dh;
382 double J1_0 = - 12.0 *
LAL_PI * r0 * r0 * r0 *
pc * dh - 3.0 * r0 * dh * dh * dh;
399 gsl_odeiv_evolve_apply(evolv, ctrl, step, &sys, &h, h1, &dh,
y);
400 if (
s != GSL_SUCCESS)
402 "Error encountered in GSL's ODE integrator\n");
406 for (
int w = 0 ;
w < 1 ; ++
w){
409 y[
i] += dy[
i] * (0.0 - h1);
413 c = vars->
m / vars->
r;
414 yy = vars->
r * vars->
b / vars->
H;
416 *int3 = (1.0 - vars->
m / vars->
r) * pow((1.0 - 2.0 * vars->
m / vars->
r), (-0.5)) - 1.0;
417 *int6 = vars->
r * (*int3);
430 gsl_odeiv_evolve_free(evolv);
431 gsl_odeiv_control_free(ctrl);
432 gsl_odeiv_step_free(step);
439 const double epsrel = 1
e-6;
444 double *int1,
double *int2,
double *int3,
double *int4,
double *int5,
double *int6,
445 double *love_number_k2,
double central_pressure_si,
448 const double epsrel = 1
e-6;
450 int1, int2, int3, int4, int5, int6, love_number_k2, central_pressure_si, eos, epsrel);
static double beta(const double a, const double b, const sysq *system)
Internal function that computes the spin-orbit couplings.
static REAL8 UNUSED C1(REAL8 Mtotal)
static REAL8 UNUSED C0(REAL8 eta)
struct tagLALSimNeutronStarEOS LALSimNeutronStarEOS
Incomplete type for the neutron star Equation of State (EOS).
#define LAL_G_C4_SI
Factor to convert pressure in Pa to geometerized units of m^-2.
double XLALSimNeutronStarEOSEnergyDensityOfPressureGeometerized(double p, LALSimNeutronStarEOS *eos)
Returns the energy density in geometerized units (m^-2) at a given pressure in geometerized units (m^...
double XLALSimNeutronStarEOSEnergyDensityDerivOfPressureGeometerized(double p, LALSimNeutronStarEOS *eos)
Returns the gradient of the energy density with respect to the pressure (dimensionless) at a given va...
double XLALSimNeutronStarEOSEnergyDensityOfPseudoEnthalpyGeometerized(double h, LALSimNeutronStarEOS *eos)
Returns the energy density in geometerized units (m^-2) at a given value of the dimensionless pseudo-...
double XLALSimNeutronStarEOSPressureOfPseudoEnthalpyGeometerized(double h, LALSimNeutronStarEOS *eos)
Returns the pressure in geometerized units (m^-2) at a given value of the dimensionless pseudo-enthal...
double XLALSimNeutronStarEOSPseudoEnthalpyOfPressureGeometerized(double p, LALSimNeutronStarEOS *eos)
Returns the dimensionless pseudo-enthalpy at a given pressure in geometerized units (m^-2).
int XLALSimNeutronStarTOVODEIntegrate(double *radius, double *mass, double *love_number_k2, double central_pressure_si, LALSimNeutronStarEOS *eos)
static int tov_ode(double h, const double *y, double *dy, void *params)
int XLALSimNeutronStarVirialODEIntegrateWithTolerance(double *radius, double *mass, double *int1, double *int2, double *int3, double *int4, double *int5, double *int6, double *love_number_k2, double central_pressure_si, LALSimNeutronStarEOS *eos, double epsrel)
Integrates the Tolman-Oppenheimer-Volkov stellar structure equations and the Virial Equations.
static int tov_virial_ode(double h, const double *y, double *dy, void *params)
#define TOV_VIRIAL_ODE_VARS_DIM
static struct tov_virial_ode_vars * tov_virial_ode_vars_cast(const double *y)
int XLALSimNeutronStarTOVODEIntegrateWithTolerance(double *radius, double *mass, double *love_number_k2, double central_pressure_si, LALSimNeutronStarEOS *eos, double epsrel)
Integrates the Tolman-Oppenheimer-Volkov stellar structure equations.
int XLALSimNeutronStarVirialODEIntegrate(double *radius, double *mass, double *int1, double *int2, double *int3, double *int4, double *int5, double *int6, double *love_number_k2, double central_pressure_si, LALSimNeutronStarEOS *eos)
static double tidal_Love_number_k2(double c, double y)
static struct tov_ode_vars * tov_ode_vars_cast(const double *y)