converts 3d density into 2d density parameter
        :param amp:
        :param sigma_x:
        :param sigma_y:
        :return:
498994
def _amp2d_to_3d(self, amp, sigma_x, sigma_y):
        
        return amp / (np.sqrt(np.pi) * np.sqrt(sigma_x * sigma_y * 2))
1198059
def MakeGaussianPmf(mu, sigma, num_sigmas, n=201):
    
    pmf = Pmf()
    low = mu - num_sigmas * sigma
    high = mu + num_sigmas * sigma

    for x in numpy.linspace(low, high, n):
        p = EvalGaussianPdf(x, mu, sigma)
        pmf.Set(x, p)
    pmf.Normalize()
    return pmf
1046152
def _gaussian_membership_sigma(smoothness, eps = 0.0005): # 275us @ smothness=10
    
    """
    error = 0
    deltas = [0.1, 0.01, 0.001, 0.0001]
    sigma = smoothness * 0.3
    point = -1. * (smoothness + 0.5)
    for delta in deltas:
        while error < eps:
            sigma += delta
            error = scipy.stats.norm.cdf(0.5, point, sigma) - scipy.stats.norm.cdf(-0.5, point, sigma) # x, mu, sigma
        sigma -= delta
    return sigma
700889
def ae_core_density(self):
        """"""
        mesh, values, attrib = self._parse_radfunc("ae_core_density")
        return RadialFunction(mesh, values)
607121
def gauss_pdf(x, mu, sigma):
    """"""
    return 1 / np.sqrt(2 * np.pi) / sigma * np.exp(-(x - mu) ** 2 / 2. / sigma ** 2)
245326
def _calc_density(
    x: np.ndarray,
    y: np.ndarray,
):
        

    # Calculate the point density
    xy = np.vstack([x,y])
    z = gaussian_kde(xy)(xy)

    min_z = np.min(z)
    max_z = np.max(z)

    # Scale between 0 and 1
    scaled_z = (z-min_z)/(max_z-min_z)

    return(scaled_z)
480154
def derivatives(self, x, y, amp, sigma, center_x=0, center_y=0):
        
        x_ = x - center_x
        y_ = y - center_y
        R = np.sqrt(x_**2 + y_**2)
        sigma_x, sigma_y = sigma, sigma
        if isinstance(R, int) or isinstance(R, float):
            R = max(R, self.ds)
        else:
            R[R <= self.ds] = self.ds
        alpha = self.alpha_abs(R, amp, sigma)
        return alpha / R * x_, alpha / R * y_
340892
def function(self, x, y, amp, sigma, center_x=0, center_y=0):
        
        x_ = x - center_x
        y_ = y - center_y
        r = np.sqrt(x_**2 + y_**2)
        sigma_x, sigma_y = sigma, sigma
        c = 1. / (2 * sigma_x * sigma_y)
        num_int = self._num_integral(r, c)
        amp_density = self._amp2d_to_3d(amp, sigma_x, sigma_y)
        amp2d = amp_density / (np.sqrt(np.pi) * np.sqrt(sigma_x * sigma_y * 2))
        amp2d *= 2 * 1. / (2 * c)
        return num_int * amp2d
256431
def d2gauss(x,m_y,sigma):
    ''''''

    return gauss(x,m_y,sigma)*[-1/sigma**2 + (n-m_y)**2/sigma**4 for n in x]
322526
def mass_enclosed_2d(self, r, kwargs_profile):
        
        kwargs = copy.deepcopy(kwargs_profile)
        try:
            del kwargs['center_x']
            del kwargs['center_y']
        except:
            pass
        # integral of self.density_2d(x)* 2*np.pi * x *dx, 0, r
        out = integrate.quad(lambda x: self.density_2d(x, kwargs)*2*np.pi*x, 0, r)
        return out[0]
315114
def sigma_r2(self, r, kwargs_profile, kwargs_anisotropy, kwargs_light):
        
        if self._mass_profile == 'power_law':
            if self._anisotropy_type == 'r_ani':
                if self._light_profile == 'Hernquist':
                    sigma_r = self.power_law_anisotropy(r, kwargs_profile, kwargs_anisotropy, kwargs_light)
                else:
                    raise ValueError('light profile %s not supported for Jeans solver' % self._light_profile)
            else:
                raise ValueError('anisotropy type %s not implemented in Jeans equation modelling' % self._anisotropy_type)
        else:
            raise ValueError('mass profile type %s not implemented in Jeans solver' % self._mass_profile)
        return sigma_r
621622
def _M2(Xvar, Yvar, mask_X=None, mask_Y=None, xsum=0, xconst=0, ysum=0, yconst=0, weights=None, diag_only=False):
    
    if mask_X is None and mask_Y is None:
        return _M2_dense(Xvar, Yvar, weights=weights, diag_only=diag_only)
    else:
        # Check if one of the masks is not None, modify it and also adjust the constant columns:
        if mask_X is None:
            mask_X = np.ones(Xvar.shape[1], dtype=np.bool)
            xconst = np.ones(0, dtype=float)
        if mask_Y is None:
            mask_Y = np.ones(Yvar.shape[1], dtype=np.bool)
            yconst = np.ones(0, dtype=float)
    if _is_zero(xsum) and _is_zero(ysum) or _is_zero(xconst) and _is_zero(yconst):
        return _M2_sparse(Xvar, mask_X, Yvar, mask_Y, weights=weights)
    else:
        return _M2_const(Xvar, mask_X, xsum[mask_X], xconst, Yvar, mask_Y, ysum[mask_Y], yconst, weights=weights)
588263
def density_2d(self, x, y, rho0, Rs, center_x=0, center_y=0):
        
        x_ = x - center_x
        y_ = y - center_y
        r = np.sqrt(x_**2 + y_**2)
        X = r/Rs
        sigma0 = self.rho2sigma(rho0, Rs)
        if isinstance(X, int) or isinstance(X, float):
            if X == 1:
                X = 1.000001
        else:
            X[X == 1] = 1.000001
        sigma = sigma0 / (X**2-1)**2 * (-3 + (2+X**2)*self._F(X))
        return sigma
913529
def make_gaussian_kernel(sigma, npix=501, cdelt=0.01, xpix=None, ypix=None):
    

    sigma /= cdelt

    def fn(t, s): return 1. / (2 * np.pi * s ** 2) * np.exp(
        -t ** 2 / (s ** 2 * 2.0))
    dxy = make_pixel_distance(npix, xpix, ypix)
    k = fn(dxy, sigma)
    k /= (np.sum(k) * np.radians(cdelt) ** 2)

    return k
1103057
def mass_enclosed_3d(self, r, kwargs_profile):
        
        kwargs = copy.deepcopy(kwargs_profile)
        try:
            del kwargs['center_x']
            del kwargs['center_y']
        except:
            pass
        # integral of self._profile.density(x)* 4*np.pi * x^2 *dx, 0,r
        out = integrate.quad(lambda x: self._profile.density(x, **kwargs)*4*np.pi*x**2, 0, r)
        return out[0]
1050491
def _jmomentdensity(self,R,z,n,m,o,nsigma=None,mc=True,nmc=10000,
                       _returnmc=False,_vrs=None,_vts=None,_vzs=None,
                       **kwargs):
        """"""
        if nsigma == None:
            nsigma= _NSIGMA
        sigmaR1= self._sr*numpy.exp((self._refr-R)/self._hsr)
        sigmaz1= self._sz*numpy.exp((self._refr-R)/self._hsz)
        thisvc= potential.vcirc(self._pot,R,use_physical=False)
        #Use the asymmetric drift equation to estimate va
        gamma= numpy.sqrt(0.5)
        va= sigmaR1**2./2./thisvc\
            *(gamma**2.-1. #Assume close to flat rotation curve, sigphi2/sigR2 =~ 0.5
               +R*(1./self._hr+2./self._hsr))
        if math.fabs(va) > sigmaR1: va = 0.#To avoid craziness near the center
        if mc:
            mvT= (thisvc-va)/gamma/sigmaR1
            if _vrs is None:
                vrs= numpy.random.normal(size=nmc)
            else:
                vrs= _vrs
            if _vts is None:
                vts= numpy.random.normal(size=nmc)+mvT
            else:
                vts= _vts
            if _vzs is None:
                vzs= numpy.random.normal(size=nmc)
            else:
                vzs= _vzs
            Is= _jmomentsurfaceMCIntegrand(vzs,vrs,vts,numpy.ones(nmc)*R,numpy.ones(nmc)*z,self,sigmaR1,gamma,sigmaz1,mvT,n,m,o)
            if _returnmc:
                return (numpy.mean(Is)*sigmaR1**2.*gamma*sigmaz1,
                        vrs,vts,vzs)
            else:
                return numpy.mean(Is)*sigmaR1**2.*gamma*sigmaz1
        else: #pragma: no cover because this is too slow; a warning is shown
            warnings.warn("Calculations using direct numerical integration using tplquad is not recommended and extremely slow; it has also not been carefully tested",galpyWarning)
            return integrate.tplquad(_jmomentsurfaceIntegrand,
                                     1./gamma*(thisvc-va)/sigmaR1-nsigma,
                                     1./gamma*(thisvc-va)/sigmaR1+nsigma,
                                     lambda x: 0., lambda x: nsigma,
                                     lambda x,y: 0., lambda x,y: nsigma,
                                     (R,z,self,sigmaR1,gamma,sigmaz1,n,m,o),
                                     **kwargs)[0]*sigmaR1**2.*gamma*sigmaz1
934869
def hessian(self, x, y, Rs, theta_Rs, r_core, center_x=0, center_y=0):

        

        if Rs < 0.0001:
            Rs = 0.0001
        x_ = x - center_x
        y_ = y - center_y
        R = np.sqrt(x_ ** 2 + y_ ** 2)

        rho0 = self._alpha2rho0(theta_Rs=theta_Rs, Rs=Rs, r_core=r_core)

        kappa = self.density_2d(x_, y_, Rs, rho0, r_core)

        gamma1, gamma2 = self.cBurkGamma(R, Rs, rho0, r_core, x_, y_)
        f_xx = kappa + gamma1
        f_yy = kappa - gamma1
        f_xy = gamma2
        return f_xx, f_yy, f_xy