Compute the second term of the equation described on p. 199:

        ``(b4 + b5 * M) * log(sqrt(Rjb ** 2 + b6 ** 2))``
427514
def _compute_distance(self, rup, dists, imt, C):
        
        return (((C['b4'] + C['b5'] * rup.mag)
                 * np.log10((np.sqrt(dists.rjb ** 2.0 + C['b6'] ** 2.0)))))
1085985
def _compute_mean(self, C, mag, r):
        
        mean = (C['c1'] +
                self._compute_term1(C, mag) +
                self._compute_term2(C, mag, r))
        return mean
816928
def _compute_mean(self, C, mag, rjb):
        
        d1 = np.sqrt(50. ** 2 + 6. ** 2)
        d = np.sqrt(rjb ** 2 + 6 ** 2)

        mean = np.zeros_like(rjb)

        mean += (
            C['a1'] + C['a2'] * (mag - 6.4) +
            C['a7'] * (8.5 - mag) ** 2
        )

        idx = rjb < 50.
        mean[idx] += (
            C['a3'] * np.log(d[idx]) +
            C['a4'] * (mag - 6.4) * np.log(d[idx]) +
            C['a5'] * rjb[idx]
        )

        idx = rjb >=50.
        mean[idx] += (
            C['a3'] * np.log(d1) +
            C['a4'] * (mag - 6.4) * np.log(d[idx]) +
            C['a5'] * rjb[idx] + C['a6'] * (np.log(d[idx]) - np.log(d1))
        )

        return mean
827602
def _get_magnitude_term(self, C, mag):
        
        if mag >= self.CONSTS["Mh"]:
            return C["e1"] + C["b3"] * (mag - self.CONSTS["Mh"])
        else:
            return C["e1"] + (C["b1"] * (mag - self.CONSTS["Mh"])) +\
                (C["b2"] * (mag - self.CONSTS["Mh"]) ** 2.)
576122
def _get_magnitude_term(self, C, mag):
        
        f_mag = C["c0"] + C["c1"] * mag
        if (mag > 4.5) and (mag <= 5.5):
            return f_mag + (C["c2"] * (mag - 4.5))
        elif (mag > 5.5) and (mag <= 6.5):
            return f_mag + (C["c2"] * (mag - 4.5)) + (C["c3"] * (mag - 5.5))
        elif mag > 6.5:
            return f_mag + (C["c2"] * (mag - 4.5)) + (C["c3"] * (mag - 5.5)) +\
                (C["c4"] * (mag - 6.5))
        else:
            return f_mag
965000
def _fit_s(D, a0, logp, tol=1e-7, maxiter=1000):
    ''''''
    N, K = D.shape
    s1 = a0.sum()
    m = a0 / s1
    mlogp = (m*logp).sum()
    for i in xrange(maxiter):
        s0 = s1
        g = psi(s1) - (m*psi(s1*m)).sum() + mlogp
        h = _trigamma(s1) - ((m**2)*_trigamma(s1*m)).sum()

        if g + s1 * h < 0:
            s1 = 1/(1/s0 + g/h/(s0**2))
        if s1 <= 0:
            s1 = s0 * exp(-g/(s0*h + g)) # Newton on log s
        if s1 <= 0:
            s1 = 1/(1/s0 + g/((s0**2)*h + 2*s0*g)) # Newton on 1/s
        if s1 <= 0:
            s1 = s0 - g/h # Newton
        if s1 <= 0:
            raise Exception('Unable to update s from {}'.format(s0))

        a = s1 * m
        if abs(s1 - s0) < tol:
            return a

    raise Exception('Failed to converge after {} iterations, s is {}'
            .format(maxiter, s1))
369270
def _get_basic_term(self, C, rup, dists):
        
        # Fictitious depth calculation
        if rup.mag > 5.:
            c4m = C['c4']
        elif rup.mag > 4.:
            c4m = C['c4'] - (C['c4']-1.) * (5. - rup.mag)
        else:
            c4m = 1.
        R = np.sqrt(dists.rrup**2. + c4m**2.)
        # basic form
        base_term = C['a1'] * np.ones_like(dists.rrup) + C['a17'] * dists.rrup
        # equation 2 at page 1030
        if rup.mag >= C['m1']:
            base_term += (C['a5'] * (rup.mag - C['m1']) +
                          C['a8'] * (8.5 - rup.mag)**2. +
                          (C['a2'] + C['a3'] * (rup.mag - C['m1'])) *
                          np.log(R))
        elif rup.mag >= self.CONSTS['m2']:
            base_term += (C['a4'] * (rup.mag - C['m1']) +
                          C['a8'] * (8.5 - rup.mag)**2. +
                          (C['a2'] + C['a3'] * (rup.mag - C['m1'])) *
                          np.log(R))
        else:
            base_term += (C['a4'] * (self.CONSTS['m2'] - C['m1']) +
                          C['a8'] * (8.5 - self.CONSTS['m2'])**2. +
                          C['a6'] * (rup.mag - self.CONSTS['m2']) +
                          C['a7'] * (rup.mag - self.CONSTS['m2'])**2. +
                          (C['a2'] + C['a3'] * (self.CONSTS['m2'] - C['m1'])) *
                          np.log(R))
        return base_term
1284626
def compute(self, x, y):
        '''
        
        '''
        return np.sqrt(np.sum((x-y)**2))
246238
def _compute_mean(self, C, mag, rjb):
        
        m1 = 6.4
        r1 = 50.
        h = 6.
        R = np.sqrt(rjb ** 2 + h ** 2)
        R1 = np.sqrt(r1 ** 2 + h ** 2)
        less_r1 = rjb < r1
        ge_r1 = rjb >= r1

        mean = (C['c1'] + C['c4'] * (mag - m1) * np.log(R) + C['c5'] * rjb +
                C['c8'] * (8.5 - mag) ** 2)

        mean[less_r1] += C['c3'] * np.log(R[less_r1])
        mean[ge_r1] += (C['c3'] * np.log(R1) +
                        C['c6'] * (np.log(R[ge_r1]) - np.log(R1)))

        if mag < m1:
            mean += C['c2'] * (mag - m1)
        else:
            mean += C['c7'] * (mag - m1)

        return mean
1065723
def _compute_f4(self, C, mag, rrup):
        
        fhw_m = 0
        fhw_r = np.zeros_like(rrup)

        if mag <= 5.5:
            fhw_m = 0
        elif 5.5 < mag < 6.5:
            fhw_m = mag - 5.5
        else:
            fhw_m = 1

        idx = (rrup > 4) & (rrup <= 8)
        fhw_r[idx] = C['a9'] * (rrup[idx] - 4.) / 4.

        idx = (rrup > 8) & (rrup <=18)
        fhw_r[idx] = C['a9']

        idx = (rrup > 18) & (rrup <= 24)
        fhw_r[idx] = C['a9'] * (1 - (rrup[idx] - 18.) / 7.)

        return fhw_m * fhw_r
475784
def _compute_magnitude(self, rup, C):
        
        return C['b1'] + (C['b2'] * rup.mag) + (C['b3'] * (rup.mag ** 2))
142996
def _compute_distance_term(self, C, rup, dists):
        
        return (C['c4'] + C['c5'] * rup.mag) * \
            np.log(np.sqrt(dists.rrup ** 2. + C['c6'] ** 2.))
867563
function loggamma(x) {
        var v = 1.0,
            w = 1.0 / (x * x),
            ret;
        while (x < N) {
            v = v * x;
            x = x + 1;
        }
        ret = B16 / (16 * 15);
        ret = ret * w + B14 / (14 * 13);
        ret = ret * w + B12 / (12 * 11);
        ret = ret * w + B10 / (10 *  9);
        ret = ret * w + B8  / (8 *  7);
        ret = ret * w + B6  / (6 *  5);
        ret = ret * w + B4  / (4 *  3);
        ret = ret * w + B2  / (2 *  1);
        ret = ret / x + 0.5 * LOG_2PI - Math.log(v) - x + (x - 0.5) * Math.log(x);
        return ret;
    }
1270101
def j0_4(f=1, rho=0.3, z=50):
    

    """

    gam = np.sqrt(2j*np.pi*mu_0*f/rho)

    def lhs(x):
        beta = np.sqrt(x**2 + gam**2)
        return x*np.exp(-beta*np.abs(z))/beta

    def rhs(b):
        R = np.sqrt(b**2 + z**2)
        return np.exp(-gam*R)/R

    return Ghosh('j0', lhs, rhs)
343163
def _g(self, x, tau):
        
        return tau ** 2 * (tau ** 2 + 1) ** -2 * (
                (tau ** 2 + 1 + 2 * (x ** 2 - 1)) * self.F(x) + tau * np.pi + (tau ** 2 - 1) * np.log(tau) +
                np.sqrt(tau ** 2 + x ** 2) * (-np.pi + self.L(x, tau) * (tau ** 2 - 1) * tau ** -1))
172139
def _F(x, gam):
        
        beta = np.sqrt(1 - gam ** -2)
        B = 1 + 0.5 * (gam ** 2 - 1)
        C = 10 * x * gam * beta * (2 + gam * beta)
        C /= 1 + x ** 2 * (gam ** 2 - 1)

        F_1 = (17 - 3 * x ** 2 / (2 - x) ** 2 - C) * np.sqrt(1 - x)
        F_2 = 12 * (2 - x) - 7 * x ** 2 / (2 - x) - 3 * x ** 4 / (2 - x) ** 3
        F_3 = np.log((1 + np.sqrt(1 - x)) / np.sqrt(x))

        return B * F_1 + F_2 * F_3
1222606
def xyz2rgb(x__, y__, z__):
    
    x2_ = x__ / 100.0
    y2_ = y__ / 100.0
    z2_ = z__ / 100.0

    r__ = x2_ *  3.2406 + y2_ * -1.5372 + z2_ * -0.4986
    g__ = x2_ * -0.9689 + y2_ *  1.8758 + z2_ *  0.0415
    b__ = x2_ *  0.0557 + y2_ * -0.2040 + z2_ *  1.0570

    def finv(arr):
        """Inverse
        """
        return np.where(arr > 0.0031308,
                        1.055 * (arr ** (1.0 / 2.4)) - 0.055,
                        12.92 * arr)

    return finv(r__) * 255, finv(g__) * 255, finv(b__) * 255