Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALSimulation 6.2.0.1-5e288d3
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
LALSimNeutronStarTOV.c
Go to the documentation of this file.
1/*
2 * Copyright (C) 2013 J. Creighton, B. Lackey
3 *
4 * This program is free software; you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License as published by
6 * the Free Software Foundation; either version 2 of the License, or
7 * (at your option) any later version.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with with program; see the file COPYING. If not, write to the
16 * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
17 * MA 02110-1301 USA
18 */
19/**
20 * @author Jolien Creighton, Benjamin Lackey
21 * @addtogroup LALSimNeutronStarTOV_c
22 * @brief Provides routines for solving the Tolman-Oppenheimer-Volkov equation.
23 * @{
24 */
25
26#include <math.h>
27#include <gsl/gsl_errno.h>
28#include <gsl/gsl_odeiv.h>
29
30#include <lal/LALStdlib.h>
31#include <lal/LALConstants.h>
32#include <lal/LALSimNeutronStar.h>
33
34/* Implements Eq. (50) of Damour & Nagar, Phys. Rev. D 80 084035 (2009).
35 * See also Eq. (14) of Hinderer et al. Phys. Rev. D 81 123016 (2010). */
36static double tidal_Love_number_k2(double c, double y)
37{
38 double num;
39 double den;
40
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));
47
48 return num / den;
49}
50
51/* For convenience, use a structure that provides a dictionary between
52 * a vector of the ode variables and what they represent. */
54 double r; /* radial coordinate, m */
55 double m; /* mass within r in geometric units, m */
56 double H; /* stellar perturbation in arbitrary units */
57 double b; /* derivative of metric pertubation in arbitrary units */
58};
59
60#define TOV_ODE_VARS_DIM (sizeof(struct tov_ode_vars)/sizeof(double))
61
62/* Casts an array of doubles to a structure with named parameter. */
63static struct tov_ode_vars *tov_ode_vars_cast(const double *y)
64{
65 union {
66 const double *y;
67 struct tov_ode_vars *v;
68 } u = {
69 y};
70 return u.v;
71}
72
73/* ODE integrand for TOV equations with pseudo-enthalpy independent variable.
74 * Implements Eqs. (5) and (6) of Lindblom, Astrophys. J. 398, 569 (1992).
75 * Also uses Eqs. (7) and (8) [ibid] for inner boundary data, and
76 * Eqs. (18), (27), (28) of Damour & Nagar, Phys. Rev. D 80, 084035 (2009)
77 * [See also: Eqs. (11) & (12) Hinderer et al. Phys. Rev. D 81 123016 (2010)]
78 * for the metric perturbation used to obtain the Love number. */
79static int tov_ode(double h, const double *y, double *dy, void *params)
80{
81 struct tov_ode_vars *vars = tov_ode_vars_cast(y);
82 struct tov_ode_vars *derivs = tov_ode_vars_cast(dy);
84
85 double r = vars->r;
86 double m = vars->m;
87 double H = vars->H;
88 double b = vars->b;
89 double p =
91 double e =
93 eos);
94 double dedp =
96 /* Eq. (18) of Damour & Nagar PRD 80 084035 (2009). */
97 double A = 1.0 / (1.0 - 2.0 * m / r);
98 /* Eq. (28) of Damour & Nagar PRD 80 084035 (2009). */
99 double C1 = 2.0 / r + A * (2.0 * m / (r * r) + 4.0 * LAL_PI * r * (p - e));
100 /* Eq. (29) of Damour & Nagar PRD 80 084035 (2009). */
101 double C0 =
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 +
104 4.0 * LAL_PI * r * r * r * p) / (r * (r - 2.0 * m)), 2.0);
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;
107 double dH = b * dr;
108 double db = -(C0 * H + C1 * b) * dr;
109
110 derivs->r = dr;
111 derivs->m = dm;
112 derivs->H = dH;
113 derivs->b = db;
114 return 0;
115}
116
117/**
118 * @brief Integrates the Tolman-Oppenheimer-Volkov stellar structure equations.
119 * @details
120 * Solves the Tolman-Oppenheimer-Volkov stellar structure equations using the
121 * pseudo-enthalpy formalism introduced in:
122 * Lindblom (1992) "Determining the Nuclear Equation of State from Neutron-Star
123 * Masses and Radii", Astrophys. J. 398 569.
124 * @param[out] radius The radius of the star in m.
125 * @param[out] mass The mass of the star in kg.
126 * @param[out] love_number_k2 The k_2 tidal love number of the star.
127 * @param[in] central_pressure_si The central pressure of the star in Pa.
128 * @param eos Pointer to the Equation of State structure.
129 * @param[in] epsrel The relative error for the TOV solver routine
130 * @retval 0 Success.
131 * @retval <0 Failure.
132 */
133int XLALSimNeutronStarTOVODEIntegrateWithTolerance(double *radius, double *mass,
134 double *love_number_k2, double central_pressure_si,
135 LALSimNeutronStarEOS * eos, double epsrel)
136{
137 /* ode integration variables */
138 const double epsabs = 0.0;
139 double y[TOV_ODE_VARS_DIM];
140 double dy[TOV_ODE_VARS_DIM];
141 struct tov_ode_vars *vars = tov_ode_vars_cast(y);
142 gsl_odeiv_system sys = { tov_ode, NULL, TOV_ODE_VARS_DIM, eos };
143 gsl_odeiv_step *step =
144 gsl_odeiv_step_alloc(gsl_odeiv_step_rk8pd, TOV_ODE_VARS_DIM);
145 gsl_odeiv_control *ctrl = gsl_odeiv_control_y_new(epsabs, epsrel);
146 gsl_odeiv_evolve *evolv = gsl_odeiv_evolve_alloc(TOV_ODE_VARS_DIM);
147
148 /* central values */
149 /* note: will be updated with Lindblom's series expansion */
150 /* geometrisized units for variables in length (m) */
151 double pc = central_pressure_si * LAL_G_C4_SI;
152 double ec =
154 double hc =
156 double dedp_c =
158 double dhdp_c = 1.0 / (ec + pc);
159 double dedh_c = dedp_c / dhdp_c;
160 double dh = -1e-12 * hc;
161 double h0 = hc + dh;
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;
165 double H0 = r0 * r0;
166 double b0 = 2.0 * r0;
167
168 double yy;
169 double c;
170 double h;
171 size_t i;
172
173 /* series expansion for the initial core */
174
175 /* second factor of Eq. (7) of Lindblom (1992) */
176 r0 *= 1.0 + 0.25 * dh * (ec - 3.0 * pc - 0.6 * dedh_c) / (ec + 3.0 * pc);
177 /* second factor of Eq. (8) of Lindblom (1992) */
178 m0 *= 1.0 + 0.6 * dh * dedh_c / ec;
179
180 /* perform integration */
181 vars->r = r0;
182 vars->m = m0;
183 vars->H = H0;
184 vars->b = b0;
185 h = h0;
186 while (h > h1) {
187 int s =
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");
192 }
193
194 /* take one final Euler step to get to surface */
195 tov_ode(h, y, dy, eos);
196 for (i = 0; i < TOV_ODE_VARS_DIM; ++i)
197 y[i] += dy[i] * (0.0 - h1);
198
199 /* compute tidal Love number k2 */
200 c = vars->m / vars->r; /* compactness */
201 yy = vars->r * vars->b / vars->H;
202
203 /* convert from geometric units to SI units */
204 *radius = vars->r;
205 *mass = vars->m * LAL_MSUN_SI / LAL_MRSUN_SI;
206 *love_number_k2 = tidal_Love_number_k2(c, yy);
207
208 /* free ode memory */
209 gsl_odeiv_evolve_free(evolv);
210 gsl_odeiv_control_free(ctrl);
211 gsl_odeiv_step_free(step);
212 return 0;
213}
214
215/* For convenience, use a structure that provides a dictionary between
216 * a vector of the ode variables and what they represent. */
218 double r; /* radial coordinate, m */
219 double m; /* mass within r in geometric units, m */
220 double H; /* stellar perturbation in arbitrary units */
221 double b; /* derivative of metric pertubation in arbitrary units */
222 double I1; /* dependent variable for Virial ODEs */
223 double I2; /* dependent variable for Virial ODEs */
224 double J1; /* dependent variable for Virial ODEs */
225 double J2; /* dependent variable for Virial ODEs */
226};
227
228#define TOV_VIRIAL_ODE_VARS_DIM (sizeof(struct tov_virial_ode_vars)/sizeof(double))
229
230/* Casts an array of doubles to a structure with named parameter. */
231static struct tov_virial_ode_vars *tov_virial_ode_vars_cast(const double *y)
232{
233 union {
234 const double *y;
235 struct tov_virial_ode_vars *v;
236 } u = {
237 y};
238 return u.v;
239}
240
241/* ODE integrand for TOV equations and Virial equations with pseudo-enthalpy independent variable.
242 * Implements Eqs. (5) and (6) of Lindblom, Astrophys. J. 398, 569 (1992).
243 * Also uses Eqs. (7) and (8) [ibid] for inner boundary data, and
244 * Eqs. (18), (27), (28) of Damour & Nagar, Phys. Rev. D 80, 084035 (2009)
245 * [See also: Eqs. (11) & (12) Hinderer et al. Phys. Rev. D 81 123016 (2010)]
246 * for the metric perturbation used to obtain the Love number.
247 * For the Virial portion, implements equations provided by A. Nikolaidis, N. Stergioulas, H. Markakis. */
248static int tov_virial_ode(double h, const double *y, double *dy, void *params)
249{
253
254 double r = vars->r;
255 double m = vars->m;
256 double H = vars->H;
257 double b = vars->b;
258 double p =
260 double e =
262 eos);
263 double dedp =
265 /* Eq. (18) of Damour & Nagar PRD 80 084035 (2009). */
266 double A = 1.0 / (1.0 - 2.0 * m / r);
267 /* Eq. (28) of Damour & Nagar PRD 80 084035 (2009). */
268 double C1 = 2.0 / r + A * (2.0 * m / (r * r) + 4.0 * LAL_PI * r * (p - e));
269 /* Eq. (29) of Damour & Nagar PRD 80 084035 (2009). */
270 double C0 =
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 +
273 4.0 * LAL_PI * r * r * r * p) / (r * (r - 2.0 * m)), 2.0);
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;
276 double dH = b * dr;
277 double db = -(C0 * H + C1 * b) * dr;
278
279 double alpha = 1.0 - 2.0 * m / r;
280 double beta = (m + 4.0 * LAL_PI * r * r * r * p) / (r * r);
281
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;
286
287
288 derivs->r = dr;
289 derivs->m = dm;
290 derivs->H = dH;
291 derivs->b = db;
292 derivs->I1 = dI1;
293 derivs->I2 = dI2;
294 derivs->J1 = dJ1;
295 derivs->J2 = dJ2;
296 return 0;
297}
298
299/**
300 * @brief Integrates the Tolman-Oppenheimer-Volkov stellar structure equations and the Virial Equations.
301 * @details
302 * Solves the Tolman-Oppenheimer-Volkov stellar structure equations using the
303 * pseudo-enthalpy formalism introduced in:
304 * Lindblom (1992) "Determining the Nuclear Equation of State from Neutron-Star
305 * Masses and Radii", Astrophys. J. 398 569.
306 * @param[out] radius The radius of the star in m.
307 * @param[out] mass The mass of the star in kg.
308 * @param[out] int1 Virial parameter.
309 * @param[out] int2 Virial parameter.
310 * @param[out] int3 Virial parameter.
311 * @param[out] int4 Virial parameter.
312 * @param[out] int5 Virial parameter.
313 * @param[out] int6 Virial parameter.
314 * @param[out] love_number_k2 The k_2 tidal love number of the star.
315 * @param[in] central_pressure_si The central pressure of the star in Pa.
316 * @param eos Pointer to the Equation of State structure.
317 * @param[in] epsrel The relative error in the TOV solver routine.
318 * @retval 0 Success.
319 * @retval <0 Failure.
320 */
322 double *int1, double *int2, double *int3, double *int4, double *int5, double *int6,
323 double *love_number_k2, double central_pressure_si,
324 LALSimNeutronStarEOS * eos, double epsrel)
325{
326 /* ode integration variables */
327 const double epsabs = 0.0;
329 double dy[TOV_VIRIAL_ODE_VARS_DIM];
331 gsl_odeiv_system sys = { tov_virial_ode, NULL, TOV_VIRIAL_ODE_VARS_DIM, eos };
332 gsl_odeiv_step *step =
333 gsl_odeiv_step_alloc(gsl_odeiv_step_rk8pd, TOV_VIRIAL_ODE_VARS_DIM);
334 gsl_odeiv_control *ctrl = gsl_odeiv_control_y_new(epsabs, epsrel);
335 gsl_odeiv_evolve *evolv = gsl_odeiv_evolve_alloc(TOV_VIRIAL_ODE_VARS_DIM);
336
337 /* central values */
338 /* note: will be updated with Lindblom's series expansion */
339 /* geometrisized units for variables in length (m) */
340 double pc = central_pressure_si * LAL_G_C4_SI;
341 double ec =
343 double hc =
345 double dedp_c =
347 double dhdp_c = 1.0 / (ec + pc);
348 double dedh_c = dedp_c / dhdp_c;
349 double dh = -1e-12 * hc;
350 double h0 = hc + dh;
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;
354 double H0 = r0 * r0;
355 double b0 = 2.0 * r0;
356
357 double yy;
358 double c;
359 double h;
360 size_t i;
361
362 /* series expansion for the initial core */
363
364 /* second factor of Eq. (7) of Lindblom (1992) */
365 r0 *= 1.0 + 0.25 * dh * (ec - 3.0 * pc - 0.6 * dedh_c) / (ec + 3.0 * pc);
366 /* second factor of Eq. (8) of Lindblom (1992) */
367 m0 *= 1.0 + 0.6 * dh * dedh_c / ec;
368
369 // double Gamma_c = (ec + pc)/(pc * dedp_c);
370 // double r1 = sqrt(3.0/(2.0 * LAL_PI * (ec + 3.0 * pc)));
371 // double r3 = - (r1 / (4.0 * (ec + 3.0 * pc))) * (ec - 3.0 * pc - 3.0 * (ec + pc) * (ec + pc) / (5.0 * pc * Gamma_c));
372 // double m3 = 4.0 * LAL_PI * ec * r1 * r1 * r1 / 3.0;
373 // double m5 = 4.0 * LAL_PI * r1 * r1 * r1 * (r3 * ec / r1 - (ec + pc) * (ec + pc) / (5.0 * pc * Gamma_c));
374
375 // r0 = r1 * sqrt(fabs(dh)) + r3 * pow(sqrt(fabs(dh)), 3.0);
376 // m0 = m3 * pow(sqrt(fabs(dh)), 3.0) + m5 * pow(sqrt(fabs(dh)), 5.0);
377
378 /* Virial ODEs starting points */
379
380 double I1_0 = - 8.0 * LAL_PI * r0 * r0 * pc * dh - dh * dh;
381 double I2_0 = - 16.0 * LAL_PI * LAL_PI * r0 * r0 * r0 * r0 * pc * pc * dh - 6.0 * LAL_PI * r0 * r0 * pc * dh * dh;
382 double J1_0 = - 12.0 * LAL_PI * r0 * r0 * r0 * pc * dh - 3.0 * r0 * dh * dh * dh;
383 double J2_0 = - 16.0 * LAL_PI * LAL_PI * r0 * r0 * r0 * r0 * r0 * pc * pc * dh + 8.0 * LAL_PI * LAL_PI * r0 * r0 * r0 * r0 * r0 * pc * pc * dh * dh;
384
385
386 /* perform integration */
387 vars->r = r0;
388 vars->m = m0;
389 vars->H = H0;
390 vars->b = b0;
391 vars->I1 = I1_0;
392 vars->I2 = I2_0;
393 vars->J1 = J1_0;
394 vars->J2 = J2_0;
395
396 h = h0;
397 while (h > h1) {
398 int s =
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");
403 }
404
405 /*take one final Euler step to get to surface*/
406 for (int w = 0 ; w < 1 ; ++w){
407 tov_virial_ode(h, y, dy, eos);
408 for (i = 0; i < TOV_ODE_VARS_DIM; ++i)
409 y[i] += dy[i] * (0.0 - h1);
410 }
411
412 /* compute tidal Love number k2 */
413 c = vars->m / vars->r; /* compactness */
414 yy = vars->r * vars->b / vars->H;
415
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);
418
419 *int1 = vars->I1;
420 *int2 = vars->I2;
421 *int4 = vars->J1;
422 *int5 = vars->J2;
423
424 /* convert from geometric units to SI units */
425 *radius = vars->r;
426 *mass = vars->m * LAL_MSUN_SI / LAL_MRSUN_SI;
427 *love_number_k2 = tidal_Love_number_k2(c, yy);
428
429 /* free ode memory */
430 gsl_odeiv_evolve_free(evolv);
431 gsl_odeiv_control_free(ctrl);
432 gsl_odeiv_step_free(step);
433 return 0;
434}
435
436int XLALSimNeutronStarTOVODEIntegrate(double *radius, double *mass,
437 double *love_number_k2, double central_pressure_si, LALSimNeutronStarEOS * eos)
438{
439 const double epsrel = 1e-6;
440 return XLALSimNeutronStarTOVODEIntegrateWithTolerance(radius, mass, love_number_k2, central_pressure_si, eos, epsrel);
441}
442
443int XLALSimNeutronStarVirialODEIntegrate(double *radius, double *mass,
444 double *int1, double *int2, double *int3, double *int4, double *int5, double *int6,
445 double *love_number_k2, double central_pressure_si,
447{
448 const double epsrel = 1e-6;
450 int1, int2, int3, int4, int5, int6, love_number_k2, central_pressure_si, eos, epsrel);
451}
452
453/** @} */
const double b0
static double beta(const double a, const double b, const sysq *system)
Internal function that computes the spin-orbit couplings.
#define pc
#define c
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.
int s
Definition: bh_qnmode.c:137
double i
Definition: bh_ringdown.c:118
double e
Definition: bh_ringdown.c:117
const double H
const double u
const double w
#define LAL_MSUN_SI
#define LAL_PI
#define LAL_MRSUN_SI
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)
#define TOV_ODE_VARS_DIM
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)
static const INT4 r
static const INT4 m
#define XLAL_ERROR(...)
XLAL_EERR
list y
p
double alpha
Definition: sgwb.c:106
Definition: burst.c:245