computes one realisation of the velocity dispersion realized in the slit
        :param gamma:
        :param rho0_r0_gamma:
        :param r_eff:
        :param r_ani:
        :param R_slit:
        :param dR_slit:
        :param FWHM:
        :return:
179903
def _vel_disp_one(self, kwargs_profile, kwargs_aperture, kwargs_light, kwargs_anisotropy):
        

        while True:
            r = self.lightProfile.draw_light(kwargs_light)  # draw r
            R, x, y = util.R_r(r)  # draw projected R
            x_, y_ = util.displace_PSF(x, y, self.FWHM)  # displace via PSF
            bool = self.aperture.aperture_select(x_, y_, kwargs_aperture)
            if bool is True:
                break
        sigma_s2 = self.sigma_s2(r, R, kwargs_profile, kwargs_anisotropy, kwargs_light)
        return sigma_s2
48284
def _phi2deriv(self,R,phi=0.,t=0.):
        
        return self._Pot.phi2deriv(R,0.,phi=phi,t=t,use_physical=False)
94219
def function(self, x, y, amp, R_sersic, Re, n_sersic, gamma, e1, e2, center_x=0, center_y=0, alpha=3.):
        
        phi_G, q = param_util.ellipticity2phi_q(e1, e2)
        Rb = R_sersic
        x_shift = x - center_x
        y_shift = y - center_y

        cos_phi = np.cos(phi_G)
        sin_phi = np.sin(phi_G)

        xt1 = cos_phi*x_shift+sin_phi*y_shift
        xt2 = -sin_phi*x_shift+cos_phi*y_shift
        xt2difq2 = xt2/(q*q)
        R_ = np.sqrt(xt1*xt1+xt2*xt2difq2)
        #R_ = R_.astype(np.float32)
        if isinstance(R_, int) or isinstance(R_, float):
            R_ = max(self._smoothing, R_)
        else:
            R_[R_ < self._smoothing] = self._smoothing
        if isinstance(R_, int) or isinstance(R_, float):
            R = max(self._smoothing, R_)
        else:
            R=np.empty_like(R_)
            _R = R_[R_ > self._smoothing]  #in the SIS regime
            R[R_ <= self._smoothing] = self._smoothing
            R[R_ > self._smoothing] = _R

        k, bn = self.k_bn(n_sersic, Re)
        result = amp * (1 + (Rb / R) ** alpha) ** (gamma / alpha) * np.exp(-bn * (((R ** alpha + Rb ** alpha) / Re ** alpha) ** (1. / (alpha * n_sersic)) - 1.))
        return np.nan_to_num(result)
121512
def rho_hv(scatterer):
    
    Z = scatterer.get_Z()
    a = (Z[2,2] + Z[3,3])**2 + (Z[3,2] - Z[2,3])**2
    b = (Z[0,0] - Z[0,1] - Z[1,0] + Z[1,1])
    c = (Z[0,0] + Z[0,1] + Z[1,0] + Z[1,1])
    return np.sqrt(a / (b*c))
707982
def F_t(X, Y, S, M_E, E, m0, rho):
    '''  '''
    r = X.shape[1]
    out1 = (((np.dot(np.dot(X, S), Y.T) - M_E) * E)**2).sum() / 2
    out2 = rho * G(Y, m0, r)
    out3 = rho * G(X, m0, r)

    return out1 + out2 + out3
948400
def sim_combi(self):
        
        n_samples = 500
        sigma_glob = 1.8

        X = np.zeros((n_samples,3))

        X[:,0] = np.random.uniform(-sigma_glob,sigma_glob,n_samples)
        X[:,1] = np.random.uniform(-sigma_glob,sigma_glob,n_samples)

        func = self.funcs['tanh']

        # XOR type
#         X[:,2] = (func(X[:,0])*sp.stats.norm.pdf(X[:,1],0,0.2)
#                   + func(X[:,1])*sp.stats.norm.pdf(X[:,0],0,0.2))
        # AND type / diagonal
#         X[:,2] = (func(X[:,0]+X[:,1])*sp.stats.norm.pdf(X[:,1]-X[:,0],0,0.2))
        # AND type / horizontal
        X[:,2] = (func(X[:,0])*sp.stats.norm.cdf(X[:,1],1,0.2))

        pl.scatter(X[:,0],X[:,1],c=X[:,2],edgecolor='face')
        pl.show()

        pl.plot(X[:,1],X[:,2],'.')
        pl.show()

        return X
480117
def _R2deriv(self,R,phi=0.,t=0.):
        
        return self._Pot.R2deriv(R,0.,phi=phi,t=t,use_physical=False)
38821
def power_law_anisotropy(self, r, kwargs_profile, kwargs_anisotropy, kwargs_light):
        
        # first term
        theta_E = kwargs_profile['theta_E']
        gamma = kwargs_profile['gamma']
        r_ani = kwargs_anisotropy['r_ani']
        a = 0.551 * kwargs_light['r_eff']
        rho0_r0_gamma = self._rho0_r0_gamma(theta_E, gamma)
        prefac1 = 4*np.pi * const.G * a**(-gamma) * rho0_r0_gamma / (3-gamma)
        prefac2 = r * (r + a)**3/(r**2 + r_ani**2)
        hyp1 = vel_util.hyp_2F1(a=2+gamma, b=gamma, c=3+gamma, z=1./(1+r/a))
        hyp2 = vel_util.hyp_2F1(a=3, b=gamma, c=1+gamma, z=-a/r)
        fac = r_ani**2/a**2 * hyp1 / ((2+gamma) * (r/a + 1)**(2+gamma)) + hyp2 / (gamma*(r/a)**gamma)
        sigma2_dim_less = prefac1 * prefac2 * fac
        return sigma2_dim_less * (self.cosmo.arcsec2phys_lens(1.) * const.Mpc / 1000)**2
1299877
def D_dt(self, H_0, Om0, Ode0=None):
        
        lensCosmo = self._get_cosom(H_0, Om0, Ode0)
        return lensCosmo.D_dt
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
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)
302153
def funTransSectY(theta, rho):
    
    return np.matrix([[1, rho * theta], [0, 1]], dtype=np.double)
1331154
def moduli_to_velocities(rho, K_s, G):
    
    return np.sqrt(K_s / rho), np.sqrt(G / rho)
1186889
def real_rotated_Gammas(m,l,phi1,phi2,theta1,theta2,gamma_ml):
    
    
    if m>0:
        ans=(1./sqrt(2))*(rotated_Gamma_ml(m,l,phi1,phi2,theta1,theta2,gamma_ml) + \
                          (-1)**m*rotated_Gamma_ml(-m,l,phi1,phi2,theta1,theta2,gamma_ml))
        return ans.real
    if m==0:
        return rotated_Gamma_ml(0,l,phi1,phi2,theta1,theta2,gamma_ml).real
    if m<0:
        ans=(1./sqrt(2)/complex(0.,1))*(rotated_Gamma_ml(-m,l,phi1,phi2,theta1,theta2,gamma_ml) - \
                                        (-1)**m*rotated_Gamma_ml(m,l,phi1,phi2,theta1,theta2,gamma_ml))
        return ans.real
391709
def _compute_attenuation(self, rup, dists, imt, C):
        

        vec = np.ones(len(dists.rrup))

        a1 = (np.log10(np.sqrt(dists.rrup ** 2.0 + C['c11'] ** 2.0)),
              np.log10(70. * vec))

        a = np.column_stack([a1[0], a1[1]])

        b3 = (np.log10(np.sqrt(dists.rrup ** 2.0 + C['c11'] ** 2.0) /
                      (70. * vec)),
              np.log10((140. / 70.) * vec))

        b2 = np.column_stack([b3[0], b3[1]])
        b1 = ([np.min(b2, axis=1), 0. * vec])
        b = np.column_stack([b1[0], b1[1]])

        c1 = (np.log10(np.sqrt(dists.rrup ** 2.0 + C['c11'] ** 2.0) /
              (140.) * vec), 0. * vec)
        c = np.column_stack([c1[0], c1[1]])

        return (((C['c4'] + C['c5'] * rup.mag) * np.min(a, axis=1)) +
                ((C['c6'] + C['c7'] * rup.mag) * np.max(b, axis=1)) +
                ((C['c8'] + C['c9'] * rup.mag) * np.max(c, axis=1)))
128702
def fRD( a, M):
    """"""
    f = (lal.C_SI**3.0 / (2.0*lal.PI*lal.G_SI*M*lal.MSUN_SI)) * (1.5251 - 1.1568*(1.0-a)**0.1292)
    return f
789882
def meff_e_Gamma(self, **kwargs):
        '''
        
        '''
        Eg = self.Eg_Gamma(**kwargs)
        Delta_SO = self.Delta_SO(**kwargs)
        Ep = self.Ep(**kwargs)
        F = self.F(**kwargs)
        return 1./((1.+2.*F)+(Ep*(Eg+2.*Delta_SO/3.))/(Eg*(Eg+Delta_SO)))