Compute distance scaling term (eq.3, page 319).

        The distance scaling assumes the near-source effect of local site
        conditions due to 50% very firm soil and soft rock and 50% firm rock.
485284
def _compute_distance_scaling(self, C, mag, rrup):
        
        g = C['c5'] + C['c6'] * 0.5 + C['c7'] * 0.5

        return (
            rrup ** 2 +
            (np.exp(C['c8'] * mag + C['c9'] * (8.5 - mag) ** 2) * g) ** 2
        )
1037588
def _compute_non_linear_term(self, C, pga_only, sites):
        
        Vref = 750.0
        Vcon = 1000.0
        lnS = np.zeros_like(sites.vs30)

        # equation (3a)
        idx = sites.vs30 < Vref
        lnS[idx] = (
            C['b1'] * np.log(sites.vs30[idx] / Vref) +
            C['b2'] * np.log(
                (pga_only[idx] + C['c'] * (sites.vs30[idx] / Vref) ** C['n']) /
                ((pga_only[idx] + C['c']) * (sites.vs30[idx] / Vref) ** C['n'])
            )
        )

        # equation (3b)
        idx = (sites.vs30 >= Vref) & (sites.vs30 <= Vcon)
        lnS[idx] = C['b1'] * np.log(sites.vs30[idx]/Vref)

        # equation (3c)
        idx = sites.vs30 > Vcon
        lnS[idx] = C['b1'] * np.log(Vcon/Vref)

        return lnS
562491
def get_magnitude_scaling_term(self, C, rup):
        
        if rup.ztor > 25.0:
            # Deep interface events
            c_int = C["cint"]
        else:
            c_int = C["cintS"]

        if rup.mag <= self.CONSTANTS["m_c"]:
            return c_int * rup.mag
        else:
            return (c_int * self.CONSTANTS["m_c"]) +\
                (C["dint"] * (rup.mag - self.CONSTANTS["m_c"]))
438733
def _get_mean_rock(self, mag, _rake, rrup, is_reverse, imt):
        
        if mag <= self.NEAR_FIELD_SATURATION_MAG:
            C = self.COEFFS_ROCK_LOWMAG[imt]
        else:
            C = self.COEFFS_ROCK_HIMAG[imt]
        # clip mag if greater than 8.5. This is to avoid
        # ValueError: negative number cannot be raised to a fractional power
        mag = 8.5 if mag > 8.5 else mag
        mean = (
            C['c1'] + C['c2'] * mag + C['c3'] * ((8.5 - mag) ** 2.5)
            + C['c4'] * numpy.log(rrup + numpy.exp(C['c5'] + C['c6'] * mag))
            + C['c7'] * numpy.log(rrup + 2)
        )
        if is_reverse:
            # footnote in table 2 says that for reverse ruptures
            # the mean amplitude value should be multiplied by 1.2
            mean += 0.1823215567939546  # == log(1.2)
        return mean
1043306
def get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types):
        
        # extracting dictionary of coefficients specific to required
        # intensity measure type.
        C = self.COEFFS[imt]
        C_SITE = self.SITE_COEFFS[imt]
        s_c, idx = self._get_site_classification(sites.vs30)
        sa_rock = (self.get_magnitude_scaling_term(C, rup) +
                   self.get_sof_term(C, rup) +
                   self.get_depth_term(C, rup) +
                   self.get_distance_term(C, dists, rup))

        sa_soil = self.add_site_amplification(C, C_SITE, sites,
                                              sa_rock, idx, rup)

        stddevs = self.get_stddevs(C, sites.vs30.shape, idx, stddev_types)
        return sa_soil, stddevs
839774
def _compute_large_distance_term(self, C, dists, rup):
        
        # equation 15, page 79
        if rup.mag < 5.5:
            T6 = 1.0
        elif rup.mag >= 5.5 and rup.mag <= 6.5:
            T6 = 0.5 * (6.5 - rup.mag) + 0.5
        else:
            T6 = 0.5

        # equation 14, page 79
        large_distance_term = np.zeros_like(dists.rrup)
        idx = dists.rrup >= 100.0
        large_distance_term[idx] = C['a18'] * (dists.rrup[idx] - 100.0) * T6

        return large_distance_term
961447
def get_magnitude_scaling(self, C, mag):
        
        d_m = mag - self.CONSTANTS["Mh"]
        if mag < self.CONSTANTS["Mh"]:
            return C["e1"] + C["b1"] * d_m + C["b2"] * (d_m ** 2.0)
        else:
            return C["e1"] + C["b3"] * d_m
1243087
def _get_shallow_site_response_term(self, C, vs30, pga_rock):
        
        vs_mod = vs30 / C["k1"]
        # Get linear global site response term
        f_site_g = C["c11"] * np.log(vs_mod)
        idx = vs30 > C["k1"]
        f_site_g[idx] = f_site_g[idx] + (C["k2"] * self.CONSTS["n"] *
                                         np.log(vs_mod[idx]))

        # Get nonlinear site response term
        idx = np.logical_not(idx)
        if np.any(idx):
            f_site_g[idx] = f_site_g[idx] + C["k2"] * (
                np.log(pga_rock[idx] +
                       self.CONSTS["c"] * (vs_mod[idx] ** self.CONSTS["n"])) -
                np.log(pga_rock[idx] + self.CONSTS["c"])
                )

        # For Japan sites (SJ = 1) further scaling is needed (equation 19)
        if self.CONSTS["SJ"]:
            fsite_j = np.log(vs_mod)
            idx = vs30 > 200.0
            if np.any(idx):
                fsite_j[idx] = (C["c13"] + C["k2"] * self.CONSTS["n"]) *\
                    fsite_j[idx]
            idx = np.logical_not(idx)
            if np.any(idx):
                fsite_j[idx] = (C["c12"] + C["k2"] * self.CONSTS["n"]) *\
                    (fsite_j[idx] - np.log(200.0 / C["k1"]))

            return f_site_g + fsite_j
        else:
            return f_site_g
945385
def _compute_magnitude_scaling(self, mag, C):
        
        dmag = mag - 6.
        return (C['B2'] * dmag) + (C['B3'] * (dmag ** 2.))
931164
def _compute_pga_rock(self, C, dc1, sites, rup, dists):
        
        mean = (self._compute_magnitude_term(C, dc1, rup.mag) +
                self._compute_distance_term(C, rup.mag, dists) +
                self._compute_focal_depth_term(C, rup) +
                self._compute_forearc_backarc_term(C, sites, dists))
        # Apply linear site term
        site_response = ((C['theta12'] + C['b'] * self.CONSTS['n']) *
                         np.log(1000. / C['vlin']))
        return mean + site_response
951229
def _get_sa_at_1180(self, C, imt, sites, rup, dists):
        
        # reference vs30 = 1180 m/s
        vs30_1180 = np.ones_like(sites.vs30) * 1180.
        # reference shaking intensity = 0
        ref_iml = np.zeros_like(sites.vs30)
        # fake Z1.0 - Since negative it will be replaced by the default Z1.0
        # for the corresponding region
        fake_z1pt0 = np.ones_like(sites.vs30) * -1
        return (self._get_basic_term(C, rup, dists) +
                self._get_faulting_style_term(C, rup) +
                self._get_site_response_term(C, imt, vs30_1180, ref_iml) +
                self._get_hanging_wall_term(C, dists, rup) +
                self._get_top_of_rupture_depth_term(C, imt, rup) +
                self._get_soil_depth_term(C, fake_z1pt0, vs30_1180) +
                self._get_regional_term(C, imt, vs30_1180, dists.rrup)
                )
505482
def _compute_hanging_wall_effect(self, C, rjb, rrup, dip, mag):
        
        # eq. 8 (to be noticed that the USGS-NSHMP implementation defines
        # the hanging-wall term for all rjb distances, while in the original
        # manuscript, hw is computed only for rjb < 5). Again the 'firm rock'
        # is considered
        hw = np.zeros_like(rjb)
        if dip <= 70.:
            hw = (5. - rjb) / 5.

        # eq. 9
        f_m = 1 if mag > 6.5 else mag - 5.5

#        # eq. 10
        f_rrup = C['c15'] + np.zeros_like(rrup)
        idx = rrup < 8
        f_rrup[idx] *= rrup[idx] / 8

        # eq. 7 (to be noticed that the f3 factor is not included
        # while this is defined in the original manuscript)
        f_hw = hw * f_m * f_rrup

        return f_hw
708026
def _compute_mean(self, C, g, mag, hypo_depth, rrup, vs30, pga_rock, imt):
        
        if hypo_depth > 100:
            hypo_depth = 100
        delta = 0.00724 * 10 ** (0.507 * mag)
        R = np.sqrt(rrup ** 2 + delta ** 2)

        s_amp = self._compute_soil_amplification(C, vs30, pga_rock, imt)

        mean = (
            # 1st term
            C['c1'] + C['c2'] * mag +
            # 2nd term
            C['c3'] * hypo_depth +
            # 3rd term
            C['c4'] * R -
            # 4th term
            g * np.log10(R) +
            # 5th, 6th and 7th terms
            s_amp
        )

        return mean
425888
def get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types):
        
        # extracting dictionary of coefficients specific to required
        # intensity measure type.
        C = self.COEFFS[imt]
        C_PGA = self.COEFFS[PGA()]
        imt_per = 0 if imt.name == 'PGV' else imt.period
        pga_rock = self._get_pga_on_rock(C_PGA, rup, dists)
        mean = (self._get_magnitude_scaling_term(C, rup) +
                self._get_path_scaling(C, dists, rup.mag) +
                self._get_site_scaling(C, pga_rock, sites, imt_per, dists.rjb))
        stddevs = self._get_stddevs(C, rup, dists, sites, stddev_types)
        return mean, stddevs
1290161
def _compute_site_response_term(self, C, sites, pga1000):
        
        vs_star = sites.vs30.copy()
        vs_star[vs_star > 1000.0] = 1000.
        arg = vs_star / C["vlin"]
        site_resp_term = C["theta12"] * np.log(arg)
        # Get linear scaling term
        idx = sites.vs30 >= C["vlin"]
        site_resp_term[idx] += (C["b"] * self.CONSTS["n"] * np.log(arg[idx]))
        # Get nonlinear scaling term
        idx = np.logical_not(idx)
        site_resp_term[idx] += (
            -C["b"] * np.log(pga1000[idx] + self.CONSTS["c"]) +
            C["b"] * np.log(pga1000[idx] + self.CONSTS["c"] *
                            (arg[idx] ** self.CONSTS["n"])))
        return site_resp_term
17025
def _get_distance_scaling_term(self, C, rval, mag):
        
        r_adj = np.sqrt(rval ** 2.0 + C["h"] ** 2.0)
        return (
            (C["c1"] + C["c2"] * (mag - self.CONSTS["Mref"])) *
            np.log10(r_adj / self.CONSTS["Rref"]) -
            (C["c3"] * (r_adj - self.CONSTS["Rref"])))
855805
def _get_pga_on_rock(self, rup, dists, _C):
        
        # Median PGA in g for Vref = 760.0, without site amplification,
        # that is equation (1) pag 106, without the third and fourth terms
        # Mref and Rref values are given in the caption to table 6, pag 119
        # Note that in the original paper, the caption reads:
        # "Distance-scaling coefficients (Mref=4.5 and Rref=1.0 km for all
        # periods, except Rref=5.0 km for pga4nl)". However this is a mistake
        # as reported in http://www.daveboore.com/pubs_online.php:
        # ERRATUM: 27 August 2008. Tom Blake pointed out that the caption to
        # Table 6 should read "Distance-scaling coefficients (Mref=4.5 and
        # Rref=1.0 km for all periods)".
        C_pga = self.COEFFS[PGA()]
        pga4nl = np.exp(self._compute_magnitude_scaling(rup, C_pga) +
                        self._compute_distance_scaling(rup, dists, C_pga))

        return pga4nl