converts projected density parameter (in units of deflection) into 3d density parameter
        :param theta_E:
        :param gamma:
        :return:
482825
def theta2rho(self, theta_E, gamma):
        
        fac1 = np.sqrt(np.pi) * special.gamma(1. / 2 * (-1 + gamma)) / special.gamma(gamma / 2.) * 2 / (3 - gamma)
        fac2 = theta_E**(gamma - 1)
        rho0 = fac2 / fac1
        return rho0
717998
def spherical_to_cartesian(r,theta,phi):
    
    x = r * np.sin(phi) * np.cos(theta)
    y = r * np.sin(phi) * np.sin(theta)
    z = r * np.cos(phi)
    return (x,y,z)
762531
def gamma(phi1,phi2,theta1,theta2):
    
    
    if phi1 == phi2 and theta1 == theta2:
        gamma = 0
    else:
        gamma = atan( sin(theta2)*sin(phi2-phi1) / \
                      (cos(theta1)*sin(theta2)*cos(phi1-phi2) - \
                       sin(theta1)*cos(theta2)) )

    dummy_arg = (cos(gamma)*cos(theta1)*sin(theta2)*cos(phi1-phi2) + \
                 sin(gamma)*sin(theta2)*sin(phi2-phi1) - \
                 cos(gamma)*sin(theta1)*cos(theta2))
        
    if dummy_arg >= 0:
        return gamma
    else:
        return pi + gamma
727606
def light_3d(self, r, kwargs_list, k=None):
        
        r = np.array(r, dtype=float)
        flux = np.zeros_like(r)
        for i, func in enumerate(self.func_list):
            if k is None or k == i:
                kwargs = {k: v for k, v in kwargs_list[i].items() if not k in ['center_x', 'center_y']}
                if self.profile_type_list[i] in ['HERNQUIST', 'HERNQUIST_ELLIPSE', 'PJAFFE', 'PJAFFE_ELLIPSE',
                                                     'GAUSSIAN', 'GAUSSIAN_ELLIPSE', 'MULTI_GAUSSIAN',
                                                     'MULTI_GAUSSIAN_ELLIPSE', 'POWER_LAW']:
                    flux += func.light_3d(r, **kwargs)
                else:
                    raise ValueError('Light model %s does not support a 3d light distribution!'
                                         % self.profile_type_list[i])
        return flux
1076073
def dists(node):
    
    hd = tuple((node['probability'], node['depth'])
               for node in node.hypoDepthDist)
    npd = tuple(
        ((node['probability'], node['rake'], node['strike'], node['dip']))
        for node in node.nodalPlaneDist)
    return hd, npd, ~node.magScaleRel
1141866
def gammalnStirling(z):
    
    return (0.5 * (np.log(2. * np.pi) - np.log(z))) \
           + (z * (np.log(z + (1. / ((12. * z) - (1. / (10. * z))))) - 1.))
1046587
def _theta_E_convert(self, theta_E, w_c, w_t):
        
        if not w_t > w_c:
            w_t, w_c = w_c, w_t
        s_scale_1 = w_c
        s_scale_2 = w_t
        f_x_1, f_y_1 = self.nie.derivatives(1, 0, theta_E=1, e1=0, e2=0, s_scale=s_scale_1)
        f_x_2, f_y_2 = self.nie.derivatives(1, 0, theta_E=1, e1=0, e2=0, s_scale=s_scale_2)
        f_x = f_x_1 - f_x_2
        theta_E_convert = theta_E / f_x
        return theta_E_convert, w_c, w_t
1094820
def lnprob(self,theta):
        
        global niter
        params,priors,loglike = self.params,self.priors,self.loglike
        # Avoid extra likelihood calls with bad priors
        _lnprior = self.lnprior(theta)
        if np.isfinite(_lnprior):
            _lnlike = self.lnlike(theta)
        else:
            _lnprior = -np.inf
            _lnlike = -np.inf

        _lnprob = _lnprior + _lnlike
     
        if (niter%100==0):
            msg = "%i function calls ...\n"%niter
            msg+= ', '.join('%s: %.3f'%(k,v) for k,v in zip(params,theta))
            msg+= '\nlog(like): %.3f, log(prior): %.3f'%(_lnprior,_lnlike)
            logger.debug(msg)
        niter+=1
        return _lnprob
731456
def to_cartesian(r, theta, theta_units="radians"):
    
    assert theta_units in ['radians', 'degrees'],\
        "kwarg theta_units must specified in radians or degrees"

    # Convert to radians
    if theta_units == "degrees":
        theta = to_radians(theta)

    theta = to_proper_radians(theta)
    x = r * cos(theta)
    y = r * sin(theta)

    return x, y
866214
def mass_in_theta_E(self, theta_E):
        
        mass = self.arcsec2phys_lens(theta_E) ** 2 * np.pi * self.epsilon_crit
        return mass
569741
def sample_lonlat(self, n):
        
        # From http://en.wikipedia.org/wiki/Ellipse#General_parametric_form
        # However, Martin et al. (2009) use PA theta "from North to East"
        # Definition of phi (position angle) is offset by pi/4
        # Definition of t (eccentric anamoly) remains the same (x,y-frame usual)
        # In the end, everything is trouble because we use glon, glat...

        radius = self.sample_radius(n)
        a = radius; b = self.jacobian * radius

        t = 2. * np.pi * np.random.rand(n)
        cost,sint = np.cos(t),np.sin(t)
        phi = np.pi/2. - np.deg2rad(self.theta)
        cosphi,sinphi = np.cos(phi),np.sin(phi)
        x = a*cost*cosphi - b*sint*sinphi
        y = a*cost*sinphi + b*sint*cosphi
        
        if self.projector  is None:
            logger.debug("Creating AITOFF projector for sampling")
            projector = Projector(self.lon,self.lat,'ait')
        else:
            projector = self.projector
        lon, lat = projector.imageToSphere(x, y)
        return lon, lat
611127
def sis_theta_E2sigma_v(self, theta_E):
        
        v_sigma_c2 = theta_E * const.arcsec / (4*np.pi) * self.D_s / self.D_ds
        return np.sqrt(v_sigma_c2)*const.c / 1000
901307
def d_distributionparams(self):
        """"""
        if not self._d_distributionparams:
            dx = 1.0e-3
            def f_alpha(alpha):
                return DiscreteGamma(alpha, self.beta_lambda, self.ncats)
            def f_beta(beta):
                return DiscreteGamma(self.alpha_lambda, beta, self.ncats)
            assert set(self.distributionparams) == {'alpha_lambda', 'beta_lambda'}
            for (param, f) in [('alpha_lambda', f_alpha), ('beta_lambda', f_beta)]:
                pvalue = getattr(self, param)
                dparam = scipy.misc.derivative(f, pvalue, dx, n=1, order=5)
                assert dparam.shape == (self.ncats,)
                for stepchange in [0.5, 2]: # make sure robust to step size
                    dparam2 = scipy.misc.derivative(f, pvalue, stepchange * dx,
                            n=1, order=5)
                    assert scipy.allclose(dparam, dparam2, atol=1e-5, rtol=1e-4), (
                            "Numerical derivative of {0} at {1} "
                            "differs for step {2} and {3}: {4} and {5}"
                            ", respectively.").format(param, pvalue,
                            dx, dx * stepchange, dparam, dparam2)
                self._d_distributionparams[param] = dparam
        return self._d_distributionparams
1204850
def probability_density(self, X):
        
        self.check_fit()

        U, V = self.split_matrix(X)

        if self.theta == 0:
            return np.multiply(U, V)

        else:
            num = np.multiply(np.multiply(-self.theta, self._g(1)), 1 + self._g(np.add(U, V)))
            aux = np.multiply(self._g(U), self._g(V)) + self._g(1)
            den = np.power(aux, 2)
            return num / den
172657
def density(self, x, y, n_sersic, r_eff, k_eff, center_x=0, center_y=0):
        
        raise ValueError("not implemented! Use a Multi-Gaussian-component decomposition.")
809216
def mass_3d_lens(self, R, Rs, theta_Rs):
        
        rho0 = self._alpha2rho0(theta_Rs, Rs)
        m_3d = self.mass_3d(R, Rs, rho0)
        return m_3d
773777
def density_2d(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(np.sqrt(x^2+r^2))* dx, 0, infty
        out = integrate.quad(lambda x: 2*self._profile.density(np.sqrt(x**2+r**2), **kwargs), 0, 100)
        return out[0]