Skip to content

Commit 075f4c0

Browse files
ChunHuangPhyChunHuangPhy
authored
Merge pull request #132 from ChunHuangPhy/maximum_rho
Fix names and definitions
2 parents 8454de0 + bdbcece commit 075f4c0

File tree

5 files changed

+212
-238
lines changed

5 files changed

+212
-238
lines changed

‎.gitignore

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,5 @@
11
__pycache__/
22
.DS_Store
33
downloads/
4-
output/
4+
output/
5+
CompactObject_TOV.egg-info/

‎TOVsolver/maximum_central_density.py

Lines changed: 13 additions & 54 deletions
Original file line numberDiff line numberDiff line change
@@ -3,50 +3,9 @@
33
from TOVsolver import unit
44
from TOVsolver.main import OutputMR
55

6-
def maxium_central_density(energy_density, pressure, central_densitys=np.logspace(14.3, 15.6, 5) * unit.g_cm_3, num2=30):
7-
"""Outputs the maxium central density of a stable EoS
8-
Args:
9-
energy_density (numpy 1Darray): Density of EoS
10-
pressure (numpy 1Darray): Pressure of EoS
11-
central_densitys (numpy 1Darray): The range of central density
12-
num2 (int): The number of segments in the density interval of second search. At least 5 points
13-
Returns:
14-
density (float): The maxium central density, in unit.g_cm_3
15-
"""
16-
############## Below is the main part of two searches for peaks
17-
######## First search
18-
Ms = [-1,-2] # Meaningless initialization, only to ensure that the following 'if (i>0) and (Ms [-1]<=Ms [-2])' statements can run properly
19-
store_d_range = [central_densitys[-2], central_densitys[-1]] # When the following loop does not output a result, initialization here will have its meaning
20-
# Find the extremum point within the predetermined range of central density and return the central density near the extremum point
21-
for i, rho in enumerate(central_densitys):
22-
M, R = OutputMR('', energy_density, pressure, [rho])[0]
23-
Ms.append(M)
24-
if (i>0) and (Ms[-1] <= Ms[-2]):
25-
store_d_range = [central_densitys[i-2], central_densitys[i]]
26-
break
27-
Ms_larg = Ms[-1] # Used to store the mass corresponding to the maximum density during the first peak search
28-
29-
######## Second search
30-
Ms = [-1,-2] # Reinitialize
31-
# Update and refine the central density range, and ultimately return the central density of the extremum points
32-
store_d_range = np.geomspace(store_d_range[0], store_d_range[1], num2) # At least 5 points
33-
# Note that the first and last points have already been calculated, so there is no need to traverse them
34-
for i, rho in enumerate(store_d_range[1:-1]):
35-
M, R = OutputMR('', energy_density, pressure, [rho])[0]
36-
Ms.append(M)
37-
if Ms[-1] <= Ms[-2]: # Note that due to the second peak search refining the interval, the result is generally not obtained at the first point.
38-
# Therefore, initializing Ms=[-1, -2] is acceptable
39-
return store_d_range[1:][i-1]
40-
# When the above peak search fails, continue comparing the last point
41-
if Ms_larg < Ms[-1]:
42-
return store_d_range[-2]
43-
else:
44-
return store_d_range[-1]
45-
46-
47-
def maxium_central_density_2(energy_density, pressure,
6+
def maximum_central_density(energy_density, pressure,
487
rho_lo=10**14.3 * unit.g_cm_3, rho_up=10**15.6 * unit.g_cm_3,
49-
tol_for_mass=1e-2, tol_for_rho=1e-5, threshold_R=40*unit.km):
8+
tol_for_mass=1e-3, tol_for_rho=1e-7, threshold_R=40*unit.km):
509
"""
5110
Find the center density of a stable point using dichotomy approach
5211
Args:
@@ -58,7 +17,7 @@ def maxium_central_density_2(energy_density, pressure,
5817
tol_for_rho (float): Relative tolerance of the density
5918
threshold_R (float): Threshold of Radius. Once the minmum Radius of the star is lager than the value of threshold, it means the results are very unphysical
6019
Returns:
61-
density (float): The maxium central density, in unit.g_cm_3
20+
density (float): The maxium central density, in g_cm_3
6221
mass_large (float): The maximum mass, in mass of sun
6322
radius_small (float): The corresponding radius of the maximum mass, in km
6423
"""
@@ -125,7 +84,7 @@ def maxium_central_density_2(energy_density, pressure,
12584
# If R_min > threshold_R, then the results are definitely not physical.
12685
if R_min>threshold_R:
12786
# print('R',R_min/unit.km)
128-
return rhos[Catch], Ms_larg/unit.Msun, R_min/unit.km # density in natural unit, mass in sun
87+
return rhos[Catch]/unit.g_cm_3, Ms_larg/unit.Msun, R_min/unit.km # density in g_cm3, mass in sun, radius in km
12988

13089
###### Start looping code *****
13190
# Update and refine the central density range, and ultimately return the central density of the maximum points
@@ -142,7 +101,7 @@ def maxium_central_density_2(energy_density, pressure,
142101
rho = rhos[1]
143102
# Note that the useful accumulation amounts here are 0 and 2 points, and the point to be solved is in the middle, so we take the value of Cumul_rho2s
144103
result = OutputMR('', energy_density, pressure, [rho])[0]
145-
Ms.insert(1, result["M"]) # Insert points into the Ms of the original 3 points
104+
Ms.insert(1, result[0]) # Insert points into the Ms of the original 3 points
146105

147106
# Firstly, we store the data, and secondly, we judge the data
148107
if Ms[1] >= Ms[2]:
@@ -153,7 +112,7 @@ def maxium_central_density_2(energy_density, pressure,
153112
# Note that the useful accumulation amounts here are 0, 1, and 2 points, and the point to be solved is on the right side, so we take Cumul_rho2s [: 3]
154113
result = OutputMR('', energy_density, pressure, [rho])[0]
155114
# Since the refinement calculation is carried out on the basis of the previous stage, if d2_min is continuously performed, The operation of d2_max=query_nearest-values (* paras_for_strink_0) is actually cumbersome
156-
Ms.insert(3, result["M"]) # Insert a point into Ms with 4 points
115+
Ms.insert(3, result[0]) # Insert a point into Ms with 4 points
157116

158117
if Ms[3] >= Ms[2]:
159118
Catch = 3
@@ -166,13 +125,13 @@ def maxium_central_density_2(energy_density, pressure,
166125
rho = rhos[1]
167126
# Note that the useful accumulation amounts here are 0 and 4 points, and the point to be calculated is on the left side of the middle, so we take Cumul_rho2s
168127
result = OutputMR('', energy_density, pressure, [rho])[0]
169-
Ms.insert(1, result["M"])
128+
Ms.insert(1, result[0])
170129

171130
elif i == 2:
172131
rho = rhos[2]
173132
# Note that the useful accumulation amounts here are 0 and 1 points, and the point to be solved is on the right side, so we take Cumul_rho2s [: 2]
174133
result = OutputMR('', energy_density, pressure, [rho])[0]
175-
Ms.insert(2, result["M"])
134+
Ms.insert(2, result[0])
176135

177136
if Ms[1] >= Ms[2]:
178137
Peak == True
@@ -182,7 +141,7 @@ def maxium_central_density_2(energy_density, pressure,
182141
rho = rhos[3]
183142
# Note that the useful accumulation amounts here are 0, 1, and 2 points, and the point to be solved is on the right side, so we take Cumul_rho2s [: 3]
184143
result = OutputMR('', energy_density, pressure, [rho])[0]
185-
Ms.insert(3, result["M"])
144+
Ms.insert(3, result[0])
186145

187146
if Ms[2] >= Ms[3]:
188147
Peak == True
@@ -204,13 +163,13 @@ def maxium_central_density_2(energy_density, pressure,
204163
rho = rhos[1]
205164
# Note that the useful accumulation amounts here are 0 and 4 points, and the point to be calculated is on the left side of the middle, so we take Cumul_rho2s
206165
result = OutputMR('', energy_density, pressure, [rho])[0]
207-
Ms.insert(1, result["M"])
166+
Ms.insert(1, result[0])
208167

209168
elif i == 2:
210169
rho = rhos[2]
211170
# Note that the useful accumulation amounts here are 0 and 1 points, and the point to be solved is on the right side, so we take Cumul_rho2s [: 2]
212171
result = OutputMR('', energy_density, pressure, [rho])[0]
213-
Ms.insert(2, result["M"])
172+
Ms.insert(2, result[0])
214173

215174
if (Ms[1] >= Ms[2]) and (Ms[1] > Ms[0]):
216175
Peak == True
@@ -220,7 +179,7 @@ def maxium_central_density_2(energy_density, pressure,
220179
rho = rhos[3]
221180
# Note that the useful accumulation amounts here are 0, 1, and 2 points, and the point to be solved is on the right side, so we take Cumul_rho2s [: 3]
222181
result = OutputMR('', energy_density, pressure, [rho])[0]
223-
Ms.insert(3, result["M"])
182+
Ms.insert(3, result[0])
224183

225184
if (Ms[2] >= Ms[3]) and (Ms[2] > Ms[1]):
226185
Peak == True
@@ -252,4 +211,4 @@ def maxium_central_density_2(energy_density, pressure,
252211
# print('Peak:',Peak)
253212
# print('Ms',np.array(Ms)/unit.Msun)
254213

255-
return rhos[Catch], Ms_larg/unit.Msun, R_min/unit.km
214+
return rhos[Catch]/unit.g_cm_3, Ms_larg/unit.Msun, R_min/unit.km

‎TOVsolver/solver_code.py

Lines changed: 15 additions & 35 deletions
Original file line numberDiff line numberDiff line change
@@ -11,11 +11,10 @@ def m1_from_mc_m2(mc, m2):
1111
1212
Args:
1313
mc (float): chrip mass of a GW event, unit in solar mass.
14-
m2 (float or numpy array): the determined mass for one of the star, this
15-
is computed from sampling of EoS.
14+
m2 (float or numpy array): the determined mass for one of the star, this is computed from sampling of EoS.
1615
1716
Returns:
18-
m1 (float or numpy array): the companion star mass in solar mass.
17+
m1 (array): the companion star mass in solar mass.
1918
"""
2019
m2 = np.array(m2)
2120
num1 = (2.0 / 3.0) ** (1.0 / 3.0) * mc**5.0
@@ -88,17 +87,13 @@ def TOV_def(r, y, inveos, ad_index):
8887
8988
Args:
9089
r (float): raius as integrate varible
91-
y (psudo-varible): containing pressure, mass, h and b as intergarte varibles
92-
to solve out the TOV equation
93-
inveos: the invert of the eos, pressure and energy density relation to integrate
94-
and interpolate.
90+
y (psudo-varible): containing pressure, mass, h and b as intergarte varibles to solve out the TOV equation
91+
inveos: the invert of the eos, pressure and energy density relation to integrate and interpolate.
9592
9693
Returns:
97-
Mass (array): The array that contains all the Stars' masses, in M_sun as a
98-
Units.
94+
Mass (array): The array that contains all the Stars' masses, in M_sun as a Units.
9995
Radius (array): The array that contains all the Stars's radius, in km.
100-
Tidal Deformability (array): The array that contains correpsonding Tidal property,
101-
These are dimension-less.
96+
Tidal Deformability (array): The array that contains correpsonding Tidal property, These are dimension-less.
10297
"""
10398
pres, m, h, b = y
10499

@@ -155,27 +150,17 @@ def tidal_deformability(y2, Mns, Rns):
155150

156151
# Function solves the TOV equation, returning mass and radius
157152
def solveTOV_tidal(center_rho, energy_density, pressure):
158-
"""Solve TOV equation from given Equation of state in the neutron star
159-
core density range
153+
"""Solve TOV equation from given Equation of state in the neutron star core density range
160154
161155
Args:
162-
center_rho(array): This is the energy density here is fixed in main
163-
that is np.logspace(14.3, 15.6, 50)
164-
energy_density (array): Desity array of the neutron star EoS, in MeV/fm^{-3}
165-
Notice here for simiplicity, we omitted G/c**4 magnitude, so
166-
(value in MeV/fm^{-3})*G/c**4, could convert to the energy density we are
167-
using, please check the Test_EOS.csv to double check the order of magnitude.
168-
169-
pressure (array): Pressure array of neutron star EoS, also in nautral unit
170-
with MeV/fm^{-3}, still please check the Test_EOS.csv, the conversion is
171-
(value in dyn/cm3)*G/c**4.
156+
center_rho(array): This is the energy density here is fixed in main that is np.logspace(14.3, 15.6, 50)
157+
energy_density (array): Desity array of the neutron star EoS, in MeV/fm^{-3}. Notice here for simiplicity, we omitted G/c**4 magnitude, so (value in MeV/fm^{-3})*G/c**4, could convert to the energy density we are using, please check the Test_EOS.csv to double check the order of magnitude.
158+
pressure (array): Pressure array of neutron star EoS, also in nautral unit with MeV/fm^{-3}, still please check the Test_EOS.csv, the conversion is (value in dyn/cm3)*G/c**4.
172159
173160
Returns:
174-
Mass (array): The array that contains all the Stars' masses, in M_sun as a
175-
Units.
161+
Mass (array): The array that contains all the Stars' masses, in M_sun as a Units.
176162
Radius (array): The array that contains all the Stars's radius, in km.
177-
Tidal Deformability (array): The array that contains correpsonding Tidal property,
178-
These are dimension-less.
163+
Tidal Deformability (array): The array that contains correpsonding Tidal property, These are dimension-less.
179164
"""
180165
# eos = UnivariateSpline(np.log10(eps), np.log10(pres), k=1, s=0)
181166
# inveos = UnivariateSpline(np.log10(pressure), np.log10(energy_density), k=1, s=0)
@@ -242,21 +227,16 @@ def solveTOV_tidal(center_rho, energy_density, pressure):
242227

243228

244229
def solveTOV(center_rho, Pmin, eos, inveos):
245-
"""Solve TOV equation from given Equation of state in the neutron star
246-
core density range
230+
"""Solve TOV equation from given Equation of state in the neutron star core density range
247231
248232
Args:
249-
center_rho(array): This is the energy density here is fixed in main
250-
that is range [10**14.3, 10**15.6] * unit.g_cm_3
251-
233+
center_rho(float): The energy density in unit unit.g_cm_3
252234
Pmin (float): In unit.G / unit.c**4
253-
254235
eos (function): pressure vs. energy density, energy density in unit.G / unit.c**2, pressure in unit.G / unit.c**4
255-
256236
inveos (function): energy density vs. pressure
257237
258238
Returns:
259-
Mass (float): The Stars' masses
239+
Mass (float): The Stars' mass
260240
Radius (float): The Stars's radius
261241
"""
262242
r = 4.441e-16 * unit.cm

‎Test_Case/test_Inference_polytrope.ipynb

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -24,7 +24,7 @@
2424
"import ultranest.stepsampler\n",
2525
"\n",
2626
"import EOSgenerators.Polytrope_EOS as Polytrope\n",
27-
"from TOVsolver.maximum_central_density import maxium_central_density\n",
27+
"from TOVsolver.maximum_central_density import maximum_central_density\n",
2828
"from TOVsolver.solver_code import solveTOV\n",
2929
"from TOVsolver.unit import g_cm_3, dyn_cm_2, km, Msun, e0, MeV, fm\n",
3030
"import TOVsolver.main as main"
@@ -416,7 +416,7 @@
416416
"name": "python",
417417
"nbconvert_exporter": "python",
418418
"pygments_lexer": "ipython3",
419-
"version": "3.7.4"
419+
"version": "3.12.4"
420420
}
421421
},
422422
"nbformat": 4,

0 commit comments

Comments
 (0)