Coverage for bilby/gw/eos/tov_solver.py: 100%

1# Monica Rizzo, 2019 

2 

3import numpy as np 

4 

5 

6class IntegrateTOV: 

7 """Class that given an initial pressure a mass radius value and a k2-love number 

8 """ 

9 

10 def __init__(self, eos, eps_0): 

11 self.eos = eos 

12 # determine central values 

13 pseudo_enthalpy0 = self.eos.pseudo_enthalpy_from_energy_density(eps_0) 

14 

15 # evolve first step analytically 

16 self.pseudo_enthalpy = pseudo_enthalpy0 - 1e-3 * pseudo_enthalpy0 

17 

18 mass0, radius0 = self.__mass_radius_cent(pseudo_enthalpy0, self.pseudo_enthalpy) 

19 

20 # k2 integration starting vals 

21 H0 = radius0 ** 2 

22 B0 = 2. * radius0 

23 

24 self.y = np.array([mass0, radius0, H0, B0]) 

25 

26 def __mass_radius_cent(self, pseudo_enthalpy0, pseudo_enthalpy): 

27 """ 

28 Calculate approximate mass/radius central values (for starting integration) 

29 

30 Eqns 7 + 8 from http://www.ccom.ucsd.edu/~lindblom/Publications/50_1992ApJ...398..569L.pdf 

31 """ 

32 

33 eps_c = self.eos.energy_density_from_pseudo_enthalpy(pseudo_enthalpy0) 

34 p_c = self.eos.pressure_from_pseudo_enthalpy(pseudo_enthalpy0) 

35 depsdh_c = self.eos.dedh(pseudo_enthalpy0) 

36 

37 radius = ((3. * (pseudo_enthalpy0 - pseudo_enthalpy)) / (2. * np.pi * (eps_c + 3. * p_c))) ** 0.5 \ 

38 * (1. - 0.25 * (eps_c - 3. * p_c - (3. / 5.) * depsdh_c) * 

39 ((pseudo_enthalpy0 - pseudo_enthalpy) / (eps_c + 3. * p_c))) 

40 

41 mass = (4. * np.pi) / 3. * eps_c * radius ** 3 * (1. - (3. / 5.) * 

42 depsdh_c * (pseudo_enthalpy0 - pseudo_enthalpy) / eps_c) 

43 

44 return mass, radius 

45 

46 def __tov_eqns(self, h, y): 

47 """ 

48 Taken from http://www.ccom.ucsd.edu/~lindblom/Publications/50_1992ApJ...398..569L.pdf 

49 

50 TOV equations in terms of pseudo-enthalpy; evolves in mass and radius 

51 H and B equations taken from: 

52 

53 https://journals.aps.org/prd/pdf/10.1103/PhysRevD.81.123016 

54 

55 Transformed (hopefully to evolve in H(h) and B(h) 

56 """ 

57 

58 m = y[0] 

59 r = y[1] 

60 H = y[2] 

61 B = y[3] 

62 eps = self.eos.energy_density_from_pseudo_enthalpy(h, interp_type='CubicSpline') 

63 p = self.eos.pressure_from_pseudo_enthalpy(h, interp_type='CubicSpline') 

64 depsdp = self.eos.dedp(p, interp_type='CubicSpline') 

65 

66 dmdh = (- (4. * np.pi * eps * r ** 3 * (r - 2. * m)) / 

67 (m + 4. * np.pi * r ** 3 * p)) 

68 

69 drdh = -(r * (r - 2. * m)) / (m + 4. * np.pi * r ** 3 * p) 

70 

71 dHdh = B * drdh 

72 

73 # taken from Damour & Nagar 

74 e_lam = (1. - 2. * m / r) ** (-1) 

75 

76 C1 = 2. / r + e_lam * (2. * m / r ** 2. + 4. * np.pi * r * (p - eps)) 

77 C0 = (e_lam * (- 6. / r ** 2. + 4. * np.pi * (eps + p) * 

78 depsdp + 4. * np.pi * (5. * eps + 9. * p)) - 

79 (2. * (m + 4. * np.pi * r ** 3. * p) / 

80 (r ** 2. - 2. * m * r)) ** 2.) 

81 

82 dBdh = -(C1 * B + C0 * H) * drdh 

83 

84 y_dot = np.array([dmdh, drdh, dHdh, dBdh]) 

85 

86 return y_dot 

87 

88 def __calc_k2(self, R, Beta, H, C): 

89 """ 

90 Using evolved quantities, calculate the k2 love 

91 

92 https://journals.aps.org/prd/pdf/10.1103/PhysRevD.81.123016 

93 

94 https://journals.aps.org/prd/pdf/10.1103/PhysRevD.80.084035 

95 """ 

96 

97 y = (R * Beta) / H 

98 

99 num = ((8. / 5.) * (1. - 2. * C) ** 2 * 

100 C ** 5 * (2. * C * (y - 1.) - y + 2.)) 

101 denom = (2. * C * (4. * (y + 1.) * C ** 4 + (6. * y - 4.) * C ** 3 + 

102 (26. - 22. * y) * C ** 2 + 3. * (5. * y - 8.) * 

103 C - 3. * y + 6.) - 3. * (1. - 2 * C) ** 2 * 

104 (2. * C * (y - 1.) - y + 2.) * 

105 np.log(1. / (1. - 2. * C))) 

106 

107 return num / denom 

108 

109 def integrate_TOV(self): 

110 """ 

111 Evolves TOV+k2 equations and returns final quantities 

112 """ 

113 from scipy.integrate import solve_ivp 

114 

115 # integration settings the same as in lalsimulation 

116 rel_err = 1e-4 

117 abs_err = 0.0 

118 

119 result = solve_ivp(self.__tov_eqns, (self.pseudo_enthalpy, 1e-16), self.y, rtol=rel_err, 

120 atol=abs_err) 

121 m_fin = result.y[0, -1] 

122 r_fin = result.y[1, -1] 

123 H_fin = result.y[2, -1] 

124 B_fin = result.y[3, -1] 

125 

126 k_2 = self.__calc_k2(r_fin, B_fin, H_fin, m_fin / r_fin) 

127 

128 return m_fin, r_fin, k_2