Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALSimulation 6.2.0.1-b246709
•All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
LALSimNeutronStarEOSTabular.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 * @addtogroup LALSimNeutronStarEOS_c
21 * @{
22 */
23/**
24 * @name Creation routines for tabulated equations of state
25 * @{
26 */
27
28#include <lal/LALSimReadData.h>
29#include <gsl/gsl_interp.h>
30#include <lal/LALSimNeutronStar.h>
31
32/** @cond */
33
34/* Contents of the tabular equation of state data structure. */
35struct tagLALSimNeutronStarEOSDataTabular {
36 double *nbdat;
37 double *log_edat;
38 double *log_pdat;
39 double *mubdat;
40 double *muedat;
41 double *log_hdat;
42 double *yedat;
43 double *log_cs2dat;
44 double *log_rhodat;
45 size_t ncol;
46 size_t ndat;
47
48 gsl_interp *log_e_of_log_p_interp;
49 gsl_interp *log_h_of_log_p_interp;
50 gsl_interp *log_e_of_log_h_interp;
51 gsl_interp *log_p_of_log_h_interp;
52 gsl_interp *log_rho_of_log_h_interp;
53 gsl_interp *log_p_of_log_e_interp;
54 gsl_interp *log_p_of_log_rho_interp;
55 gsl_interp *log_cs2_of_log_h_interp;
56 gsl_interp_accel *log_e_of_log_p_acc;
57 gsl_interp_accel *log_h_of_log_p_acc;
58 gsl_interp_accel *log_e_of_log_h_acc;
59 gsl_interp_accel *log_p_of_log_h_acc;
60 gsl_interp_accel *log_rho_of_log_h_acc;
61 gsl_interp_accel *log_p_of_log_e_acc;
62 gsl_interp_accel *log_p_of_log_rho_acc;
63 gsl_interp_accel *log_cs2_of_log_h_acc;
64};
65
66static double eos_p_of_e_tabular(double e, LALSimNeutronStarEOS * eos)
67{
68 double log_e;
69 double log_p;
70 if (e == 0.0)
71 return 0.0;
72 log_e = log(e);
73 if (log_e < eos->data.tabular->log_edat[0])
74 /* use non-relativistic degenerate gas, p = K * e**(5./3.) */
75 return exp(eos->data.tabular->log_pdat[0] + (5.0 / 3.0) * (log_e - eos->data.tabular->log_edat[0]));
76 log_p = gsl_interp_eval(eos->data.tabular->log_p_of_log_e_interp,
77 eos->data.tabular->log_edat, eos->data.tabular->log_pdat, log_e,
78 eos->data.tabular->log_p_of_log_e_acc);
79 return exp(log_p);
80}
81
82static double eos_p_of_rho_tabular(double rho, LALSimNeutronStarEOS * eos)
83{
84 double log_rho;
85 double log_p;
86 if (rho == 0.0)
87 return 0.0;
88 log_rho = log(rho);
89 if (log_rho < eos->data.tabular->log_rhodat[0])
90 /* use non-relativistic degenerate gas, p = K * rho**(5./3.) */
91 return exp(eos->data.tabular->log_pdat[0] + (5.0 / 3.0) * (log_rho - eos->data.tabular->log_rhodat[0]));
92 log_p = gsl_interp_eval(eos->data.tabular->log_p_of_log_rho_interp,
93 eos->data.tabular->log_rhodat, eos->data.tabular->log_pdat, log_rho,
94 eos->data.tabular->log_p_of_log_rho_acc);
95 return exp(log_p);
96}
97
98static double eos_e_of_p_tabular(double p, LALSimNeutronStarEOS * eos)
99{
100 double log_p;
101 double log_e;
102 if (p == 0.0)
103 return 0.0;
104 log_p = log(p);
105 if (log_p < eos->data.tabular->log_pdat[0])
106 /* use non-relativistic degenerate gas, p = K * e**(5./3.) */
107 return exp(eos->data.tabular->log_edat[0] + (3.0 / 5.0) * (log_p - eos->data.tabular->log_pdat[0]));
108 log_e = gsl_interp_eval(eos->data.tabular->log_e_of_log_p_interp,
109 eos->data.tabular->log_pdat, eos->data.tabular->log_edat, log_p,
110 eos->data.tabular->log_e_of_log_p_acc);
111 return exp(log_e);
112}
113
114static double eos_e_of_h_tabular(double h, LALSimNeutronStarEOS * eos)
115{
116 double log_h;
117 double log_e;
118 if (h == 0.0)
119 return 0.0;
120 log_h = log(h);
121 if (log_h < eos->data.tabular->log_hdat[0])
122 /* use non-relativistic degenerate gas, e = K * h**(3./2.) */
123 return exp(eos->data.tabular->log_edat[0] + 1.5 * (log_h - eos->data.tabular->log_hdat[0]));
124 log_e = gsl_interp_eval(eos->data.tabular->log_e_of_log_h_interp,
125 eos->data.tabular->log_hdat, eos->data.tabular->log_edat, log_h,
126 eos->data.tabular->log_e_of_log_h_acc);
127 return exp(log_e);
128}
129
130static double eos_p_of_h_tabular(double h, LALSimNeutronStarEOS * eos)
131{
132 double log_h;
133 double log_p;
134 if (h == 0.0)
135 return 0.0;
136 log_h = log(h);
137 if (log_h < eos->data.tabular->log_hdat[0])
138 /* use non-relativistic degenerate gas, p = K * h**(5./2.) */
139 return exp(eos->data.tabular->log_pdat[0] + 2.5 * (log_h - eos->data.tabular->log_hdat[0]));
140 log_p = gsl_interp_eval(eos->data.tabular->log_p_of_log_h_interp,
141 eos->data.tabular->log_hdat, eos->data.tabular->log_pdat, log_h,
142 eos->data.tabular->log_p_of_log_h_acc);
143 return exp(log_p);
144}
145
146static double eos_rho_of_h_tabular(double h, LALSimNeutronStarEOS * eos)
147{
148 double log_h;
149 double log_rho;
150 if (h == 0.0)
151 return 0.0;
152 log_h = log(h);
153 if (log_h < eos->data.tabular->log_hdat[0])
154 /* use non-relativistic degenerate gas, rho = K * h**(3./2.) */
155 return exp(eos->data.tabular->log_rhodat[0] + 1.5 * (log_h - eos->data.tabular->log_hdat[0]));
156 log_rho = gsl_interp_eval(eos->data.tabular->log_rho_of_log_h_interp,
157 eos->data.tabular->log_hdat, eos->data.tabular->log_rhodat, log_h,
158 eos->data.tabular->log_rho_of_log_h_acc);
159 return exp(log_rho);
160}
161
162static double eos_h_of_p_tabular(double p, LALSimNeutronStarEOS * eos)
163{
164 double log_p;
165 double log_h;
166 if (p == 0)
167 return 0.0;
168 log_p = log(p);
169 if (log_p < eos->data.tabular->log_pdat[0])
170 /* use non-relativistic degenerate gas, h = K * p**(2./5.) */
171 return exp(eos->data.tabular->log_hdat[0] + 0.4 * (log_p - eos->data.tabular->log_pdat[0]));
172 log_h = gsl_interp_eval(eos->data.tabular->log_h_of_log_p_interp,
173 eos->data.tabular->log_pdat, eos->data.tabular->log_hdat, log_p,
174 eos->data.tabular->log_h_of_log_p_acc);
175 return exp(log_h);
176}
177
178static double eos_dedp_of_p_tabular(double p, LALSimNeutronStarEOS * eos)
179{
180 double log_p;
181 double log_e;
182 double d_log_e_d_log_p;
183 if (p == 0 || (log_p = log(p)) < eos->data.tabular->log_pdat[0])
184 /* use non-relativistic degenerate gas, p = K * e**(5./3.) */
185 return (3.0 / 5.0) * exp(eos->data.tabular->log_edat[0] - eos->data.tabular->log_pdat[0]);
186 log_e = gsl_interp_eval(eos->data.tabular->log_e_of_log_p_interp,
187 eos->data.tabular->log_pdat, eos->data.tabular->log_edat, log_p,
188 eos->data.tabular->log_e_of_log_p_acc);
189 d_log_e_d_log_p = gsl_interp_eval_deriv(eos->data.tabular->log_e_of_log_p_interp,
190 eos->data.tabular->log_pdat, eos->data.tabular->log_edat, log_p,
191 eos->data.tabular->log_e_of_log_p_acc);
192 return d_log_e_d_log_p * exp(log_e - log_p);
193}
194
195static double eos_v_of_h_tabular(double h, LALSimNeutronStarEOS * eos)
196{
197 double p, dedp, log_cs2, log_h;
198 if (eos->data.tabular->ncol == 2)
199 {
200 p = eos_p_of_h_tabular(h, eos);
201 dedp = eos_dedp_of_p_tabular(p, eos);
202 return pow(dedp, -0.5);
203 }
204 log_h = log(h);
205 log_cs2 = gsl_interp_eval(eos->data.tabular->log_cs2_of_log_h_interp,
206 eos->data.tabular->log_hdat, eos->data.tabular->log_cs2dat, log_h,
207 eos->data.tabular->log_cs2_of_log_h_acc);
208
209 return pow(exp(log_cs2), -0.5);
210}
211
212//static double eos_v_of_h_tabular(double h, LALSimNeutronStarEOS *eos)
213//{
214// double dpdh, dedh;
215// printf("hi4\n");
216// dpdh = gsl_interp_eval_deriv(eos->data.tabular->log_p_of_log_h_interp, eos->data.tabular->hdat, eos->data.tabular->pdat, h, eos->data.tabular->p_of_h_acc);
217// printf("hi5\n");
218// dedh = gsl_interp_eval_deriv(eos->data.tabular->e_of_h_interp, eos->data.tabular->hdat, eos->data.tabular->edat, h, eos->data.tabular->e_of_h_acc);
219// return sqrt(dpdh/dedh);
220//}
221
222static void eos_free_tabular_data(LALSimNeutronStarEOSDataTabular * data)
223{
224 if (data) {
225 gsl_interp_free(data->log_e_of_log_p_interp);
226 gsl_interp_free(data->log_e_of_log_h_interp);
227 gsl_interp_free(data->log_p_of_log_h_interp);
228 gsl_interp_free(data->log_h_of_log_p_interp);
229 gsl_interp_free(data->log_rho_of_log_h_interp);
230 gsl_interp_free(data->log_p_of_log_e_interp);
231 gsl_interp_free(data->log_p_of_log_rho_interp);
232 gsl_interp_free(data->log_cs2_of_log_h_interp);
233 gsl_interp_accel_free(data->log_e_of_log_p_acc);
234 gsl_interp_accel_free(data->log_e_of_log_h_acc);
235 gsl_interp_accel_free(data->log_p_of_log_h_acc);
236 gsl_interp_accel_free(data->log_h_of_log_p_acc);
237 gsl_interp_accel_free(data->log_rho_of_log_h_acc);
238 gsl_interp_accel_free(data->log_p_of_log_e_acc);
239 gsl_interp_accel_free(data->log_p_of_log_rho_acc);
240 gsl_interp_accel_free(data->log_cs2_of_log_h_acc);
241 LALFree(data->log_edat);
242 LALFree(data->log_pdat);
243 LALFree(data->mubdat);
244 LALFree(data->muedat);
245 LALFree(data->log_hdat);
246 LALFree(data->yedat);
247 LALFree(data->log_cs2dat);
248 LALFree(data->log_rhodat);
249 LALFree(data);
250 }
251 return;
252}
253
254static void eos_free_tabular(LALSimNeutronStarEOS * eos)
255{
256 if (eos) {
257 eos_free_tabular_data(eos->data.tabular);
258 LALFree(eos);
259 }
260 return;
261}
262
263/* Finding density where EOS becomes acausal */
264
265/* Evaluate vSound at each tabulated point until vSound>1 or you get to last
266 * point. If vSound>1 interpolate between that point and previous point to
267 * find h where EOS first becomes acausal. (You could do root-finding to be
268 * even more precise, but the calculation of v from the tabulated EOS isn't
269 * that exact anyway.) */
270
271/* Minimum pseudo-enthalpy at which EOS becomes acausal (speed of sound > 1).
272 * If the EOS is always causal, return some large value hmax instead. */
273static double eos_min_acausal_pseudo_enthalpy_tabular(double hmax,
274 LALSimNeutronStarEOS * eos)
275{
276 size_t i;
277 double h_im1, h_i;
278 double v_im1, v_i;
279 double m; /* slope for linear interpolation */
280 double hMinAcausal = hmax; /* default large number for EOS that is always causal */
281
282 h_im1 = exp(eos->data.tabular->log_hdat[0]);
283 v_im1 = eos_v_of_h_tabular(h_im1, eos);
284 for (i = 1; i < eos->data.tabular->ndat; i++) {
285 h_i = exp(eos->data.tabular->log_hdat[i]);
286 v_i = eos_v_of_h_tabular(h_i, eos);
287 if (v_i > 1.0) {
288 /* solve vsound(h) = 1 */
289 m = (v_i - v_im1) / (h_i - h_im1);
290 hMinAcausal = h_im1 + (1.0 - v_im1) / m;
291 break;
292 }
293 h_im1 = h_i;
294 v_im1 = v_i;
295 }
296// printf("hMinAcausal = %e, v = %e\n", hMinAcausal, eos_v_of_h_tabular(hMinAcausal, eos));
297
298 /* Value of h where EOS first becomes acausal.
299 * Or, if EOS is always causal, hmax */
300 return hMinAcausal;
301}
302
303static LALSimNeutronStarEOS *eos_alloc_tabular(double *nbdat, double *edat, double *pdat,
304 double *mubdat, double *muedat, double *hdat, double *yedat, double *cs2dat, size_t ndat, size_t ncol)
305{
306 LALSimNeutronStarEOS *eos;
307 LALSimNeutronStarEOSDataTabular *data;
308 size_t i;
309
310 eos = LALCalloc(1, sizeof(*eos));
311 data = LALCalloc(1, sizeof(*data));
312
313 eos->datatype = LALSIM_NEUTRON_STAR_EOS_DATA_TYPE_TABULAR;
314 eos->data.tabular = data;
315
316 /* setup function pointers */
317 eos->free = eos_free_tabular;
318 eos->e_of_p = eos_e_of_p_tabular;
319 eos->h_of_p = eos_h_of_p_tabular;
320 eos->e_of_h = eos_e_of_h_tabular;
321 eos->p_of_h = eos_p_of_h_tabular;
322 eos->rho_of_h = eos_rho_of_h_tabular;
323 eos->p_of_e = eos_p_of_e_tabular;
324 eos->p_of_rho = eos_p_of_rho_tabular;
325 eos->dedp_of_p = eos_dedp_of_p_tabular;
326 eos->v_of_h = eos_v_of_h_tabular;
327
328 data->log_rhodat = XLALMalloc(ndat * sizeof(*data->log_rhodat));
329
330 if(ncol == 2) {
331 /* allocate memory for eos data; ignore first points if 0 */
332 while (*pdat == 0.0 || *edat == 0.0) {
333 ++pdat;
334 ++edat;
335 --ndat;
336 }
337
338 data->ncol = ncol;
339 data->ndat = ndat;
340 data->log_pdat = XLALMalloc(ndat * sizeof(*data->log_pdat));
341 data->log_edat = XLALMalloc(ndat * sizeof(*data->log_edat));
342 data->log_hdat = XLALMalloc(ndat * sizeof(*data->log_hdat));
343
344 /* take log of eos data */
345 for (i = 0; i < ndat; ++i) {
346 data->log_pdat[i] = log(pdat[i]);
347 data->log_edat[i] = log(edat[i]);
348 }
349 /* compute pseudo-enthalpy h from dhdp */
350 /* Integrate in log space:
351 dhdp = 1 / [e(p) + p]
352 h(p) = h(p0) + \int_p0^p dhdp dp
353 h(p) = h(p0) + \int_ln(p0)^ln(p) exp[ln(p) + ln(dhdp)] dln(p)
354 First point is
355 h(p0) = p0 / [e(p0) + p0]
356 */
357 double *integrand;
358 integrand = LALMalloc(ndat * sizeof(*integrand));
359 for (i = 0; i < ndat; ++i)
360 integrand[i] = exp(data->log_pdat[i] + log(1.0 / (edat[i] + pdat[i])));
361
362 gsl_interp_accel * dhdp_of_p_acc_temp = gsl_interp_accel_alloc();
363 gsl_interp * dhdp_of_p_interp_temp = gsl_interp_alloc(gsl_interp_linear, ndat);
364 gsl_interp_init(dhdp_of_p_interp_temp, data->log_pdat, integrand, ndat);
365
366 data->log_hdat[0] = log(pdat[0] / (edat[0] + pdat[0]));
367 for (i = 1; i < ndat; ++i)
368 data->log_hdat[i] = log(exp(data->log_hdat[0]) + gsl_interp_eval_integ(dhdp_of_p_interp_temp, data->log_pdat, integrand, data->log_pdat[0], data->log_pdat[i], dhdp_of_p_acc_temp));
369
370 gsl_interp_free(dhdp_of_p_interp_temp);
371 gsl_interp_accel_free(dhdp_of_p_acc_temp);
372 LALFree(integrand);
373 }
374 else if (ncol > 2) {
375 /* allocate memory for eos data; ignore first points if 0 */
376 while (*pdat == 0.0 || *edat == 0.0 || *hdat == 0.0) {
377 ++pdat;
378 ++edat;
379 ++hdat;
380 --ndat;
381 }
382
383 data->ndat = ndat;
384 data->ncol = ncol - 1;
385 data->nbdat = XLALMalloc(ndat * sizeof(*data->nbdat));
386 data->log_pdat = XLALMalloc(ndat * sizeof(*data->log_pdat));
387 data->log_edat = XLALMalloc(ndat * sizeof(*data->log_edat));
388 data->mubdat = XLALMalloc(ndat * sizeof(*data->mubdat));
389 data->muedat = XLALMalloc(ndat * sizeof(*data->muedat));
390 data->log_hdat = XLALMalloc(ndat * sizeof(*data->log_hdat));
391 data->yedat = XLALMalloc(ndat * sizeof(*data->yedat));
392 data->log_cs2dat = XLALMalloc(ndat * sizeof(*data->log_cs2dat));
393
394 /* take log of eos data */
395 for (i = 0; i < ndat; ++i) {
396 data->nbdat[i] = nbdat[i];
397 data->log_pdat[i] = log(pdat[i]);
398 data->log_edat[i] = log(edat[i]);
399 data->mubdat[i] = mubdat[i];
400 data->muedat[i] = muedat[i];
401 data->log_hdat[i] = log(hdat[i]);
402 data->yedat[i] = yedat[i];
403 data->log_cs2dat[i] = log(cs2dat[i]);
404 }
405
406 /* these can be set up only using the new eos tables; so they are set up here */
407 data->log_cs2_of_log_h_acc = gsl_interp_accel_alloc();
408 data->log_cs2_of_log_h_interp = gsl_interp_alloc(gsl_interp_cspline, ndat);
409 gsl_interp_init(data->log_cs2_of_log_h_interp, data->log_hdat, data->log_cs2dat, ndat);
410 }
411
412 // Find rho from e, p, and h: rho = (e+p)/exp(h)
413 for (i = 0; i < ndat; i++)
414 data->log_rhodat[i] = log(edat[i] + pdat[i]) - exp(data->log_hdat[i]);
415
416 eos->pmax = exp(data->log_pdat[ndat - 1]);
417 eos->hmax = exp(data->log_hdat[ndat - 1]);
418
419 /* setup interpolation tables */
420
421 data->log_e_of_log_p_acc = gsl_interp_accel_alloc();
422 data->log_h_of_log_p_acc = gsl_interp_accel_alloc();
423 data->log_e_of_log_h_acc = gsl_interp_accel_alloc();
424 data->log_p_of_log_h_acc = gsl_interp_accel_alloc();
425 data->log_rho_of_log_h_acc = gsl_interp_accel_alloc();
426 data->log_p_of_log_e_acc = gsl_interp_accel_alloc();
427 data->log_p_of_log_rho_acc = gsl_interp_accel_alloc();
428
429 data->log_e_of_log_p_interp = gsl_interp_alloc(gsl_interp_cspline, ndat);
430 data->log_h_of_log_p_interp = gsl_interp_alloc(gsl_interp_cspline, ndat);
431 data->log_e_of_log_h_interp = gsl_interp_alloc(gsl_interp_cspline, ndat);
432 data->log_p_of_log_h_interp = gsl_interp_alloc(gsl_interp_cspline, ndat);
433 data->log_rho_of_log_h_interp = gsl_interp_alloc(gsl_interp_cspline, ndat);
434 data->log_p_of_log_e_interp = gsl_interp_alloc(gsl_interp_cspline, ndat);
435 data->log_p_of_log_rho_interp = gsl_interp_alloc(gsl_interp_cspline, ndat);
436
437 gsl_interp_init(data->log_e_of_log_p_interp, data->log_pdat, data->log_edat, ndat);
438 gsl_interp_init(data->log_h_of_log_p_interp, data->log_pdat, data->log_hdat, ndat);
439 gsl_interp_init(data->log_e_of_log_h_interp, data->log_hdat, data->log_edat, ndat);
440 gsl_interp_init(data->log_p_of_log_h_interp, data->log_hdat, data->log_pdat, ndat);
441 gsl_interp_init(data->log_rho_of_log_h_interp, data->log_hdat, data->log_rhodat, ndat);
442 gsl_interp_init(data->log_p_of_log_e_interp, data->log_edat, data->log_pdat, ndat);
443 gsl_interp_init(data->log_p_of_log_rho_interp, data->log_rhodat, data->log_pdat, ndat);
444
445 eos->hMinAcausal =
446 eos_min_acausal_pseudo_enthalpy_tabular(eos->hmax, eos);
447
448// printf("%e\n", XLALSimNeutronStarEOSEnergyDensityOfPressureGeometerized(eos->pmax, eos));
449//
450// printf("datatype = %d\n", eos->datatype);
451// printf("pmax = %e\n", eos->pmax);
452// printf("hmax = %e\n", eos->hmax);
453// printf("hMinAcausal = %e\n", eos->hMinAcausal);
454
455 return eos;
456}
457
458
459/** @endcond */
460
461/**
462 * @brief Reads a data file containing a tabulated equation of state.
463 * @details Read a data file specified by a path fname that contains two
464 * whitespace separated columns of equation of state data. The first column
465 * contains the pressure in Pa and the second column contains the energy
466 * density in J/m^3. Every line beginning with the character '#' then it is
467 * ignored. If the path is an absolute path then this specific file is opened;
468 * otherwise, search for the file in paths given in the environment variable
469 * LALSIM_DATA_PATH, and finally search in the installed PKG_DATA_DIR path.
470 * @param[in] fname The path of the file to open.
471 * @return A pointer to neutron star equation of state structure.
472 */
473LALSimNeutronStarEOS *XLALSimNeutronStarEOSFromFile(const char *fname)
474{
475 LALSimNeutronStarEOS *eos;
476 double *f_dat;
477 size_t ncol;
478 size_t ndat;
479 LALFILE *fp;
480
481 double *nbdat;
482 double *edat;
483 double *pdat;
484 double *mubdat;
485 double *muedat;
486 double *hdat;
487 double *yedat;
488 double *cs2dat;
489
490 fp = XLALSimReadDataFileOpen(fname);
491 if (!fp)
492 XLAL_ERROR_NULL(XLAL_EFUNC);
493
494 ndat = XLALSimReadDataFileNCol(&f_dat, &ncol, fp);
495 XLALFileClose(fp);
496
497 if (ndat == (size_t) (-1))
498 XLAL_ERROR_NULL(XLAL_EFUNC);
499
500 nbdat = LALMalloc(ndat * sizeof(*nbdat));
501 edat = LALMalloc(ndat * sizeof(*edat));
502 pdat = LALMalloc(ndat * sizeof(*pdat));
503 mubdat = LALMalloc(ndat * sizeof(*mubdat));
504 muedat = LALMalloc(ndat * sizeof(*muedat));
505 hdat = LALMalloc(ndat * sizeof(*hdat));
506 yedat = LALMalloc(ndat * sizeof(*yedat));
507 cs2dat = LALMalloc(ndat * sizeof(*cs2dat));
508
509
510 if (ncol > 2)
511 {
512 for (size_t i = 0 ; i < ndat ; i++) {
513 nbdat[i] = f_dat[i * ncol + 1];
514 edat[i] = f_dat[i * ncol + 2] * 1e3 * LAL_G_C2_SI; /* transform from CGS to SI and then to Geometrized units */
515 pdat[i] = f_dat[i * ncol + 3] * 1e-1 * LAL_G_C4_SI; /* transform from CGS to SI and then to Geometrized units */
516 mubdat[i] = f_dat[i * ncol + 4];
517 muedat[i] = f_dat[i * ncol + 5];
518 hdat[i] = f_dat[i * ncol + 6];
519 yedat[i] = f_dat[i * ncol + 7];
520 cs2dat[i] = f_dat[i * ncol + 8];
521 }
522 }
523 else if (ncol == 2)
524 {
525 for (size_t i = 0 ; i < ndat ; i++) {
526 pdat[i] = f_dat[i * ncol];
527 edat[i] = f_dat[i * ncol + 1];
528 }
529
530 nbdat = NULL;
531 mubdat = NULL;
532 muedat = NULL;
533 yedat = NULL;
534 cs2dat = NULL;
535 }
536 else if (ncol < 2)
537 {
538 fprintf(stderr, "error: equation of state files must have at least 2 columns, ncol >= 2\n");
539 exit(1);
540 }
541
542
543 eos = eos_alloc_tabular(nbdat, edat, pdat, mubdat, muedat, hdat, yedat, cs2dat, ndat, ncol);
544
545 XLALFree(f_dat);
546 LALFree(nbdat);
547 LALFree(edat);
548 LALFree(pdat);
549 LALFree(mubdat);
550 LALFree(muedat);
551 LALFree(hdat);
552 LALFree(yedat);
553 LALFree(cs2dat);
554
555 snprintf(eos->name, sizeof(eos->name), "%s", fname);
556 return eos;
557}
558
559/**
560 * @brief Creates an equation of state structure from tabulated equation
561 * of state data of a known name. The name of the tabulated equation of state
562 * must belong to the sample of equations of state from the old frame work or
563 * added for the new framework.
564 * @details A known, installed, named tabulated equation of state data file, whose name
565 * is included in the old EOS framework names or the new ones, is read and then used to
566 * create the equation of state structure.
567 * The equations of state for the OLD framework available are the representative sample drawn from
568 * http://xtreme.as.arizona.edu/NeutronStars/ they are:
569 * - ALF1
570 * - ALF2
571 * - ALF3
572 * - ALF4
573 * - AP1
574 * - AP2
575 * - AP3
576 * - AP4
577 * - APR4_EPP
578 * - BBB2
579 * - BGN1H1
580 * - BPAL12
581 * - BSK19
582 * - BSK20
583 * - BSK21
584 * - ENG
585 * - FPS
586 * - GNH3
587 * - GS1
588 * - GS2
589 * - H1
590 * - H2
591 * - H3
592 * - H4
593 * - H5
594 * - H6
595 * - H7
596 * - MPA1
597 * - MS1B
598 * - MS1B_PP
599 * - MS1_PP
600 * - MS1
601 * - MS2
602 * - PAL6
603 * - PCL2
604 * - PS
605 * - QMC700
606 * - SLY4
607 * - SLY
608 * - SQM1
609 * - SQM2
610 * - SQM3
611 * - WFF1
612 * - WFF2
613 * - WFF3
614 * We also include more modern equations from the CompOSE website
615 * https://compose.obspm.fr/ downloaded on 18 June 2018. These EOSs are:
616 * - APR
617 * - BHF_BBB2
618 * - KDE0V
619 * - KDE0V1
620 * - RS
621 * - SK255
622 * - SK272
623 * - SKA
624 * - SKB
625 * - SKI2
626 * - SKI3
627 * - SKI4
628 * - SKI5
629 * - SKI6
630 * - SKMP
631 * - SKOP
632 * - SLY2
633 * - SLY230A
634 * - SLY9
635 * And we include HQC18 from http://user.numazu-ct.ac.jp/~sumi/eos/HQC18_submit
636 * - HQC18
637 *
638 *
639 * @param[in] name The name of the equation of state.
640 * @return A pointer to neutron star equation of state structure.
641 */
642LALSimNeutronStarEOS *XLALSimNeutronStarEOSByName(const char *name)
643{
644 static const char fname_base[] = "LALSimNeutronStarEOS_";
645 static const char fname_extn[] = ".dat";
646
647 size_t n = XLAL_NUM_ELEM(lalSimNeutronStarEOSNames);
648 size_t i;
649 char fname[FILENAME_MAX];
650
651 for (i = 0; i < n ; ++i)
652 if (XLALStringCaseCompare(name, lalSimNeutronStarEOSNames[i]) == 0) {
653 LALSimNeutronStarEOS *eos;
654 snprintf(fname, sizeof(fname), "%s%s%s", fname_base, lalSimNeutronStarEOSNames[i],
655 fname_extn);
656 eos = XLALSimNeutronStarEOSFromFile(fname);
657 if (!eos)
658 XLAL_ERROR_NULL(XLAL_EFUNC);
659 snprintf(eos->name, sizeof(eos->name), "%s", lalSimNeutronStarEOSNames[i]);
660 return eos;
661 }
662
663 XLAL_PRINT_ERROR("Unrecognized EOS name %s...", name);
664 XLALPrintError("\tRecognised EOS names are: %s", lalSimNeutronStarEOSNames[0]);
665 for (i = 1; i < n ; ++i)
666 XLALPrintError(", %s", lalSimNeutronStarEOSNames[i]);
667 XLALPrintError("\n");
668 XLAL_ERROR_NULL(XLAL_ENAME);
669}
670
671/** @} */
672/** @} */
LALCalloc
#define LALCalloc(m, n)
LALMalloc
#define LALMalloc(n)
LALFree
#define LALFree(p)
name
const char * name
Definition: LALSimInspiralWaveformParams_common.c:6
LAL_G_C2_SI
#define LAL_G_C2_SI
Factor to convert density in kg/m^3 to geometerized units of m^-2.
Definition: LALSimNeutronStar.h:41
LALSimNeutronStarEOS
struct tagLALSimNeutronStarEOS LALSimNeutronStarEOS
Incomplete type for the neutron star Equation of State (EOS).
Definition: LALSimNeutronStar.h:56
LAL_G_C4_SI
#define LAL_G_C4_SI
Factor to convert pressure in Pa to geometerized units of m^-2.
Definition: LALSimNeutronStar.h:44
XLALSimReadDataFileOpen
LALFILE * XLALSimReadDataFileOpen(const char *fname)
Opens a specified data file, searching default path if necessary.
Definition: LALSimReadData.c:59
XLALSimReadDataFileNCol
size_t XLALSimReadDataFileNCol(double **data, size_t *ncol, LALFILE *fp)
Read a multi-column data file.
Definition: LALSimReadData.c:133
fprintf
#define fprintf
i
double i
Definition: bh_ringdown.c:118
e
double e
Definition: bh_ringdown.c:117
data
sigmaKerr data[0]
Definition: exact_derivatives-HrealHeader.c:4
XLALFileClose
int XLALFileClose(LALFILE *file)
XLAL_NUM_ELEM
#define XLAL_NUM_ELEM(x)
XLALMalloc
void * XLALMalloc(size_t n)
XLALFree
void XLALFree(void *p)
XLALSimNeutronStarEOSByName
LALSimNeutronStarEOS * XLALSimNeutronStarEOSByName(const char *name)
Creates an equation of state structure from tabulated equation of state data of a known name.
Definition: LALSimNeutronStarEOSTabular.c:642
lalSimNeutronStarEOSNames
const char *const lalSimNeutronStarEOSNames[111]
Recognised names of equations of state.
Definition: LALSimNeutronStarEOS.c:83
XLALSimNeutronStarEOSFromFile
LALSimNeutronStarEOS * XLALSimNeutronStarEOSFromFile(const char *fname)
Reads a data file containing a tabulated equation of state.
Definition: LALSimNeutronStarEOSTabular.c:473
XLALStringCaseCompare
int XLALStringCaseCompare(const char *s1, const char *s2)
m
static const INT4 m
XLAL_ERROR_NULL
#define XLAL_ERROR_NULL(...)
XLALPrintError
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
XLAL_PRINT_ERROR
#define XLAL_PRINT_ERROR(...)
XLAL_ENAME
XLAL_ENAME
XLAL_EFUNC
XLAL_EFUNC
p
p
fp
FILE * fp