Update makemodel.py

Author: NGC

EXPtools/basis_builder/makemodel.py

@@ -2,125 +2,6 @@
 import scipy 
 from scipy.interpolate import interp1d, BSpline, splrep
 
-def empirical_density_profile(rbins, pos, mass):
-    """
-    Computes the number density radial profile assuming all particles have the same mass.
-
-    Args:
-        pos (ndarray): array of particle positions in cartesian coordinates with shape (n,3).
-        mass (ndarray): array of particle masses with shape (n,).
-        nbins (int, optional): number of bins in the radial profile. Default is 500.
-        rmin (float, optional): minimum radius of the radial profile. Default is 0.
-        rmax (float, optional): maximum radius of the radial profile. Default is 600.
-        log_space (bool, optional): whether to use logarithmic binning. Default is False.
-
-    Returns:
-        tuple: a tuple containing the arrays of radius and density with shapes (nbins,) and (nbins,), respectively.
-
-    Raises:
-        ValueError: if pos and mass arrays have different lengths or if nbins is not a positive integer.
-    """
-    if len(pos) != len(mass):
-        raise ValueError("pos and mass arrays must have the same length")
-    if not isinstance(len(rbins), int) or len(rbins) <= 0:
-        raise ValueError("rbins must be a positive integer")
-
-    # Compute radial distances
-    r_p = np.sqrt(np.sum(pos**2, axis=1))
-
-    V_shells = 4/3 * np.pi * (rbins[1:]**3 - rbins[:-1]**3)
-
-    # Compute density profile
-    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])
-
-    tck_s = splrep(radius, np.log10(density), s=0.15)
-    dens2 = BSpline(*tck_s)(radius)
-
-    return radius, 10**dens2
-
-
-
-def powerhalo(r, rs=1.,rc=0.,alpha=1.,beta=1.e-7):
-    """return generic twopower law distribution
-    
-    inputs
-    ----------
-    r      : (float) radius values
-    rs     : (float, default=1.) scale radius 
-    rc     : (float, default=0. i.e. no core) core radius
-    alpha  : (float, default=1.) inner halo slope
-    beta   : (float, default=1.e-7) outer halo slope
-    
-    returns
-    ----------
-    densities evaluated at r
-    
-    notes
-    ----------
-    different combinations are known distributions.
-    alpha=1,beta=2 is NFW
-    alpha=1,beta=3 is Hernquist
-    alpha=2.5,beta=0 is a typical single power law halo
-    
-    
-    """
-    ra = r/rs
-    return 1./(((ra+rc/rs)**alpha)*((1+ra)**beta))
-    
- 
-def powerhalorolloff(r, rs=1., rc=0., alpha=1., beta=1.e-7):
-    """return generic twopower law distribution with an erf rolloff
-    
-    inputs
-    ----------
-    r      : (float) radius values
-    rs     : (float, default=1.) scale radius 
-    rc     : (float, default=0. i.e. no core) core radius
-    alpha  : (float, default=1.) inner halo slope
-    beta   : (float, default=1.e-7) outer halo slope
-    
-    returns
-    ----------
-    densities evaluated at r
-    
-    notes
-    ----------
-    different combinations are known distributions.
-    alpha=1,beta=2 is NFW
-    alpha=1,beta=3 is Hernquist
-    alpha=2.5,beta=0 is a typical single power law halo
-    
-    
-    """
-    ra = r/rs
-    dens = 1./(((ra+rc/rs)**alpha)*((1+ra)**beta))
-    rtrunc = 25*rs
-    wtrunc = rtrunc*0.2
-    rolloff = 0.5 - 0.5*scipy.special.erf((r-rtrunc)/wtrunc)
-    return dens*rolloff
-
-
-def plummer_density(radius, scale_radius=1.0, mass=1.0, astronomicalG=1.0):
-    """basic plummer density profile"""
-    return ((3.0*mass)/(4*np.pi))*(scale_radius**2.)*((scale_radius**2 + radius**2)**(-2.5))
-
-def twopower_density_withrolloff(r, a, alpha, beta, rcen, wcen):
-    """a twopower density profile"""
-    ra = r/a
-    prefac = 0.5*(1.-scipy.special.erf((ra-rcen)/wcen))
-    return prefac*(ra**-alpha)*(1+ra)**(-beta+alpha)
-
-def hernquist_halo(r, a):
-    return 1 / ( 2*np.pi * (r/a) * (1 + r/a)**3)
-
-
-
-## Teporary Makemodel function that will be integrated into pyEXP
-
 def makemodel(func, M, funcargs, rvals = 10.**np.linspace(-2.,4.,2000), pfile='', plabel = '',verbose=True):
     """make an EXP-compatible spherical basis function table