|
1 | 1 | import numpy as np |
2 | 2 | import scipy |
| 3 | +from scipy.interpolate import interp1d, BSpline, splrep |
3 | 4 |
|
4 | 5 | def empirical_density_profile(rbins, pos, mass): |
5 | 6 | """ |
@@ -30,13 +31,16 @@ def empirical_density_profile(rbins, pos, mass): |
30 | 31 | V_shells = 4/3 * np.pi * (rbins[1:]**3 - rbins[:-1]**3) |
31 | 32 |
|
32 | 33 | # Compute density profile |
33 | | - density, _ = np.histogram(r_p, bins=rbins+1, weights=mass) |
| 34 | + density, _ = np.histogram(r_p, bins=rbins, weights=mass) |
34 | 35 | density /= V_shells |
35 | 36 |
|
36 | | - # Compute bin centers and return profile |
37 | | - #radius = 0.5 * (rbins[1:] + rbins[:-1]) |
| 37 | + #compute bin centers and return profile |
| 38 | + radius = 0.5 * (rbins[1:] + rbins[:-1]) |
38 | 39 |
|
39 | | - return density |
| 40 | + tck_s = splrep(radius, np.log10(density), s=0.15) |
| 41 | + dens2 = BSpline(*tck_s)(radius) |
| 42 | + |
| 43 | + return radius, 10**dens2 |
40 | 44 |
|
41 | 45 |
|
42 | 46 |
|
@@ -142,7 +146,10 @@ def makemodel(func, M, funcargs, rvals = 10.**np.linspace(-2.,4.,2000), pfile='' |
142 | 146 | R = np.nanmax(rvals) |
143 | 147 |
|
144 | 148 | # query out the density values |
145 | | - dvals = func(rvals,*funcargs) |
| 149 | + if func == empirical_density_profile: |
| 150 | + rvals, dvals = func(rvals,*funcargs) |
| 151 | + else: |
| 152 | + dvals = func(rvals,*funcargs) |
146 | 153 |
|
147 | 154 | # make the mass and potential arrays |
148 | 155 | mvals = np.zeros(dvals.size) |
|
0 commit comments