computes the projected density along the line-of-sight
        :param r: radius (arcsec)
        :param kwargs_profile: keyword argument list with lens model parameters
        :return: 2d projected density at projected radius r
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]
67728
def density(self, r, rho0, Rs):
        
        rho = rho0 / (r/Rs * (1 + (r/Rs))**3)
        return rho
360074
def _fit_radec(self):
        

        # call fitradec with mpcfile, abgfile, resfile
        self.orbfit.fitradec.restype = ctypes.c_int
        self.orbfit.fitradec.argtypes = [ ctypes.c_char_p, ctypes.c_char_p, ctypes.c_char_p ]

        mpc_file = tempfile.NamedTemporaryFile(suffix='.mpc')
        for observation in self.observations:
            mpc_file.write("{}\n".format(str(observation)))
        mpc_file.seek(0)

        abg_file = tempfile.NamedTemporaryFile()
        res_file = tempfile.NamedTemporaryFile()

        self.orbfit.fitradec(ctypes.c_char_p(mpc_file.name),
                             ctypes.c_char_p(abg_file.name),
                             ctypes.c_char_p(res_file.name))

        self.abg = abg_file
        self.abg.seek(0)
        self.residuals = res_file
        self.residuals.seek(0)
1293319
def K(self, r, R, kwargs):
        
        u = r / R
        if np.min(u) < 1:
            raise ValueError("3d radius is smaller than projected radius! Does not make sense.")

        if self._type == 'const_wrong':
            beta = kwargs['beta']
            k = 1./2. * u**(2*beta - 1.) * ((3./2 - beta) * np.sqrt(np.pi) * special.gamma(beta - 1./2)/special.gamma(beta)
                        + beta * self._B(x=1./u**2, a=beta+1./2, b=1./2) - self._B(x=1./u**2, a=beta-1./2, b=1./2))
        elif self._type == 'const':
            beta = kwargs['beta']
            k = np.sqrt(1 - 1./u**2) / (1. - 2*beta) + np.sqrt(np.pi)/2 * special.gamma(beta - 1./2)/special.gamma(beta)\
                * (3./2 - beta) * u**(2*beta - 1.) * (1 - special.betainc(beta+1./2, 1./2, 1./u**2))
        elif self._type == 'isotropic':
            k = np.sqrt(1 - 1./u**2)
        elif self._type == 'radial':
            k = np.pi/4 * u - 1./2 * np.sqrt(1 - 1./u**2) - u/2. * np.arcsin(1./u)
        elif self._type == 'OsipkovMerritt':
            r_ani = kwargs['r_ani']
            ua = r_ani / R
            k = (ua**2 + 1./2) / (ua**2 + 1)**(3./2) * (u**2 + ua**2) / u * np.arctan(np.sqrt((u**2 - 1) / (ua**2 + 1))) \
                - 1./2/(ua**2 + 1) * np.sqrt(1 -1./u**2)
        elif self._type == 'Colin':
            r_ani = kwargs['r_ani']
            ua = r_ani / R
            if ua == 1:
                k = (1 + 1./u) * np.arccosh(u) - 1./6 * (8./u + 7) * np.sqrt((u-1.)/(u+1.))
            else:
                k = 0.5 / (ua**2 - 1) * np.sqrt(1 - 1./u**2) + (1. + ua/u) * np.cosh(u) - np.sign(ua - 1) * ua * \
                            (ua**2 - 0.5) / np.abs(ua**2-1)**(3./2) * (1. + ua/u) * np.arccosh((ua*u + 1)/(u + ua))
        else:
            raise ValueError('anisotropy type %s not supported!' % self._type)
        return k
494598
def nfwParam_physical(self, M, c):
        
        r200 = self.nfw_param.r200_M(M * self.h) / self.h * self.a_z(self.z_lens)  # physical radius r200
        rho0 = self.nfw_param.rho0_c(c) * self.h**2 / self.a_z(self.z_lens)**3 # physical density in M_sun/Mpc**3
        Rs = r200/c
        return rho0, Rs, r200
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
535278
def observe(self, dataset, kind, components=None, distance=1.0, l3=0.0, **kwargs):
        

        meshes = self.meshes
        if kind=='lp':
            def sv(p, p0, w):
                # Subsidiary variable:
                return (p0-p)/(w/2)

            def lorentzian(sv):
                return 1-1./(1+sv**2)

            def gaussian(sv):
                return 1-np.exp(-np.log(2)*sv**2)

            profile_func = kwargs.get('profile_func')
            profile_rest = kwargs.get('profile_rest')
            profile_sv = kwargs.get('profile_sv')
            wavelengths = kwargs.get('wavelengths')
            if profile_func == 'gaussian':
                func = gaussian
            elif profile_func == 'lorentzian':
                func = lorentzian
            else:
                raise NotImplementedError("profile_func='{}' not supported".format(profile_func))

            visibilities = meshes.get_column_flat('visibilities', components)

            abs_intensities = meshes.get_column_flat('abs_intensities:{}'.format(dataset), components)
            # mus here will be from the tnormals of the triangle and will not
            # be weighted by the visibility of the triangle
            mus = meshes.get_column_flat('mus', components)
            areas = meshes.get_column_flat('areas_si', components)
            ldint = meshes.get_column_flat('ldint:{}'.format(dataset), components)

            rvs = (meshes.get_column_flat("rvs:{}".format(dataset), components)*u.solRad/u.d).to(u.m/u.s).value
            dls = rvs*profile_rest/c.c.si.value

            line = func(sv(wavelengths, profile_rest, profile_sv))
            lines = np.array([np.interp(wavelengths, wavelengths+dl, line) for dl in dls])
            if not np.any(visibilities):
                avg_line = np.full_like(wavelengths, np.nan)
            else:
                avg_line = np.average(lines, axis=0, weights=abs_intensities*areas*mus*ldint*visibilities)

            return {'flux_densities': avg_line}


        elif kind=='rv':
            visibilities = meshes.get_column_flat('visibilities', components)

            if np.all(visibilities==0):
                # then no triangles are visible, so we should return nan
                return {'rv': np.nan}

            rvs = meshes.get_column_flat("rvs:{}".format(dataset), components)
            abs_intensities = meshes.get_column_flat('abs_intensities:{}'.format(dataset), components)
            # mus here will be from the tnormals of the triangle and will not
            # be weighted by the visibility of the triangle
            mus = meshes.get_column_flat('mus', components)
            areas = meshes.get_column_flat('areas_si', components)
            ldint = meshes.get_column_flat('ldint:{}'.format(dataset), components)
            # NOTE: don't need ptfarea because its a float (same for all
            # elements, regardless of component)

            # NOTE: the intensities are already projected but are per unit area
            # so we need to multiply by the /projected/ area of each triangle (thus the extra mu)
            return {'rv': np.average(rvs, weights=abs_intensities*areas*mus*ldint*visibilities)}

        elif kind=='lc':
            visibilities = meshes.get_column_flat('visibilities')

            if np.all(visibilities==0):
                # then no triangles are visible, so we should return nan -
                # probably shouldn't ever happen for lcs
                return {'flux': np.nan}

            intensities = meshes.get_column_flat("intensities:{}".format(dataset), components)
            mus = meshes.get_column_flat('mus', components)
            areas = meshes.get_column_flat('areas_si', components)
            ldint = meshes.get_column_flat('ldint:{}'.format(dataset), components)

            # assume that all bodies are using the same passband and therefore
            # will have the same ptfarea.  If this assumption is ever a problem -
            # then we will need to build a flat column based on the component
            # of each element so that ptfarea is an array with the same shape
            # as those above
            if isinstance(self.bodies[0], Envelope):
                # for envelopes, we'll make the same assumption and just grab
                # that value stored in the first "half"
                ptfarea = self.bodies[0]._halves[0].get_ptfarea(dataset)
            else:
                ptfarea = self.bodies[0].get_ptfarea(dataset)

            # intens_proj is the intensity in the direction of the observer per unit surface area of the triangle
            # areas is the area of each triangle
            # areas*mus is the area of each triangle projected in the direction of the observer
            # visibilities is 0 for hidden, 0.5 for partial, 1.0 for visible
            # areas*mus*visibilities is the visibile projected area of each triangle (ie half the area for a partially-visible triangle)
            # so, intens_proj*areas*mus*visibilities is the intensity in the direction of the observer per the observed projected area of that triangle
            # and the sum of these values is the observed flux

            # note that the intensities are already projected but are per unit area
            # so we need to multiply by the /projected/ area of each triangle (thus the extra mu)

            return {'flux': np.sum(intensities*areas*mus*visibilities)*ptfarea/(distance**2)+l3}

        else:
            raise NotImplementedError("observe for dataset with kind '{}' not implemented".format(kind))
457888
def multi_gaussian_lens(self, kwargs_lens, model_bool_list=None, e1=0, e2=0, n_comp=20):
        
        if 'center_x' in kwargs_lens[0]:
            center_x = kwargs_lens[0]['center_x']
            center_y = kwargs_lens[0]['center_y']
        else:
            raise ValueError('no keyword center_x defined!')
        theta_E = self._lensModelExtensions.effective_einstein_radius(kwargs_lens)
        r_array = np.logspace(-4, 2, 200) * theta_E
        x_coords, y_coords = param_util.transform_e1e2(r_array, np.zeros_like(r_array), e1=-e1, e2=-e2)
        x_coords += center_x
        y_coords += center_y
        #r_array = np.logspace(-2, 1, 50) * theta_E
        if model_bool_list is None:
            model_bool_list = [True] * len(kwargs_lens)
        kappa_s = np.zeros_like(r_array)
        for i in range(len(kwargs_lens)):
            if model_bool_list[i] is True:
                kappa_s += self.LensModel.kappa(x_coords, y_coords, kwargs_lens, k=i)
        amplitudes, sigmas, norm = mge.mge_1d(r_array, kappa_s, N=n_comp)
        return amplitudes, sigmas, center_x, center_y
1078747
def pseudo_core_density(self):
        """"""
        mesh, values, attrib = self._parse_radfunc("pseudo_core_density")
        return RadialFunction(mesh, values)
427828
def re_normalize_flux(self, kwargs_ps, norm_factor):
        
        for i, model in enumerate(self.point_source_type_list):
            if model == 'UNLENSED':
                kwargs_ps[i]['point_amp'] *= norm_factor
            elif model in ['LENSED_POSITION', 'SOURCE_POSITION']:
                if self._fixed_magnification_list[i] is True:
                    kwargs_ps[i]['source_amp'] *= norm_factor
                else:
                    kwargs_ps[i]['point_amp'] *= norm_factor
        return kwargs_ps
298547
def check_solver(self, kwargs_lens, kwargs_ps, kwargs_cosmo={}):
        
        if self._solver is True:
            image_x, image_y = kwargs_ps[0]['ra_image'], kwargs_ps[0]['dec_image']
            image_x, image_y = self.real_image_positions(image_x, image_y, kwargs_cosmo)
            dist = self._solver_module.check_solver(image_x, image_y, kwargs_lens)
            return np.max(dist)
        else:
            return 0
1222534
def vel_disp(self, kwargs_profile, kwargs_aperture, kwargs_light, kwargs_anisotropy, num=1000):
        
        sigma_s2_sum = 0
        for i in range(0, num):
            sigma_s2_draw = self._vel_disp_one(kwargs_profile, kwargs_aperture, kwargs_light, kwargs_anisotropy)
            sigma_s2_sum += sigma_s2_draw
        sigma_s2_average = sigma_s2_sum/num
        return np.sqrt(sigma_s2_average)
1162948
def arbanglemap(self):
        '''
        
        '''
        # raySpacingDeg = self.arbfov / self.nCutPix
        maxAng = self.boresightEl + self.arbfov / 2
        minAng = self.boresightEl - self.arbfov / 2
        angles = np.linspace(maxAng, minAng, num=self.ncutpix, endpoint=True)
        assert np.isclose(angles[0], maxAng) & np.isclose(angles[-1], minAng)
        return angles
743844
def density(self):
        
        r = float(self.sp_index.npoints) / float(self.sp_index.length)
        return r
700889
def ae_core_density(self):
        """"""
        mesh, values, attrib = self._parse_radfunc("ae_core_density")
        return RadialFunction(mesh, values)
311668
def _call_marginalizevperp(self,o,**kwargs):
        """"""
        #Get l, vlos
        l= o.ll(obs=[1.,0.,0.],ro=1.)*_DEGTORAD
        vlos= o.vlos(ro=1.,vo=1.,obs=[1.,0.,0.,0.,0.,0.])
        R= o.R(use_physical=False)
        phi= o.phi(use_physical=False)
        #Get local circular velocity, projected onto the los
        vcirc= R**self._beta
        vcirclos= vcirc*math.sin(phi+l)
        #Marginalize
        alphalos= phi+l
        if not 'nsigma' in kwargs or ('nsigma' in kwargs and \
                                          kwargs['nsigma'] is None):
            nsigma= _NSIGMA
        else:
            nsigma= kwargs['nsigma']
        kwargs.pop('nsigma',None)
        sigmaR2= self.targetSigma2(R,use_physical=False)
        sigmaR1= sc.sqrt(sigmaR2)
        #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 math.fabs(math.sin(alphalos)) < math.sqrt(1./2.):
            cosalphalos= math.cos(alphalos)
            tanalphalos= math.tan(alphalos)            
            return integrate.quad(_marginalizeVperpIntegrandSinAlphaSmall,
                                  -self._gamma*va/sigmaR1-nsigma,
                                  -self._gamma*va/sigmaR1+nsigma,
                                  args=(self,R,cosalphalos,tanalphalos,
                                        vlos-vcirclos,vcirc,
                                        sigmaR1/self._gamma),
                                  **kwargs)[0]/math.fabs(cosalphalos)\
                                  *sigmaR1/self._gamma
        else:
            sinalphalos= math.sin(alphalos)
            cotalphalos= 1./math.tan(alphalos)
            return integrate.quad(_marginalizeVperpIntegrandSinAlphaLarge,
                                  -nsigma,nsigma,
                                  args=(self,R,sinalphalos,cotalphalos,
                                        vlos-vcirclos,vcirc,sigmaR1),
                                  **kwargs)[0]/math.fabs(sinalphalos)*sigmaR1
1277537
def projected_area_factor(ax, original_crs):
    
    if not hasattr(ax, 'projection'):
        return 1
    if isinstance(ax.projection, ccrs.PlateCarree):
        return 1
    x1, x2, y1, y2 = ax.get_extent()
    pbounds = \
        get_projection_from_crs(original_crs).transform_points(ax.projection,
                    np.array([x1, x2]), np.array([y1, y2]))

    return np.sqrt(abs((x2 - x1) * (y2 - y1))
                   /abs((pbounds[0] - pbounds[1])[:2].prod()))