diff --git a/EXPtools/basis_builder/makemodel.py b/EXPtools/basis_builder/makemodel.py index 9d6aa01..cf5f3fc 100644 --- a/EXPtools/basis_builder/makemodel.py +++ b/EXPtools/basis_builder/makemodel.py @@ -1,5 +1,6 @@ import numpy as np import scipy +from scipy.interpolate import interp1d, BSpline, splrep def empirical_density_profile(rbins, pos, mass): """ @@ -30,13 +31,16 @@ def empirical_density_profile(rbins, pos, mass): V_shells = 4/3 * np.pi * (rbins[1:]**3 - rbins[:-1]**3) # Compute density profile - density, _ = np.histogram(r_p, bins=rbins+1, weights=mass) + density, _ = np.histogram(r_p, bins=rbins, weights=mass) density /= V_shells - # Compute bin centers and return profile - #radius = 0.5 * (rbins[1:] + rbins[:-1]) + #compute bin centers and return profile + radius = 0.5 * (rbins[1:] + rbins[:-1]) - return density + tck_s = splrep(radius, np.log10(density), s=0.15) + dens2 = BSpline(*tck_s)(radius) + + return radius, 10**dens2 @@ -142,7 +146,10 @@ def makemodel(func, M, funcargs, rvals = 10.**np.linspace(-2.,4.,2000), pfile='' R = np.nanmax(rvals) # query out the density values - dvals = func(rvals,*funcargs) + if func == empirical_density_profile: + rvals, dvals = func(rvals,*funcargs) + else: + dvals = func(rvals,*funcargs) # make the mass and potential arrays mvals = np.zeros(dvals.size) diff --git a/EXPtools/utils/coefficients.py b/EXPtools/utils/coefficients.py index d44ad9d..ee1593d 100644 --- a/EXPtools/utils/coefficients.py +++ b/EXPtools/utils/coefficients.py @@ -18,7 +18,7 @@ def remove_terms(original_coefficients, n, l, m): for i in range(len(l)): lm_idx = int(l[i]*(l[i]+1) / 2) + m[i] print(n[i],l[i],m[i], lm_idx) - try: coefs_matrix[n[i], lm_idx, t] = np.complex128(0) + try: coefs_matrix[lm_idx, n[i], t] = np.complex128(0) except IndexError: continue copy_coeffcients.setMatrix(mat=coefs_matrix[:,:, t], time=t_snaps[t])