computes the mass enclosed the projected line-of-sight
        :param r: radius (arcsec)
        :param kwargs_profile: keyword argument list with lens model parameters
        :return: projected mass enclosed radius r
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]
672117
def mass_2d(self, r, rho0, Ra, Rs):
        
        Ra, Rs = self._sort_ra_rs(Ra, Rs)
        sigma0 = self.rho2sigma(rho0, Ra, Rs)
        m_2d = 2 * np.pi * sigma0 * Ra * Rs / (Rs - Ra) * (np.sqrt(Ra ** 2 + r ** 2) - Ra - np.sqrt(Rs ** 2 + r ** 2) + Rs)
        return m_2d
1301510
def _extract_array(self, kwargs_list):
        
        if self._solver_type == 'PROFILE_SHEAR':
            e1 = kwargs_list[1]['e1']
            e2 = kwargs_list[1]['e2']
            phi_ext, gamma_ext = param_util.ellipticity2phi_gamma(e1, e2)
        else:
            phi_ext = 0
        lens_model = self._lens_mode_list[0]
        if lens_model in ['SPEP', 'SPEMD', 'SIE', 'NIE']:
            e1 = kwargs_list[0]['e1']
            e2 = kwargs_list[0]['e2']
            center_x = kwargs_list[0]['center_x']
            center_y = kwargs_list[0]['center_y']
            theta_E = kwargs_list[0]['theta_E']
            x = [theta_E, e1, e2, center_x, center_y, phi_ext]
        elif lens_model in ['NFW_ELLIPSE']:
            e1 = kwargs_list[0]['e1']
            e2 = kwargs_list[0]['e2']
            center_x = kwargs_list[0]['center_x']
            center_y = kwargs_list[0]['center_y']
            theta_Rs = kwargs_list[0]['theta_Rs']
            x = [theta_Rs, e1, e2, center_x, center_y, phi_ext]
        elif lens_model in ['SHAPELETS_CART']:
            coeffs = list(kwargs_list[0]['coeffs'])
            [c10, c01, c20, c11, c02] = coeffs[1: 6]
            x = [c10, c01, c20, c11, c02, phi_ext]
        else:
            raise ValueError("Lens model %s not supported for 4-point solver!" % lens_model)
        return x
1002568
def detectability(self,**kwargs):
        
        distance_modulus = kwargs.get('distance_modulus')
        distance = mod2dist(distance_modulus)
        stellar_mass = kwargs.get('stellar_mass')
        extension = kwargs.get('extension')

        # Normalized to 10^3 Msolar at mod=18
        norm = 10**3/mod2dist(18)**2
        detect = stellar_mass / distance**2
        detect /= norm
372599
def mass_3d_lens(self, r, theta_E, gamma, e1, e2):
        
        return self.spp.mass_3d_lens(r, theta_E, gamma)
1146634
def _get_distance_term(self, C, rrup, backarc):
        
        # Geometric attenuation function
        distance_scale = -np.log10(np.sqrt(rrup ** 2 + 3600.0))
        # Anelastic attenuation in the backarc
        distance_scale[backarc] += (C["c2"] * rrup[backarc])
        # Anelastic Attenuation in the forearc
        idx = np.logical_not(backarc)
        distance_scale[idx] += (C["c1"] * rrup[idx])
        return distance_scale
1303157
def update_lens_scaling(self, kwargs_cosmo, kwargs_lens, inverse=False):
        
        kwargs_lens_updated = copy.deepcopy(kwargs_lens)
        if self._mass_scaling is False:
            return kwargs_lens_updated
        scale_factor_list = np.array(kwargs_cosmo['scale_factor'])
        if inverse is True:
            scale_factor_list = 1. / np.array(kwargs_cosmo['scale_factor'])
        for i, kwargs in enumerate(kwargs_lens_updated):
            if self._mass_scaling_list[i] is not False:
                scale_factor = scale_factor_list[self._mass_scaling_list[i] - 1]
                if 'theta_E' in kwargs:
                    kwargs['theta_E'] *= scale_factor
                elif 'theta_Rs' in kwargs:
                    kwargs['theta_Rs'] *= scale_factor
                elif 'sigma0' in kwargs:
                    kwargs['sigma0'] *= scale_factor
                elif 'k_eff' in kwargs:
                    kwargs['k_eff'] *= scale_factor
        return kwargs_lens_updated
92819
def mass_2d_lens(self, r, sigma0, Rs):
        
        X = r/Rs
        alpha_r = 2*sigma0 * Rs * X * (1-self._F(X)) / (X**2-1)
        mass_2d = alpha_r * r * np.pi
        return mass_2d
1027489
def point_source_list(self, kwargs_ps, kwargs_lens, k=None):
        
        ra_list, dec_list = self.image_position(kwargs_ps, kwargs_lens, k=k)
        amp_list = self.image_amplitude(kwargs_ps, kwargs_lens)
        ra_array, dec_array, amp_array = [], [], []
        for i, ra in enumerate(ra_list):
            for j in range(len(ra)):
                ra_array.append(ra_list[i][j])
                dec_array.append(dec_list[i][j])
                amp_array.append(amp_list[i][j])
        return ra_array, dec_array, amp_array
428590
def kwargs2args(self, kwargs_lens=None, kwargs_source=None, kwargs_lens_light=None, kwargs_ps=None, kwargs_cosmo=None):
        

        args = self.lensParams.setParams(kwargs_lens)
        args += self.souceParams.setParams(kwargs_source)
        args += self.lensLightParams.setParams(kwargs_lens_light)
        args += self.pointSourceParams.setParams(kwargs_ps)
        args += self.cosmoParams.setParams(kwargs_cosmo)
        return args
646205
def _vmomentsurfacemass(self,R,n,m,romberg=False,nsigma=None,
                           relative=False,phi=0.,deriv=None):
        """"""
        #odd moments of vR are zero
        if isinstance(n,int) and n%2 == 1:
            return 0.
        if nsigma == None:
            nsigma= _NSIGMA
        logSigmaR= self.targetSurfacemass(R,log=True,use_physical=False)
        sigmaR2= self.targetSigma2(R,use_physical=False)
        sigmaR1= sc.sqrt(sigmaR2)
        logsigmaR2= sc.log(sigmaR2)
        if relative:
            norm= 1.
        else:
            norm= sc.exp(logSigmaR+logsigmaR2*(n+m)/2.)/self._gamma**m
        #Use the asymmetric drift equation to estimate va
        va= sigmaR2/2./R**self._beta*(1./self._gamma**2.-1.
                                      -R*self._surfaceSigmaProfile.surfacemassDerivative(R,log=True)
                                      -R*self._surfaceSigmaProfile.sigma2Derivative(R,log=True))
        if math.fabs(va) > sigmaR1: va = 0. #To avoid craziness near the center
        if deriv is None:
            if romberg:
                return sc.real(bovy_dblquad(_vmomentsurfaceIntegrand,
                                            self._gamma*(R**self._beta-va)/sigmaR1-nsigma,
                                            self._gamma*(R**self._beta-va)/sigmaR1+nsigma,
                                            lambda x: -nsigma, lambda x: nsigma,
                                            [R,self,logSigmaR,logsigmaR2,sigmaR1,
                                             self._gamma,n,m],
                                            tol=10.**-8)/sc.pi*norm/2.)
            else:
                return integrate.dblquad(_vmomentsurfaceIntegrand,
                                         self._gamma*(R**self._beta-va)/sigmaR1-nsigma,
                                         self._gamma*(R**self._beta-va)/sigmaR1+nsigma,
                                         lambda x: -nsigma, lambda x: nsigma,
                                         (R,self,logSigmaR,logsigmaR2,sigmaR1,
                                          self._gamma,n,m),
                                         epsrel=_EPSREL)[0]/sc.pi*norm/2.
        else:
            if romberg:
                return sc.real(bovy_dblquad(_vmomentderivsurfaceIntegrand,
                                            self._gamma*(R**self._beta-va)/sigmaR1-nsigma,
                                            self._gamma*(R**self._beta-va)/sigmaR1+nsigma,
                                            lambda x: -nsigma, lambda x: nsigma,
                                            [R,self,logSigmaR,logsigmaR2,sigmaR1,
                                             self._gamma,n,m,deriv],
                                            tol=10.**-8)/sc.pi*norm/2.)
            else:
                return integrate.dblquad(_vmomentderivsurfaceIntegrand,
                                         self._gamma*(R**self._beta-va)/sigmaR1-nsigma,
                                         self._gamma*(R**self._beta-va)/sigmaR1+nsigma,
                                         lambda x: -nsigma, lambda x: nsigma,
                                         (R,self,logSigmaR,logsigmaR2,sigmaR1,
                                          self._gamma,n,m,deriv),
                                         epsrel=_EPSREL)[0]/sc.pi*norm/2.
319053
def external_lensing_effect(self, kwargs_lens, lens_model_internal_bool=None):
        
        alpha0_x, alpha0_y = 0, 0
        kappa_ext = 0
        shear1, shear2 = 0, 0
        if lens_model_internal_bool is None:
            lens_model_internal_bool = [True] * len(kwargs_lens)
        for i, kwargs in enumerate(kwargs_lens):
            if not lens_model_internal_bool[i] is True:
                f_x, f_y = self._lensModel.alpha(0, 0, kwargs_lens, k=i)
                f_xx, f_xy, f_yx, f_yy = self._lensModel.hessian(0, 0, kwargs_lens, k=i)
                alpha0_x += f_x
                alpha0_y += f_y
                kappa_ext += (f_xx + f_yy)/2.
                shear1 += 1./2 * (f_xx - f_yy)
                shear2 += f_xy
        return alpha0_x, alpha0_y, kappa_ext, shear1, shear2
828577
def logL(self, args):
        
        #extract parameters
        kwargs_lens, kwargs_source, kwargs_lens_light, kwargs_ps, kwargs_cosmo = self.param.args2kwargs(args)
        #generate image and computes likelihood
        self._reset_point_source_cache(bool=True)
        logL = 0
        if self._check_bounds is True:
            penalty, bound_hit = self.check_bounds(args, self._lower_limit, self._upper_limit)
            logL -= penalty
            if bound_hit:
                return logL, None
        if self._image_likelihood is True:
            logL += self.image_likelihood.logL(kwargs_lens, kwargs_source, kwargs_lens_light, kwargs_ps)
        if self._time_delay_likelihood is True:
            logL += self.time_delay_likelihood.logL(kwargs_lens, kwargs_ps, kwargs_cosmo)
        if self._check_positive_flux is True:
            bool = self.param.check_positive_flux(kwargs_source, kwargs_lens_light, kwargs_ps)
            if bool is False:
                logL -= 10**10
        if self._flux_ratio_likelihood is True:
            logL += self.flux_ratio_likelihood.logL(kwargs_lens, kwargs_ps, kwargs_cosmo)
        logL += self._position_likelihood.logL(kwargs_lens, kwargs_ps, kwargs_cosmo)
        logL += self._prior_likelihood.logL(kwargs_lens, kwargs_source, kwargs_lens_light, kwargs_ps, kwargs_cosmo)
        self._reset_point_source_cache(bool=False)
        return logL, None
1312174
def effectiv_num_data_points(self, kwargs_lens, kwargs_source, kwargs_lens_light, kwargs_ps):
        
        num_linear = 0
        if self._image_likelihood is True:
            num_linear = self.image_likelihood.num_param_linear(kwargs_lens, kwargs_source, kwargs_lens_light, kwargs_ps)
        num_param, _ = self.param.num_param()
        return self.num_data - num_param - num_linear
615621
def cars_profile(filename, doy, latitude, longitude, depth):
    
    assert np.size(doy) == 1
    assert np.size(latitude) == 1
    assert np.size(longitude) == 1
    #assert np.size(depth) == 1

    assert (longitude >= 0) & (longitude <= 360)
    assert depth >= 0

    nc = netCDF4.Dataset(filename)

    t = 2 * np.pi * doy/366

    # Improve this. Optimize to get only necessary Z
    Z = slice(0, nc['depth'].size)
    I = np.absolute(nc['lat'][:] - latitude).argmin()
    J = np.absolute(nc['lon'][:] - longitude).argmin()

    # Not efficient, but it works
    assert (nc['depth'][:64] == nc['depth_ann'][:]).all()
    assert (nc['depth'][:55] == nc['depth_semiann'][:]).all()
    value = nc['mean'][:, I, J]
    value[:64] += nc['an_cos'][Z, I, J] * np.cos(t) + \
            nc['an_sin'][:, I, J] * np.sin(t)
    value[:55] += nc['sa_cos'][Z, I, J] * np.cos(2*t) + \
            nc['sa_sin'][:, I, J] * np.sin(2*t)
    value = value

    output = {'depth': np.asanyarray(depth)}
    from scipy.interpolate import griddata
    output['value'] = griddata(nc['depth'][Z], value[Z], depth)

    for v in ['std_dev']:
        output[v] = griddata(nc['depth'][Z], nc[v][Z, I, J], depth)

    return output
574642
def beta_r(self, r, kwargs):
        
        if self._type == 'const':
            return self.const_beta(kwargs)
        elif self._type == 'r_ani':
            return self.beta_r_ani(r, kwargs)
        else:
            raise ValueError('anisotropy type %s not supported!' % self._type)
132513
def mass_3d(self, r, rho0, gamma):
        
        mass_3d = 4 * np.pi * rho0 /(-gamma + 3) * r ** (-gamma + 3)
        return mass_3d