equation A20 in Eliasdottir (2013)
        :param r_a: r/Ra
        :param r_s: r/Rs
        :return:
771954
def _f_A20(self, r_a, r_s):
        
        return r_a/(1+np.sqrt(1 + r_a**2)) - r_s/(1+np.sqrt(1 + r_s**2))
1049434
def _h_ars(x,params):
    """"""
    mO, sO2= params
    return -0.5*(x-mO)**2./sO2+numpy.log(x)
428380
def a2s(a):
    
    s = np.zeros((6,), 'f')  # make the a matrix
    for i in range(3):
        s[i] = a[i][i]
    s[3] = a[0][1]
    s[4] = a[1][2]
    s[5] = a[0][2]
    return s
1277988
def rsdl_r(self, AX, Y):
        

        if not hasattr(self, '_cnst_c0'):
            self._cnst_c0 = self.cnst_c0()
        if not hasattr(self, '_cnst_c1'):
            self._cnst_c1 = self.cnst_c1()
        return AX - self.block_cat(self.var_y0() + self._cnst_c0,
                                   self.var_y1() + self._cnst_c1)
391824
def _erfc(a):
    
    # approximation for abs(a) < 8 and abs(a) >= 1
    P = [
        2.46196981473530512524E-10,
        5.64189564831068821977E-1,
        7.46321056442269912687E0,
        4.86371970985681366614E1,
        1.96520832956077098242E2,
        5.26445194995477358631E2,
        9.34528527171957607540E2,
        1.02755188689515710272E3,
        5.57535335369399327526E2,
    ]

    Q = [
        1.32281951154744992508E1,
        8.67072140885989742329E1,
        3.54937778887819891062E2,
        9.75708501743205489753E2,
        1.82390916687909736289E3,
        2.24633760818710981792E3,
        1.65666309194161350182E3,
        5.57535340817727675546E2,
    ]

    # approximation for abs(a) >= 8
    R = [
        5.64189583547755073984E-1,
        1.27536670759978104416E0,
        5.01905042251180477414E0,
        6.16021097993053585195E0,
        7.40974269950448939160E0,
        2.97886665372100240670E0,
    ]

    S = [
        2.26052863220117276590E0,
        9.39603524938001434673E0,
        1.20489539808096656605E1,
        1.70814450747565897222E1,
        9.60896809063285878198E0,
        3.36907645100081516050E0,
    ]

    # Shortcut special cases
    if a == 0:
        return 1
    if a >= MAXVAL:
        return 0
    if a <= -MAXVAL:
        return 2

    x = a
    if a < 0:
        x = -a

    # computationally cheaper to calculate erf for small values, I guess.
    if x < 1:
        return 1 - erf(a)

    z = -a * a

    z = math.exp(z)

    if x < 8:
        p = _polevl(x, P, 8)
        q = _p1evl(x, Q, 8)
    else:
        p = _polevl(x, R, 5)
        q = _p1evl(x, S, 6)

    y = (z * p) / q

    if a < 0:
        y = 2 - y

    return y
244980
def DSYR_blas(A, x, alpha=1.):
    
    At = blas.dsyr(lower=0, x=x, a=A, alpha=alpha, overwrite_a=False) #See https://github.com/scipy/scipy/issues/8155
    A[:] = At
    symmetrify(A, upper=True)
1146156
def _compute_mean_on_rock(self, C, mag, rrup, rvol, hypo_depth, CN, CR,
                              f4HW):
        

        lnSA_AB = (
            # line 1 of equation 1
            C['c1'] + C['c4as'] * (mag - 6) +
            # line 2
            C['c3as'] * (8.5 - mag) ** 2 +
            # line 3
            C['c5'] * rrup +
            # line 3 and 4
            (C['c8'] + C['c6as'] * (mag - 6)) *
            np.log((rrup ** 2 + C['c10as'] ** 2) ** 0.5) +
            # line 5
            C['c46'] * rvol +
            # line 6
            C['c32'] * CN + C['c33as'] * CR + f4HW
        )

        return lnSA_AB
205648
def MatrixSolve(a, rhs, adj):
    
    return np.linalg.solve(a if not adj else _adjoint(a), rhs),
927852
def _compute_mean(self, C, mag, rrup, hypo_depth, delta_R, delta_S,
                      delta_V, delta_I, vs30):

        
        # mean is calculated for all the 4 classes using the same equation.
        # For DowrickRhoades2005SSlab, the coefficients which don't appear in
        # Model 3 equationare assigned to zero

        mean = (C['A1'] + (C['A2'] + C['A2R'] * delta_R + C['A2V'] * delta_V) *
                mag + (C['A3'] + C['A3S'] * delta_S + C['A3V'] * delta_V) *
                np.log10(np.power((rrup**3 + C['d']**3), 1.0 / 3.0)) +
                C['A4'] * hypo_depth + C['A5'] * delta_I)

        # Get S site class term
        S = self._get_site_class(vs30, mean)

        # Add S amplification term to mean value
        mean = mean + S

        return mean
1027584
def _update_rho(self, k, r, s):
    

    if self.opt['AutoRho', 'Enabled']:
        tau = self.rho_tau
        mu = self.rho_mu
        xi = self.rho_xi
        if k != 0 and cp.mod(k + 1, self.opt['AutoRho', 'Period']) == 0:
            if self.opt['AutoRho', 'AutoScaling']:
                if s == 0.0 or r == 0.0:
                    rhomlt = tau
                else:
                    rhomlt = cp.sqrt(r / (s * xi) if r > s * xi
                                     else (s * xi) / r)
                    if rhomlt > tau:
                        rhomlt = tau
            else:
                rhomlt = tau
            rsf = 1.0
            if r > xi * mu * s:
                rsf = rhomlt
            elif s > (mu / xi) * r:
                rsf = 1.0 / rhomlt
            self.rho *= float(rsf)
            self.U /= rsf
            if rsf != 1.0:
                self.rhochange()
871138
def mass_3d(self, r, rho0, Rs):
        
        mass_3d = 2*np.pi*Rs**3*rho0 * r**2/(r + Rs)**2
        return mass_3d
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
656918
def mass_tot(self, rho0, Rs):
        
        m_tot = 2*np.pi*rho0*Rs**3
        return m_tot
142694
def sigma_r2(self, r, kwargs_profile, kwargs_anisotropy, kwargs_light):
        
        return self.jeans_solver.sigma_r2(r, kwargs_profile, kwargs_anisotropy, kwargs_light)
494509
def rarlognormal(a, sigma, rho, size=1):
    
    """
    f = utils.ar1
    if np.isscalar(a):
        r = f(rho, 0, sigma, size)
    else:
        n = len(a)
        r = [f(rho, 0, sigma, n) for i in range(size)]
        if size == 1:
            r = r[0]
    return a * np.exp(r)
836551
def j0_2(a=1):
    """

    def lhs(x):
        return np.exp(-a*x)

    def rhs(b):
        return 1/np.sqrt(b**2 + a**2)

    return Ghosh('j0', lhs, rhs)
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