Compute and return e2 factor, equation 19, page 80.
434716
def _compute_e2_factor(self, imt, vs30):
        
        e2 = np.zeros_like(vs30)

        if imt.name == "PGV":
            period = 1
        elif imt.name == "PGA":
            period = 0
        else:
            period = imt.period

        if period < 0.35:
            return e2
        else:
            idx = vs30 <= 1000
            if period >= 0.35 and period <= 2.0:
                e2[idx] = (-0.25 * np.log(vs30[idx] / 1000) *
                           np.log(period / 0.35))
            elif period > 2.0:
                e2[idx] = (-0.25 * np.log(vs30[idx] / 1000) *
                           np.log(2.0 / 0.35))
            return e2
166040
def _t_of_e(self, a0=None, t_start=None, f0=None, ef=None, t_obs=5.0):
        
        if ef is None:
            ef = np.ones_like(self.e0)*0.0000001

        beta = 64.0/5.0*self.m1*self.m2*(self.m1+self.m2)

        e_vals = np.asarray([np.linspace(ef[i], self.e0[i], self.num_points)
                            for i in range(len(self.e0))])
        integrand = self._find_integrand(e_vals)
        integral = np.asarray([np.trapz(integrand[:, i:], x=e_vals[:, i:])
                              for i in range(e_vals.shape[1])]).T

        if a0 is None and f0 is None:

            a0 = (19./12.*t_start*beta*1/integral[:, 0])**(1./4.) * self._f_e(e_vals[:, -1])

        elif a0 is None:
            a0 = ((self.m1 + self.m2)/self.f0**2)**(1./3.)

        c0 = self._c0_func(a0, self.e0)

        a_vals = c0[:, np.newaxis]*self._f_e(e_vals)

        delta_t = 12./19*c0[:, np.newaxis]**4/beta[:, np.newaxis]*integral

        return e_vals, a_vals, delta_t
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
154168
function e7() {
  const d0 = (Math.random() * 0xffffffff) | 0;
  const d1 = (Math.random() * 0xffffffff) | 0;
  const d2 = (Math.random() * 0xffffffff) | 0;
  const d3 = (Math.random() * 0xffffffff) | 0;

  return `${lut[d0 & 0xff] +
  lut[(d0 >> 8) & 0xff] +
  lut[(d0 >> 16) & 0xff] +
  lut[(d0 >> 24) & 0xff]}-${lut[d1 & 0xff]}${lut[(d1 >> 8) & 0xff]}-${lut[
    ((d1 >> 16) & 0x0f) | 0x40
  ]}${lut[(d1 >> 24) & 0xff]}-${lut[(d2 & 0x3f) | 0x80]}${lut[
    (d2 >> 8) & 0xff
  ]}-${lut[(d2 >> 16) & 0xff]}${lut[(d2 >> 24) & 0xff]}${lut[d3 & 0xff]}${lut[
    (d3 >> 8) & 0xff
  ]}${lut[(d3 >> 16) & 0xff]}${lut[(d3 >> 24) & 0xff]}`;
}
1113853
def _c0(self):
        ""
        h_df = self.df / 2
        p, S = self._p, self.S

        return h_df * (logdet(S) + p * logtwo) + lpgamma(p, h_df)
507296
def calc_cR(Q2, sigma):
    """"""
    return Q2 * np.exp(np.sum(np.log(sigma**2))/sigma.shape[0])
225101
def obfn_g0(self, Y0):
        
        """

        return np.linalg.norm(sp.proj_l2(Y0, self.epsilon,
                                         axis=self.cri.axisN) - Y0)
1094826
def get_step_f(step_f, lR2, lS2):
    
    mu, tau = 10, 2
    if lR2 > mu*lS2:
        return step_f * tau
    elif lS2 > mu*lR2:
        return step_f / tau
    return step_f
1320778
def calculate_iI_correspondence(omega):
    """
    Ne = len(omega[0])
    om = omega[0][0]
    correspondence = []
    I = 0
    for i in range(Ne):
        if omega[i][0] != om:
            om = omega[i][0]
            I += 1
        correspondence += [(i+1, I+1)]
    Nnd = I+1

    def I_nd(i):
        return correspondence[i-1][1]

    def i_d(I):
        for i in range(Ne):
            if correspondence[i][1] == I:
                return correspondence[i][0]

    return i_d, I_nd, Nnd
421392
def e2dnde_deriv(self, x, params=None):
        """"""
        params = self.params if params is None else params
        return np.squeeze(self.eval_e2dnde_deriv(x, params, self.scale,
                                                 self.extra_params))
1310726
def Exponential(x, a, tau, y0):
    
    return np.exp(x / tau) * a + y0
442663
def _g(self, h, xp, s):
        """"""

        nphi = sum(self.phi)

        return (nphi / 2.0) * log(2 * pi) + nphi * \
            log(s) + 0.5 * sum((h - xp) ** 2) / (s ** 2)
144133
def sigma_s2(self, r, R, r_ani, a, gamma, rho0_r0_gamma):
        
        beta = self._beta_ani(r, r_ani)
        return (1 - beta * R**2/r**2) * self.sigma_r2(r, a, gamma, rho0_r0_gamma, r_ani)
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
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)
768774
def _compute_v1_factor(self, imt):
        
        if imt.name == "SA":
            t = imt.period
            if t <= 0.50:
                v1 = 1500.0
            elif t > 0.50 and t <= 1.0:
                v1 = np.exp(8.0 - 0.795 * np.log(t / 0.21))
            elif t > 1.0 and t < 2.0:
                v1 = np.exp(6.76 - 0.297 * np.log(t))
            else:
                v1 = 700.0
        elif imt.name == "PGA":
            v1 = 1500.0
        else:
            # this is for PGV
            v1 = 862.0

        return v1
764295
def obfn_dfd(self):
        
        """

        Ef = self.eval_Rf(self.Xf)
        return sl.rfl2norm2(Ef, self.S.shape, axis=self.cri.axisN) / 2.0