LALSimulation  5.4.0.1-fe68b98
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. */
35 struct 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 
66 static 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 
82 static 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 
98 static 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 
114 static 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 
130 static 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 
146 static 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 
162 static 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 
178 static 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 
195 static 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 
222 static 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 
254 static 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. */
273 static 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 
303 static 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  */
473 LALSimNeutronStarEOS *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  */
642 LALSimNeutronStarEOS *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
XLALSimReadDataFileNCol
size_t XLALSimReadDataFileNCol(double **data, size_t *ncol, LALFILE *fp)
Read a multi-column data file.
Definition: LALSimReadData.c:133
XLALSimReadDataFileOpen
LALFILE * XLALSimReadDataFileOpen(const char *fname)
Opens a specified data file, searching default path if necessary.
Definition: LALSimReadData.c:59
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)
lalSimNeutronStarEOSNames
const char *const lalSimNeutronStarEOSNames[111]
Recognised names of equations of state.
Definition: LALSimNeutronStarEOS.c:83
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
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
list p
fp
FILE * fp