Coverage for bilby/gw/eos/eos.py: 64%

1import os 

2import numpy as np 

3from scipy.interpolate import interp1d, CubicSpline 

4 

5from .tov_solver import IntegrateTOV 

6from ...core import utils 

7 

8C_SI = utils.speed_of_light # m/s 

9C_CGS = C_SI * 100. 

10G_SI = utils.gravitational_constant # m^3 kg^-1 s^-2 

11MSUN_SI = utils.solar_mass # Kg 

12 

13# Stores conversions from geometerized to cgs or si unit systems 

14conversion_dict = {'pressure': {'cgs': C_SI ** 4. / G_SI * 10., 'si': C_SI ** 4. / G_SI, 'geom': 1.}, 

15 'energy_density': {'cgs': C_SI ** 4. / G_SI * 10., 'si': C_SI ** 4. / G_SI, 'geom': 1.}, 

16 'density': {'cgs': C_SI ** 2. / G_SI / 1000., 'si': C_SI ** 2. / G_SI, 'geom': 1.}, 

17 'pseudo_enthalpy': {'dimensionless': 1.}, 

18 'mass': {'g': C_SI ** 2. / G_SI * 1000, 'kg': C_SI ** 2. / G_SI, 'geom': 1., 

19 'm_sol': C_SI ** 2. / G_SI / MSUN_SI}, 

20 'radius': {'cm': 100., 'm': 1., 'km': .001}, 

21 'tidal_deformability': {'geom': 1.}} 

22 

23 

24# construct dictionary of pre-shipped EOS pressure denstity table 

25path_to_eos_tables = os.path.join(os.path.dirname(__file__), 'eos_tables') 

26list_of_eos_tables = os.listdir(path_to_eos_tables) 

27valid_eos_files = [i for i in list_of_eos_tables if 'LAL' in i] 

28valid_eos_file_paths = [os.path.join(path_to_eos_tables, filename) for filename in valid_eos_files] 

29valid_eos_names = [i.split('_', maxsplit=1)[-1].strip('.dat') for i in valid_eos_files] 

30valid_eos_dict = dict(zip(valid_eos_names, valid_eos_file_paths)) 

31 

32 

33class TabularEOS(object): 

34 """ 

35 Given a valid eos input format, such as 2-D array, an ascii file, or a string, parse, and interpolate 

36 

37 Parameters 

38 ========== 

39 eos: (numpy.ndarray, str, ASCII TABLE) 

40 if `numpy.ndarray` then user supplied pressure-density 2D numpy array. 

41 if `str` then given a valid eos name, relevant preshipped ASCII table will be loaded 

42 if ASCII TABLE then given viable file extensions, which include .txt,.dat, etc (np.loadtxt used), 

43 read in pressure density from file. 

44 sampling_flag: bool 

45 Do you plan on sampling the parameterized EOS? Highly recommended. Defaults to False. 

46 warning_flag: bool 

47 Keeps track of status of various physical checks on EoS. 

48 

49 Attributes 

50 ========== 

51 msg: str 

52 Human readable string describing the exception. 

53 code: int 

54 Exception error code. 

55 """ 

56 

57 def __init__(self, eos, sampling_flag=False, warning_flag=False): 

58 from scipy.integrate import cumulative_trapezoid 

59 

60 self.sampling_flag = sampling_flag 

61 self.warning_flag = warning_flag 

62 

63 if isinstance(eos, str): 

64 if eos in valid_eos_dict.keys(): 

65 table = np.loadtxt(valid_eos_dict[eos]) 

66 else: 

67 table = np.loadtxt(eos) 

68 elif isinstance(eos, np.ndarray): 

69 table = eos 

70 else: 

71 raise ValueError("eos provided is invalid type please supply a str name, str path to ASCII file, " 

72 "or a numpy array") 

73 

74 table = self.__remove_leading_zero(table) 

75 

76 # all have units of m^-2 

77 self.pressure = table[:, 0] 

78 self.energy_density = table[:, 1] 

79 

80 self.minimum_pressure = min(self.pressure) 

81 self.minimum_energy_density = min(self.energy_density) 

82 if (not self.check_monotonicity() and self.sampling_flag) or self.warning_flag: 

83 self.warning_flag = True 

84 else: 

85 integrand = self.pressure / (self.energy_density + self.pressure) 

86 self.pseudo_enthalpy = cumulative_trapezoid(integrand, np.log(self.pressure), initial=0) + integrand[0] 

87 

88 self.interp_energy_density_from_pressure = CubicSpline(np.log10(self.pressure), 

89 np.log10(self.energy_density), 

90 ) 

91 

92 self.interp_energy_density_from_pseudo_enthalpy = CubicSpline(np.log10(self.pseudo_enthalpy), 

93 np.log10(self.energy_density)) 

94 

95 self.interp_pressure_from_pseudo_enthalpy = CubicSpline(np.log10(self.pseudo_enthalpy), 

96 np.log10(self.pressure)) 

97 

98 self.interp_pseudo_enthalpy_from_energy_density = CubicSpline(np.log10(self.energy_density), 

99 np.log10(self.pseudo_enthalpy)) 

100 

101 self.__construct_all_tables() 

102 

103 self.minimum_pseudo_enthalpy = min(self.pseudo_enthalpy) 

104 if not self.check_causality() and self.sampling_flag: 

105 self.warning_flag = True 

106 

107 def __remove_leading_zero(self, table): 

108 """ 

109 For interpolation of lalsimulation tables; 

110 loglog interpolation breaks if the first entries are 0s 

111 """ 

112 

113 if table[0, 0] == 0. or table[0, 1] == 0.: 

114 return table[1:, :] 

115 

116 else: 

117 return table 

118 

119 def energy_from_pressure(self, pressure, interp_type='CubicSpline'): 

120 """ 

121 Find value of energy_from_pressure 

122 as in lalsimulation, return e = K * p**(3./5.) below min pressure 

123 

124 Parameters 

125 ========== 

126 pressure: float 

127 pressure in geometerized units. 

128 interp_type: str 

129 String specifying which interpolation type to use. 

130 Currently implemented: 'CubicSpline', 'linear'. 

131 energy_density: float 

132 energy-density in geometerized units. 

133 """ 

134 pressure = np.atleast_1d(pressure) 

135 energy_returned = np.zeros(pressure.size) 

136 indices_less_than_min = np.nonzero(pressure < self.minimum_pressure) 

137 indices_greater_than_min = np.nonzero(pressure >= self.minimum_pressure) 

138 

139 # We do this special for less than min pressure 

140 energy_returned[indices_less_than_min] = 10 ** (np.log10(self.energy_density[0]) + 

141 (3. / 5.) * (np.log10(pressure[indices_less_than_min]) - 

142 np.log10(self.pressure[0]))) 

143 

144 if interp_type == 'CubicSpline': 

145 energy_returned[indices_greater_than_min] = ( 

146 10. ** self.interp_energy_density_from_pressure(np.log10(pressure[indices_greater_than_min]))) 

147 elif interp_type == 'linear': 

148 energy_returned[indices_greater_than_min] = np.interp(pressure[indices_greater_than_min], 

149 self.pressure, 

150 self.energy_density) 

151 else: 

152 raise ValueError('Interpolation scheme must be linear or CubicSpline') 

153 

154 if energy_returned.size == 1: 

155 return energy_returned[0] 

156 else: 

157 return energy_returned 

158 

159 def pressure_from_pseudo_enthalpy(self, pseudo_enthalpy, interp_type='CubicSpline'): 

160 """ 

161 Find p(h) 

162 as in lalsimulation, return p = K * h**(5./2.) below min enthalpy 

163 

164 :param pseudo_enthalpy (`float`): Dimensionless pseudo-enthalpy. 

165 :interp_type (`str`): String specifying interpolation type. 

166 Current implementations are 'CubicSpline', 'linear'. 

167 

168 :return pressure (`float`): pressure in geometerized units. 

169 """ 

170 pseudo_enthalpy = np.atleast_1d(pseudo_enthalpy) 

171 pressure_returned = np.zeros(pseudo_enthalpy.size) 

172 indices_less_than_min = np.nonzero(pseudo_enthalpy < self.minimum_pseudo_enthalpy) 

173 indices_greater_than_min = np.nonzero(pseudo_enthalpy >= self.minimum_pseudo_enthalpy) 

174 

175 pressure_returned[indices_less_than_min] = 10. ** (np.log10(self.pressure[0]) + 

176 2.5 * (np.log10(pseudo_enthalpy[indices_less_than_min]) - 

177 np.log10(self.pseudo_enthalpy[0]))) 

178 

179 if interp_type == 'CubicSpline': 

180 pressure_returned[indices_greater_than_min] = ( 

181 10. ** self.interp_pressure_from_pseudo_enthalpy(np.log10(pseudo_enthalpy[indices_greater_than_min]))) 

182 elif interp_type == 'linear': 

183 pressure_returned[indices_greater_than_min] = np.interp(pseudo_enthalpy[indices_greater_than_min], 

184 self.pseudo_enthalpy, 

185 self.pressure) 

186 else: 

187 raise ValueError('Interpolation scheme must be linear or CubicSpline') 

188 

189 if pressure_returned.size == 1: 

190 return pressure_returned[0] 

191 else: 

192 return pressure_returned 

193 

194 def energy_density_from_pseudo_enthalpy(self, pseudo_enthalpy, interp_type='CubicSpline'): 

195 """ 

196 Find energy_density_from_pseudo_enthalpy(pseudo_enthalpy) 

197 as in lalsimulation, return e = K * h**(3./2.) below min enthalpy 

198 

199 :param pseudo_enthalpy (`float`): Dimensionless pseudo-enthalpy. 

200 :param interp_type (`str`): String specifying interpolation type. 

201 Current implementations are 'CubicSpline', 'linear'. 

202 

203 :return energy_density (`float`): energy-density in geometerized units. 

204 """ 

205 pseudo_enthalpy = np.atleast_1d(pseudo_enthalpy) 

206 energy_returned = np.zeros(pseudo_enthalpy.size) 

207 indices_less_than_min = np.nonzero(pseudo_enthalpy < self.minimum_pseudo_enthalpy) 

208 indices_greater_than_min = np.nonzero(pseudo_enthalpy >= self.minimum_pseudo_enthalpy) 

209 

210 energy_returned[indices_less_than_min] = 10 ** (np.log10(self.energy_density[0]) + 

211 1.5 * (np.log10(pseudo_enthalpy[indices_less_than_min]) - 

212 np.log10(self.pseudo_enthalpy[0]))) 

213 if interp_type == 'CubicSpline': 

214 x = np.log10(pseudo_enthalpy[indices_greater_than_min]) 

215 energy_returned[indices_greater_than_min] = 10 ** self.interp_energy_density_from_pseudo_enthalpy(x) 

216 elif interp_type == 'linear': 

217 energy_returned[indices_greater_than_min] = np.interp(pseudo_enthalpy[indices_greater_than_min], 

218 self.pseudo_enthalpy, 

219 self.energy_density) 

220 else: 

221 raise ValueError('Interpolation scheme must be linear or CubicSpline') 

222 

223 if energy_returned.size == 1: 

224 return energy_returned[0] 

225 else: 

226 return energy_returned 

227 

228 def pseudo_enthalpy_from_energy_density(self, energy_density, interp_type='CubicSpline'): 

229 """ 

230 Find h(epsilon) 

231 as in lalsimulation, return h = K * e**(2./3.) below min enthalpy 

232 

233 :param energy_density (`float`): energy-density in geometerized units. 

234 :param interp_type (`str`): String specifying interpolation type. 

235 Current implementations are 'CubicSpline', 'linear'. 

236 

237 :return pseudo_enthalpy (`float`): Dimensionless pseudo-enthalpy. 

238 """ 

239 energy_density = np.atleast_1d(energy_density) 

240 pseudo_enthalpy_returned = np.zeros(energy_density.size) 

241 indices_less_than_min = np.nonzero(energy_density < self.minimum_energy_density) 

242 indices_greater_than_min = np.nonzero(energy_density >= self.minimum_energy_density) 

243 

244 pseudo_enthalpy_returned[indices_less_than_min] = 10 ** (np.log10(self.pseudo_enthalpy[0]) + (2. / 3.) * 

245 (np.log10(energy_density[indices_less_than_min]) - 

246 np.log10(self.energy_density[0]))) 

247 

248 if interp_type == 'CubicSpline': 

249 x = np.log10(energy_density[indices_greater_than_min]) 

250 pseudo_enthalpy_returned[indices_greater_than_min] = 10**self.interp_pseudo_enthalpy_from_energy_density(x) 

251 elif interp_type == 'linear': 

252 pseudo_enthalpy_returned[indices_greater_than_min] = np.interp(energy_density[indices_greater_than_min], 

253 self.energy_density, 

254 self.pseudo_enthalpy) 

255 else: 

256 raise ValueError('Interpolation scheme must be linear or CubicSpline') 

257 

258 if pseudo_enthalpy_returned.size == 1: 

259 return pseudo_enthalpy_returned[0] 

260 else: 

261 return pseudo_enthalpy_returned 

262 

263 def dedh(self, pseudo_enthalpy, rel_dh=1e-5, interp_type='CubicSpline'): 

264 """ 

265 Value of [depsilon/dh](p) 

266 

267 :param pseudo_enthalpy (`float`): Dimensionless pseudo-enthalpy. 

268 :param interp_type (`str`): String specifying interpolation type. 

269 Current implementations are 'CubicSpline', 'linear'. 

270 :param rel_dh (`float`): Relative step size in pseudo-enthalpy space. 

271 

272 :return dedh (`float`): Derivative of energy-density with respect to pseudo-enthalpy 

273 evaluated at `pseudo_enthalpy` in geometerized units. 

274 """ 

275 

276 # step size=fraction of value 

277 dh = pseudo_enthalpy * rel_dh 

278 

279 eps_upper = self.energy_density_from_pseudo_enthalpy(pseudo_enthalpy + dh, interp_type=interp_type) 

280 eps_lower = self.energy_density_from_pseudo_enthalpy(pseudo_enthalpy - dh, interp_type=interp_type) 

281 

282 return (eps_upper - eps_lower) / (2. * dh) 

283 

284 def dedp(self, pressure, rel_dp=1e-5, interp_type='CubicSpline'): 

285 """ 

286 Find value of [depsilon/dp](p) 

287 

288 :param pressure (`float`): pressure in geometerized units. 

289 :param rel_dp (`float`): Relative step size in pressure space. 

290 :param interp_type (`float`): String specifying interpolation type. 

291 Current implementations are 'CubicSpline', 'linear'. 

292 

293 :return dedp (`float`): Derivative of energy-density with respect to pressure 

294 evaluated at `pressure`. 

295 """ 

296 

297 # step size=fraction of value 

298 dp = pressure * rel_dp 

299 

300 eps_upper = self.energy_from_pressure(pressure + dp, interp_type=interp_type) 

301 eps_lower = self.energy_from_pressure(pressure - dp, interp_type=interp_type) 

302 

303 return (eps_upper - eps_lower) / (2. * dp) 

304 

305 def velocity_from_pseudo_enthalpy(self, pseudo_enthalpy, interp_type='CubicSpline'): 

306 """ 

307 Returns the speed of sound in geometerized units in the 

308 neutron star at the specified pressure. 

309 

310 Assumes the equation 

311 vs = c (de/dp)^{-1/2} 

312 

313 :param pseudo_enthalpy (`float`): Dimensionless pseudo-enthalpy. 

314 :param interp_type (`str`): String specifying interpolation type. 

315 Current implementations are 'CubicSpline', 'linear'. 

316 

317 :return v_s (`float`): Speed of sound at `pseudo-enthalpy` in geometerized units. 

318 """ 

319 pressure = self.pressure_from_pseudo_enthalpy(pseudo_enthalpy, interp_type=interp_type) 

320 return self.dedp(pressure, interp_type=interp_type) ** -0.5 

321 

322 def check_causality(self): 

323 """ 

324 Checks to see if the equation of state is causal i.e. the speed 

325 of sound in the star is less than the speed of light. 

326 Returns True if causal, False if not. 

327 """ 

328 pmax = self.pressure[-1] 

329 emax = self.energy_from_pressure(pmax) 

330 hmax = self.pseudo_enthalpy_from_energy_density(emax) 

331 vsmax = self.velocity_from_pseudo_enthalpy(hmax) 

332 if vsmax < 1.1: 

333 return True 

334 else: 

335 return False 

336 

337 def check_monotonicity(self): 

338 """ 

339 Checks to see if the equation of state is monotonically increasing 

340 in energy density-pressure space. Returns True if monotonic, False if not. 

341 """ 

342 e1 = self.energy_density[1:] 

343 e2 = self.energy_density[:-1] 

344 ediff = e1 - e2 

345 e_negatives = len(np.where(ediff < 0)) 

346 

347 p1 = self.pressure[1:] 

348 p2 = self.pressure[:-1] 

349 pdiff = p1 - p2 

350 p_negatives = len(np.where(pdiff < 0)) 

351 if e_negatives > 1 or p_negatives > 1: 

352 return False 

353 else: 

354 return True 

355 

356 def __get_plot_range(self, data): 

357 """ 

358 Determines default plot range based on data provided. 

359 """ 

360 low = np.amin(data) 

361 high = np.amax(data) 

362 dx = 0.05 * (high - low) 

363 

364 xmin = low - dx 

365 xmax = high + dx 

366 xlim = [xmin, xmax] 

367 

368 return xlim 

369 

370 def __construct_all_tables(self): 

371 """Pressure and epsilon already tabular, now create array of enthalpies""" 

372 edat = self.energy_density 

373 hdat = [self.pseudo_enthalpy_from_energy_density(e) for e in edat] 

374 self.pseudo_enthalpy = np.array(hdat) 

375 

376 def plot(self, rep, xlim=None, ylim=None, units=None): 

377 """ 

378 Given a representation in the form 'energy_density-pressure', plot the EoS in that space. 

379 

380 Parameters 

381 ========== 

382 rep: str 

383 Representation to plot. For example, plotting in energy_density-pressure space, 

384 specify 'energy_density-pressure' 

385 xlim: list 

386 Plotting bounds for x-axis in the form [low, high]. 

387 Defaults to 'None' which will plot from 10% below min x value to 10% above max x value 

388 ylim: list 

389 Plotting bounds for y-axis in the form [low, high]. 

390 Defaults to 'None' which will plot from 10% below min y value to 10% above max y value 

391 units: str 

392 Specifies unit system to plot. Currently can plot in CGS:'cgs', SI:'si', or geometerized:'geom' 

393 

394 Returns 

395 ======= 

396 fig: matplotlib.figure.Figure 

397 EOS plot. 

398 """ 

399 import matplotlib.pyplot as plt 

400 

401 # Set data based on specified representation 

402 varnames = rep.split('-') 

403 

404 assert varnames[0] != varnames[ 

405 1], 'Cannot plot the same variable against itself. Please choose another representation' 

406 

407 # Correspondence of rep parameter, data, and latex symbol 

408 # rep_dict = {'energy_density': [self.epsilon, r'$\epsilon$'], 

409 # 'pressure': [self.p, r'$p$'], 'pseudo_enthalpy': [pseudo_enthalpy, r'$h$']} 

410 

411 # FIXME: The second element in these arrays should be tex labels, but tex's not working rn 

412 rep_dict = {'energy_density': [self.energy_density, 'energy_density'], 

413 'pressure': [self.pressure, 'pressure'], 

414 'pseudo_enthalpy': [self.pseudo_enthalpy, 'pseudo_enthalpy']} 

415 

416 xname = varnames[1] 

417 yname = varnames[0] 

418 

419 # Set units 

420 eos_default_units = {'pressure': 'cgs', 'energy_density': 'cgs', 'density': 'cgs', 

421 'pseudo_enthalpy': 'dimensionless'} 

422 if units is None: 

423 units = [eos_default_units[yname], eos_default_units[xname]] # Default unit system is cgs 

424 elif isinstance(units, str): 

425 units = [units, units] # If only one unit system given, use for both 

426 

427 xunits = units[1] 

428 yunits = units[0] 

429 

430 # Ensure valid units 

431 if xunits not in list(conversion_dict[xname].keys()) or yunits not in list(conversion_dict[yname].keys()): 

432 s = ''' 

433 Invalid unit system. Valid variable-unit pairs are: 

434 p: {p_units} 

435 e: {e_units} 

436 rho: {rho_units} 

437 h: {h_units}. 

438 '''.format(p_units=list(conversion_dict['pressure'].keys()), 

439 e_units=list(conversion_dict['energy_density'].keys()), 

440 rho_units=list(conversion_dict['density'].keys()), 

441 h_units=list(conversion_dict['pseudo_enthalpy'].keys())) 

442 raise ValueError(s) 

443 

444 xdat = rep_dict[xname][0] * conversion_dict[xname][xunits] 

445 ydat = rep_dict[yname][0] * conversion_dict[yname][yunits] 

446 

447 xlabel = rep_dict[varnames[1]][1].replace('_', ' ') 

448 ylabel = rep_dict[varnames[0]][1].replace('_', ' ') + '(' + xlabel + ')' 

449 

450 # Determine plot ranges. Currently shows 10% wider than actual data range. 

451 if xlim is None: 

452 xlim = self.__get_plot_range(xdat) 

453 

454 if ylim is None: 

455 ylim = self.__get_plot_range(ydat) 

456 

457 fig, ax = plt.subplots() 

458 

459 ax.loglog(xdat, ydat) 

460 ax.set_xlim(xlim) 

461 ax.set_ylim(ylim) 

462 ax.set_xlabel(xlabel) 

463 ax.set_ylabel(ylabel) 

464 

465 return fig 

466 

467 

468def spectral_adiabatic_index(gammas, x): 

469 arg = 0 

470 for i in range(len(gammas)): 

471 arg += gammas[i] * x ** i 

472 

473 return np.exp(arg) 

474 

475 

476class SpectralDecompositionEOS(TabularEOS): 

477 """ 

478 Parameterized EOS using a spectral 

479 decomposition per Lindblom 

480 arXiv: 1009.0738v2. Inherits from TabularEOS. 

481 

482 Parameters 

483 ========== 

484 gammas: list 

485 List of adiabatic expansion parameters used 

486 to construct the equation of state in various 

487 spaces. 

488 p0: float 

489 The starting point in pressure of the high-density EoS. This is stitched to 

490 the low-density portion of the SLY EoS model. The default value chosen is set to 

491 a sufficiently low pressure so that the high-density EoS will never be 

492 overconstrained. 

493 e0/c**2: float 

494 The starting point in energy-density of the high-density EoS. This is stitched to 

495 the low-density portion of the SLY EoS model. The default value chosen is set to 

496 a sufficiently low energy density so that the high-density EoS will never be 

497 overconstrained. 

498 xmax: float 

499 highest dimensionless pressure value in EoS 

500 npts: float (optional) 

501 number of points in pressure-energy density data. 

502 """ 

503 

504 def __init__(self, gammas, p0=3.01e33, e0=2.03e14, xmax=None, npts=100, sampling_flag=False, warning_flag=False): 

505 self.warning_flag = warning_flag 

506 self.gammas = gammas 

507 self.p0 = p0 

508 self.e0 = e0 

509 self.xmax = xmax 

510 self.npts = npts 

511 if self.xmax is None: 

512 self.xmax = self.__determine_xmax() 

513 self.sampling_flag = sampling_flag 

514 self.__construct_a_of_x_table() 

515 

516 # Construct pressure-energy density table and 

517 # set up interpolating functions. 

518 

519 if self.warning_flag and self.sampling_flag: 

520 # If sampling prior is enabled and adiabatic check 

521 # has failed, empty the array values 

522 self.e_pdat = np.zeros((2, 2)) 

523 else: 

524 self.e_pdat = self.__construct_e_of_p_table() 

525 

526 super().__init__(self.e_pdat, sampling_flag=self.sampling_flag, warning_flag=self.warning_flag) 

527 

528 def __determine_xmax(self, a_max=6.0): 

529 highest_order_gamma = np.abs(self.gammas[-1])[0] 

530 expansion_order = float(len(self.gammas) - 1) 

531 

532 xmax = (np.log(a_max) / highest_order_gamma) ** (1. / expansion_order) 

533 return xmax 

534 

535 def __mu_integrand(self, x): 

536 

537 return 1. / spectral_adiabatic_index(self.gammas, x) 

538 

539 def mu(self, x): 

540 from scipy.integrate import quad 

541 return np.exp(-quad(self.__mu_integrand, 0, x)[0]) 

542 

543 def __eps_integrand(self, x): 

544 

545 return np.exp(x) * self.mu(x) / spectral_adiabatic_index(self.gammas, x) 

546 

547 def energy_density(self, x, eps0): 

548 from scipy.integrate import quad 

549 quad_result, quad_err = quad(self.__eps_integrand, 0, x) 

550 eps_of_x = (eps0 * C_CGS ** 2.) / self.mu(x) + self.p0 / self.mu(x) * quad_result 

551 return eps_of_x 

552 

553 def __construct_a_of_x_table(self): 

554 

555 xdat = np.linspace(0, self.xmax, num=self.npts) 

556 

557 # Generate adiabatic index points until a point is out of prior range 

558 if self.sampling_flag: 

559 adat = np.empty(self.npts) 

560 for i in range(self.npts): 

561 if 0.6 < spectral_adiabatic_index(self.gammas, xdat[i]) < 4.6: 

562 adat[i] = spectral_adiabatic_index(self.gammas, xdat[i]) 

563 else: 

564 break 

565 

566 # Truncate arrays to last point within range 

567 adat = adat[:i] 

568 xmax_new = xdat[i - 1] 

569 

570 # If EOS is too short, set prior to 0, else resample the function and set new xmax 

571 if xmax_new < 4. or i == 0: 

572 self.warning_flag = True 

573 else: 

574 xdat = np.linspace(0, xmax_new, num=self.npts) 

575 adat = spectral_adiabatic_index(self.gammas, xdat) 

576 self.xmax = xmax_new 

577 else: 

578 adat = spectral_adiabatic_index(self.gammas, xdat) 

579 

580 self.adat = adat 

581 

582 def __construct_e_of_p_table(self): 

583 

584 """ 

585 Creates p, epsilon table for a given set of spectral parameters 

586 """ 

587 

588 # make p range 

589 # to match lalsimulation tables: array = [pressure, density] 

590 x_range = np.linspace(0, self.xmax, self.npts) 

591 p_range = self.p0 * np.exp(x_range) 

592 

593 eos_vals = np.zeros((self.npts, 2)) 

594 eos_vals[:, 0] = p_range 

595 

596 for i in range(0, len(x_range)): 

597 eos_vals[i, 1] = self.energy_density(x_range[i], self.e0) 

598 

599 # convert eos to geometrized units in *m^-2* 

600 # IMPORTANT 

601 eos_vals = eos_vals * 0.1 * G_SI / C_SI ** 4 

602 

603 # doing as those before me have done and using SLY4 as low density region 

604 # SLY4 in geometrized units 

605 low_density_path = os.path.join(os.path.dirname(__file__), 'eos_tables', 'LALSimNeutronStarEOS_SLY4.dat') 

606 low_density = np.loadtxt(low_density_path) 

607 

608 cutoff = eos_vals[0, :] 

609 

610 # Then find overlap point 

611 break_pt = len(low_density) 

612 for i in range(1, len(low_density)): 

613 if low_density[-i, 0] < cutoff[0] and low_density[-i, 1] < cutoff[1]: 

614 break_pt = len(low_density) - i + 1 

615 break 

616 

617 # stack EOS arrays 

618 eos_vals = np.vstack((low_density[0:break_pt, :], eos_vals)) 

619 

620 return eos_vals 

621 

622 

623class EOSFamily(object): 

624 """ 

625 Create a EOS family and get mass-radius information 

626 

627 Parameters 

628 ========== 

629 eos: object 

630 Supply a `TabularEOS` class (or subclass) 

631 npts: float 

632 Number of points to calculate for mass-radius relation. Default is 500. 

633 

634 Notes 

635 ===== 

636 The mass-radius and mass-k2 data should be 

637 populated here via the TOV solver upon object construction. 

638 """ 

639 def __init__(self, eos, npts=500): 

640 from scipy.optimize import minimize_scalar 

641 self.eos = eos 

642 

643 # FIXME: starting_energy_density is set somewhat arbitrarily 

644 starting_energy_density = 1.6e-10 

645 ending_energy_density = max(self.eos.energy_density) 

646 log_starting_energy_density = np.log(starting_energy_density) 

647 log_ending_energy_density = np.log(ending_energy_density) 

648 log_energy_density_grid = np.linspace(log_starting_energy_density, 

649 log_ending_energy_density, 

650 num=npts) 

651 energy_density_grid = np.exp(log_energy_density_grid) 

652 

653 # Generate m, r, and k2 lists 

654 mass = [] 

655 radius = [] 

656 k2love_number = [] 

657 for i in range(len(energy_density_grid)): 

658 tov_solver = IntegrateTOV(self.eos, energy_density_grid[i]) 

659 m, r, k2 = tov_solver.integrate_TOV() 

660 mass.append(m) 

661 radius.append(r) 

662 k2love_number.append(k2) 

663 

664 # Check if maximum mass has been found 

665 if i > 0 and mass[i] <= mass[i - 1]: 

666 break 

667 

668 # If we're not at the end of the array, determine actual maximum mass. Else, assume 

669 # last point is the maximum mass and proceed. 

670 if i < (npts - 1): 

671 # Now replace with point with interpolated maximum mass 

672 # This is done by interpolating the last three points and then 

673 # minimizing the negative of the interpolated function 

674 x = [energy_density_grid[i - 2], energy_density_grid[i - 1], energy_density_grid[i]] 

675 y = [mass[i - 2], mass[i - 1], mass[i]] 

676 

677 f = interp1d(x, y, kind='quadratic', bounds_error=False, fill_value='extrapolate') 

678 

679 res = minimize_scalar(lambda x: -f(x)) 

680 

681 # Integrate max energy density to get maximum mass 

682 tov_solver = IntegrateTOV(self.eos, res.x) 

683 mfin, rfin, k2fin = tov_solver.integrate_TOV() 

684 

685 mass[-1] = mfin 

686 radius[-1] = rfin 

687 k2love_number[-1] = k2fin 

688 

689 # Currently, everything is in geometerized units. 

690 # The mass variables have dimensions of length, k2 is dimensionless 

691 # and radii have dimensions of length. Calculate dimensionless lambda 

692 # with these quantities, then convert to SI. 

693 

694 # Calculating dimensionless lambda values from k2, radii, and mass 

695 tidal_deformability = [2. / 3. * k2 * r ** 5. / m ** 5. for k2, r, m in 

696 zip(k2love_number, radius, mass)] 

697 

698 # As a last resort, if highest mass is still smaller than second 

699 # to last point, remove the last point from each array 

700 if mass[-1] < mass[-2]: 

701 mass = mass[:-1] 

702 radius = radius[:-1] 

703 k2love_number = k2love_number[:-1] 

704 tidal_deformability = tidal_deformability[:-1] 

705 

706 self.mass = np.array(mass) 

707 self.radius = np.array(radius) 

708 self.k2love_number = np.array(k2love_number) 

709 self.tidal_deformability = np.array(tidal_deformability) 

710 self.maximum_mass = mass[-1] * conversion_dict['mass']['m_sol'] 

711 

712 def radius_from_mass(self, m): 

713 """ 

714 :param m: mass of neutron star in solar masses 

715 :return: radius of neutron star in meters 

716 """ 

717 f = CubicSpline(self.mass, self.radius, bc_type='natural', extrapolate=True) 

718 

719 mass_converted_to_geom = m * MSUN_SI * G_SI / C_SI ** 2. 

720 return f(mass_converted_to_geom) 

721 

722 def k2_from_mass(self, m): 

723 """ 

724 :param m: mass of neutron star in solar masses. 

725 :return: dimensionless second tidal love number. 

726 """ 

727 f = CubicSpline(self.mass, self.k2love_number, bc_type='natural', extrapolate=True) 

728 

729 m_geom = m * MSUN_SI * G_SI / C_SI ** 2. 

730 return f(m_geom) 

731 

732 def lambda_from_mass(self, m): 

733 """ 

734 Convert from equation of state model parameters to 

735 component tidal parameters. 

736 

737 :param m: Mass of neutron star in solar masses. 

738 :return: Tidal parameter of neutron star of mass m. 

739 """ 

740 

741 # Get lambda from mass and equation of state 

742 

743 r = self.radius_from_mass(m) 

744 k = self.k2_from_mass(m) 

745 m_geom = m * MSUN_SI * G_SI / C_SI ** 2. 

746 c = m_geom / r 

747 lmbda = (2. / 3.) * k / c ** 5. 

748 

749 return lmbda 

750 

751 def __get_plot_range(self, data): 

752 low = np.amin(data) 

753 high = np.amax(data) 

754 dx = 0.05 * (high - low) 

755 

756 xmin = low - dx 

757 xmax = high + dx 

758 xlim = [xmin, xmax] 

759 

760 return xlim 

761 

762 def plot(self, rep, xlim=None, ylim=None, units=None): 

763 """ 

764 Given a representation in the form 'm-r', plot the family in that space. 

765 

766 Parameters 

767 ========== 

768 rep: str 

769 Representation to plot. For example, plotting in mass-radius space, specify 'm-r' 

770 xlim: list 

771 Plotting bounds for x-axis in the form [low, high]. 

772 Defaults to 'None' which will plot from 10% below min x value to 10% above max x value 

773 ylim: list 

774 Plotting bounds for y-axis in the form [low, high]. 

775 Defaults to 'None' which will plot from 10% below min y value to 10% above max y value 

776 units: str 

777 Specifies unit system to plot. Currently can plot in CGS:'cgs', SI:'si', or geometerized:'geom' 

778 

779 Returns 

780 ======= 

781 fig: matplotlib.figure.Figure 

782 EOS Family plot. 

783 """ 

784 import matplotlib.pyplot as plt 

785 

786 # Set data based on specified representation 

787 varnames = rep.split('-') 

788 

789 assert varnames[0] != varnames[ 

790 1], 'Cannot plot the same variable against itself. Please choose another representation' 

791 

792 # Correspondence of rep parameter, data, and latex symbol 

793 rep_dict = {'mass': [self.mass, r'$M$'], 'radius': [self.radius, r'$R$'], 'k2': [self.k2love_number, r'$k_2$'], 

794 'tidal_deformability': [self.tidal_deformability, r'$l$']} 

795 

796 xname = varnames[1] 

797 yname = varnames[0] 

798 

799 # Set units 

800 fam_default_units = {'mass': 'm_sol', 'radius': 'km', 'tidal_deformability': 'geom'} 

801 if units is None: 

802 units = [fam_default_units[yname], fam_default_units[xname]] # Default unit system is cgs 

803 elif isinstance(units, str): 

804 units = [units, units] # If only one unit system given, use for both