diff --git a/src/main/java/gov/usgs/earthquake/nshmp/gmm/ParkerEtAl_2020.java b/src/main/java/gov/usgs/earthquake/nshmp/gmm/ParkerEtAl_2020.java index 9ba2e5505a7867142780b75bdb58f1f60ebdff4e..5d8eb2a70b5d9d3d27efdc69a2a4448387cb89a3 100644 --- a/src/main/java/gov/usgs/earthquake/nshmp/gmm/ParkerEtAl_2020.java +++ b/src/main/java/gov/usgs/earthquake/nshmp/gmm/ParkerEtAl_2020.java @@ -27,40 +27,30 @@ import gov.usgs.earthquake.nshmp.gmm.GmmInput.Constraints; import gov.usgs.earthquake.nshmp.tree.LogicTree; /** - * Implementation of the PEER NGA-Subduction GMM, Parker et al., 2020 GMM - * (PSBAH20) + * Implementation of the Parker et al. (2022) ground motion model for subduction + * regions developed as part of the <a + * href="https://www.risksciences.ucla.edu/nhr3/nga-subduction">NGA + * Subduction</a> project. * - * Notes: + * <p><b>Implementation notes:</b> * - * This model is only applicable to sites in the forearc, future residuals - * analysis is planned to determine if there are attenuation differences between - * forearc and backarc zones. + * <ul><li>Model currently uses {@code zTor} for intraslab depth scaling but + * authors recommend {@code zHyp}. * - * PSBAH20 uses zHyp (randomized source location on fault), but this can be - * replaced this event-specific zHyp with the mean depth expected for the fault - * plane, zHypBar, in Eq. 4.6 without modification to coefficients or - * uncertainty. + * <ul><li>The basin depth scaling model is that developed for Cascadia.</ul> * - * Path term (Eqs. 4.2 and 4.3 in PEER report) appears to differ in the official - * Python implementation. Compare implementations of getPathTerm() below - * - * Cascadia a0 interface coefficient is set to global values in coeff file. - * - * For this initial implementation we are supporting the Seattle basin model - * only. That is, if z2p5 != NaN, then the basin term is non-zero and uses the - * Cascadia and Seattle specific parameters. - * - * We've switched depth scaling to use zTor, but zHyp is specified by Parker; - * need to check database definitions of these. + * <p><b>Note:</b> Direct instantiation of {@code GroundMotionModel}s is + * prohibited. Use {@link Gmm#instance(Imt)} to retrieve an instance for a + * desired {@link Imt}. * * <p><b>Reference:</b> Parker, G.A, Stewart, J.P., Boore, D.M., Atkinson, G.M., - * and Hassani, B. (2020) NGA-Subduction global ground-motion models with + * and Hassani, B. (2022) NGA-Subduction global ground motion models with * regional adjustment factors, Earthquake Spectra, v. 38, n. 1, p. 456-493, * doi: <a * href="https://doi.org/10.1177/87552930211034889">10.1177/87552930211034889</a>. * * <p><b>Reference:</b> Parker, G.A, Stewart, J.P., Boore, D.M., Atkinson, G.M., - * and Hassani, B. (2020) NGA-Subduction global ground-motion models with + * and Hassani, B. (2020) NGA-Subduction global ground motion models with * regional adjustment factors. Report no. 2020/03. Berkeley, CA: PEER, 131 p. * * <p><b>Component:</b> average horizontal (RotD50) @@ -86,77 +76,82 @@ public abstract class ParkerEtAl_2020 implements GroundMotionModel { .build(); /* + * Developer notes: + * * Regions: global, Alaska, Cascadia, [CAM (Central America), Japan, SA (South * America), Taiwan] * * Saturation Regions: global, Aleutian, Alaska, Cascadia, [Central America S, * Central America N, Japan Pacific, Japan Phillipines, South America N, South * America S, Taiwan W, Taiwan E] + * + * Implementation currently uses zTor for intraslab depth scaling but authors + * regressed on zHyp. Which parameter is used will possibly depend on + * implementation choices for the 2023 NSHM. Historically, intraslab GMMs have + * used zTor (check this) and zHyp makes depths a little deeper, assuming + * pseudo faults have a measureable width, and would yield slightly higher + * ground motions). */ - static final class Basin { - static final String NONE = "None"; - static final String PACNW_NONE = "PacNW-no-basin"; - static final String PACNW = "PacNW"; - static final String SEATTLE = "Seattle"; - } - - // aleatory uncertainty coefficients are duplicated in interface and interface - // CSV files (as are some other coefficients that are also independent of - // event type) - private static final String INTERFACE_CSV = "PSBAH20_interface.csv"; - private static final String SLAB_CSV = "PSBAH20_slab.csv"; - static final CoefficientContainer COEFFS_INTERFACE; static final CoefficientContainer COEFFS_SLAB; static { - COEFFS_INTERFACE = new CoefficientContainer(INTERFACE_CSV); - COEFFS_SLAB = new CoefficientContainer(SLAB_CSV); + /* + * Aleatory uncertainty coefficients are duplicated in interface and + * interface CSV files (as are some other coefficients that are also + * independent of event type) + */ + COEFFS_INTERFACE = new CoefficientContainer("PSBAH20_interface.csv"); + COEFFS_SLAB = new CoefficientContainer("PSBAH20_slab.csv"); } - // private static final double B4 = 0.1; - // private static final double V1 = 270.0; // m/s - // private static final double VREF = 760.0; // m/s + private static final double B4 = 0.1; + private static final double V1 = 270.0; // m/s + private static final double VREF = 760.0; // m/s + private static final double DB1 = 20.0; // km private static final double DB2 = 67.0; // km /* ϕ parameters */ - private static final double R1 = 200; // km - private static final double R2 = 500; // km - private static final double V1 = 200; // m/s - private static final double V2 = 500; // m/s - /* ϕS2S and ϕSS parameters */ - private static final double R3 = 200; // km - private static final double R4 = 500; // km - private static final double R5 = 500; // km - private static final double R6 = 800; // km - private static final double V3 = 200; // m/s - private static final double V4 = 800; // m/s + private static final double ΦR1 = 200; // km + private static final double ΦR2 = 500; // km + private static final double ΦV1 = 200; // m/s + private static final double ΦV2 = 500; // m/s private static final class Coefficients { final Imt imt; final double c0; // model constant - final double c1, b4, a0; // path scaling + final double c1, a0; // path scaling final double mc; // saturation corner magnitude final double c4, c5, c6; // magnitude scaling - final double d, m, db; // source depth scaling (slab only) - final double v1, v2, vRef, s1, s2; // linear site amplification + final double d, m; // source depth scaling (slab only) + final double v2, s1, s2; // linear site amplification final double f4, f5; // non-linear site amplification - final double c_e1, c_e2, c_e3, del_None, del_Seattle; // Cascadia basin - final double τ, ϕ21, ϕ22, ϕ2v, vm, ϕ2s2s0, a1, ϕ2ss1, ϕ2ss2, a2; + final double c_e1, c_e2, c_e3; // Cascadia basin + final double τ, ϕ21, ϕ22, ϕ2v; final double σε1, σε2, t1, t2; + // inlined + // final double b4, v1, vRef; + + // unused + // final double db; + // final del_None, del_Seattle; + // final double vm, ϕ2s2s0, a1, ϕ2ss1, ϕ2ss2, a2; + Coefficients(Imt imt, CoefficientContainer cc, SubductionZone zone, boolean slab) { this.imt = imt; // Required for σε Map<String, Double> coeffs = cc.get(imt); + // Model constant c0 = coeffs.get(zone + "_c0"); + // Path c1 = coeffs.get("c1"); - b4 = coeffs.get("b4"); // period-independent + // b4 inlined as constant // a0 is region-specific for Alaska and, for slab events, Cascadia, // otherwise use Global value @@ -181,15 +176,14 @@ public abstract class ParkerEtAl_2020 implements GroundMotionModel { // Source depth scaling (slab only) d = slab ? coeffs.get("d") : 0.0; m = slab ? coeffs.get("m") : 0.0; - db = slab ? coeffs.get("db") : 0.0; // Linear site amplification - v1 = coeffs.get("V1"); + // v1 inlined as constant v2 = coeffs.get("V2");; - vRef = coeffs.get("Vref"); + // vRef inlined as constant - // s2 is region specific, except for Aleutians (not currently supported), - // which uses the Global value + // s2 is region specific, except for Aleutians (not + // currently supported), which uses the Global value s2 = coeffs.get(zone + "_s2"); s1 = s2; // s1 = s2 except for Japan and Taiwan @@ -201,20 +195,12 @@ public abstract class ParkerEtAl_2020 implements GroundMotionModel { c_e1 = coeffs.get("C_e1"); c_e2 = coeffs.get("C_e2"); c_e3 = coeffs.get("C_e3"); - del_None = coeffs.get("del_None"); - del_Seattle = coeffs.get("del_Seattle"); // aleatory model τ = coeffs.get("Tau"); ϕ21 = coeffs.get("phi21"); ϕ22 = coeffs.get("phi22"); ϕ2v = coeffs.get("phi2V"); - vm = coeffs.get("VM"); - ϕ2s2s0 = coeffs.get("phi2S2S0"); - a1 = coeffs.get("a1"); - ϕ2ss1 = coeffs.get("phi2SS1"); - ϕ2ss2 = coeffs.get("phi2SS2"); - a2 = coeffs.get("a2"); // epistemic model, Table E4 if (zone == CASCADIA) { @@ -251,32 +237,18 @@ public abstract class ParkerEtAl_2020 implements GroundMotionModel { @Override public final LogicTree<GroundMotion> calc(GmmInput in) { - /* - * Equation 4.1 μ_lny = c0 + Fp + Fm + Fd + Fs - */ - // double Mc = getSaturationMagnitude(getRegion(), isSlab()); double pgaRef = exp(calcMean(coeffsPGA, slab(), in.Mw, in.rRup, in.zTor)); - double μ = calcMean( coeffs, slab(), basin(), in.Mw, in.rRup, in.zTor, pgaRef, in.vs30, in.z2p5); - double ϕTotal = getϕTotal(coeffs, in.rRup, in.vs30); + double ϕTotal = phiTotal(coeffs, in.rRup, in.vs30); double σ = Maths.hypot(coeffs.τ, ϕTotal); - - // System.out.println("calc for T = " + coeffs.imt.toString()); - // System.out.printf(" ϕTotal = %.6f\n", ϕTotal); - // System.out.printf(" σ = %.6f\n", σ); - - // double ϕS2S = getϕS2S(coeffs, in.rRup, in.vs30); - // double ϕSS = getϕSS(coeffs, in.rRup, in.vs30); - double σε = getσε(coeffs); double[] μs = new double[] { μ + σε, μ, μ - σε }; - // return GroundMotions.createTree(μ, σ); return GroundMotions.createTree(GroundMotions.EPI_IDS, μs, GroundMotions.EPI_WTS, σ); } @@ -284,7 +256,7 @@ public abstract class ParkerEtAl_2020 implements GroundMotionModel { abstract boolean basin(); - /* calc median at 760 reference conditions */ + /* Calc median at 760 reference condition */ private static double calcMean( Coefficients c, boolean isSlab, @@ -292,23 +264,14 @@ public abstract class ParkerEtAl_2020 implements GroundMotionModel { double rRup, double zHyp) { double h = nearSourceSaturation(Mw, c.mc, isSlab); - double Fp = getPathTermPy(c, rRup, Mw, h); - double Fm = getMagnitudeTerm(c, Mw); - double Fd = getDepthTerm(c, isSlab, zHyp); - - // System.out.println("calcMean Reference for T = " + c.imt.toString()); - // System.out.printf(" mc = %.2f\n", c.mc); - // System.out.printf(" h = %.6f\n", h); - // System.out.printf(" Fp = %.6f\n", Fp); - // System.out.printf(" Fm(ref) = %.6f\n", Fm); - // System.out.printf(" Fd = %.6f\n", Fd); - // System.out.printf(" Fb = %.6f (reference calc)\n", 0.0); - // System.out.printf(" mu = %.6f\n", c.c0 + Fp + Fm + Fd); + double Fp = pathTerm(c, rRup, Mw, h); + double Fm = magnitudeTerm(c, Mw); + double Fd = depthTerm(c, isSlab, zHyp); return c.c0 + Fp + Fm + Fd; } - /* Equation 4.1 */ + /* Equation 1 */ private static double calcMean( Coefficients c, boolean slab, @@ -321,112 +284,83 @@ public abstract class ParkerEtAl_2020 implements GroundMotionModel { double z2p5) { double h = nearSourceSaturation(Mw, c.mc, slab); - double Fp = getPathTermPy(c, rRup, Mw, h); - double Fm = getMagnitudeTerm(c, Mw); - double Fd = getDepthTerm(c, slab, zHyp); + double Fp = pathTerm(c, rRup, Mw, h); + double Fm = magnitudeTerm(c, Mw); + double Fd = depthTerm(c, slab, zHyp); - /* Equation 4.7: Fs = Flin + Fnl + Fb */ - double Fslin = getSiteTermLinear(c, vs30); - double Fsnl = getSiteTermNonLinear(c, pgaRef, vs30); + /* Equation 7 */ + double Fslin = siteTermLinear(c, vs30); + double Fsnl = siteTermNonLinear(c, pgaRef, vs30); double Fsb = basin ? getSiteBasinTerm(c, vs30, z2p5) : 0.0; double Fs = Fslin + Fsnl + Fsb; - // System.out.println("calcMean for T = " + c.imt.toString()); - // System.out.printf(" Vs30 = %.2f\n", vs30); - // System.out.printf(" h = %.6f\n", h); - // System.out.printf(" Fp = %.6f\n", Fp); - // System.out.printf(" Fm = %.6f\n", Fm); - // System.out.printf(" Fd = %.6f\n", Fd); - // System.out.printf(" PGAref = %.6f\n", pgaRef); - // System.out.printf(" Fslin = %.6e\n", Fslin); - // System.out.printf(" Fsnl = %.6e\n", Fsnl); - // System.out.printf(" Fssb = %.6e\n", Fsb); - // System.out.printf(" Fs = %.6e\n", Fs); - // System.out.printf(" mu = %.6f\n", c.c0 + Fp + Fm + Fd + Fs); - // System.out.printf(" y = %.6f\n", exp(c.c0 + Fp + Fm + Fd + Fs)); - return c.c0 + Fp + Fm + Fd + Fs; } - /* Equations 4.2 and 4.3 */ - private static double getPathTerm(Coefficients c, double rRup, double Mw, double h) { - double R = Maths.hypot(rRup, h); - // PP The PEER report has just log(R) in middle term - return c.c1 * log(R) + c.b4 * Mw * log(rRup / R) + c.a0 * R; - } - - /* Equations 4.2 and 4.3, following official Python implementation */ - private static double getPathTermPy(Coefficients c, double rRup, double Mw, + /* Equation 2, 3 */ + private static double pathTerm(Coefficients c, double rRup, double Mw, double h) { double Rref = Maths.hypot(1, h); double R = Maths.hypot(rRup, h); double R_Rref = R / Rref; - // System.out.printf(" c1 = %.6e\n", c.c1); - // System.out.printf(" b4 = %.6e\n", c.b4); - // System.out.printf(" a0 = %.6e\n", c.a0); - // System.out.printf(" Rref = %.6e\n", R_Rref); - // System.out.printf(" R = %.6e\n", R); - // System.out.printf(" Mw = %.6e\n", Mw); - - return c.c1 * log(R) + c.b4 * Mw * log(R_Rref) + c.a0 * R; + return c.c1 * log(R) + B4 * Mw * log(R_Rref) + c.a0 * R; } private static double nearSourceSaturation(double Mw, double mc, boolean slab) { if (slab) { - /* Equation 4.4b */ + /* Equation 4b */ return (Mw <= mc) ? pow(10, (1.050 / (mc - 4.0)) * (Mw - mc) + 1.544) : 35; } else { - /* Equation 4.4a */ + /* Equation 4a */ return pow(10.0, -0.82 + 0.252 * Mw); } } - /* Equation 4.5 */ - private static double getMagnitudeTerm(Coefficients c, double Mw) { + /* Equation 5 */ + private static double magnitudeTerm(Coefficients c, double Mw) { double dM = Mw - c.mc; return (Mw <= c.mc) ? c.c4 * dM + c.c5 * dM * dM : c.c6 * dM; } - /* Equation 4.6 */ - private static double getDepthTerm(Coefficients c, boolean slab, double zHyp) { + /* Equation 6 */ + private static double depthTerm(Coefficients c, boolean slab, double zHyp) { if (!slab) { return 0.0; } - // zHyp (see section 4.3.3) can be replaced with mean depth (zHyp_bar) - // that depends on zTor, fault width, and dip angle + // zHyp (see section 4.3.3) could be replaced with mean depth + // (zHyp_bar) that depends on zTor, fault width, and dip angle double dTerm; - if (zHyp < DB1) { // "less than" + if (zHyp < DB1) { dTerm = c.m * (DB1 - DB2) + c.d; } else if (zHyp > DB2) { dTerm = c.d; - } else { // "db1 <= zHyp" possible typo? (written as "db1 < xHyp") + } else { dTerm = c.m * (zHyp - DB2) + c.d; } return dTerm; } - /* Equation 4.8 */ - private static double getSiteTermLinear(Coefficients c, double vs30) { + /* Equation 8 */ + private static double siteTermLinear(Coefficients c, double vs30) { double Flin; - if (vs30 <= c.v1) { - Flin = c.s1 * log(vs30 / c.v1) + c.s2 * log(c.v1 / c.vRef); + if (vs30 <= V1) { + Flin = c.s1 * log(vs30 / V1) + c.s2 * log(V1 / VREF); } else if (vs30 > c.v2) { - Flin = c.s2 * log(c.v2 / c.vRef); + Flin = c.s2 * log(c.v2 / VREF); } else { - Flin = c.s2 * log(vs30 / c.vRef); + Flin = c.s2 * log(vs30 / VREF); } return Flin; } - /* Equation 4.9 */ - private static double getSiteTermNonLinear(Coefficients c, double pgaRef, - double vs30) { + /* Equation 9 */ + private static double siteTermNonLinear(Coefficients c, double pgaRef, double vs30) { double f3 = 0.05; double f2 = c.f4 * (exp(c.f5 * (min(vs30, 760) - 200)) - exp(c.f5 * (760 - 200))); return 0.0 + f2 * log((pgaRef + f3) / f3); } - /* Equations 4.11/5.5 - 4.13/5.7 */ + /* Equation 11, 12, 13 */ private static double getSiteBasinTerm(Coefficients c, double vs30, double z2p5) { if (Double.isNaN(z2p5)) { @@ -435,172 +369,62 @@ public abstract class ParkerEtAl_2020 implements GroundMotionModel { double z2p5m = z2p5 * 1000.0; - // personal communication Parker (via Ahdi) use e2/e3 terms for general PNW - // basins - // Cascadia coefficients double θ0 = 3.94; double θ1 = -0.42; double νμ = 200.0; double νσ = 0.2; - double e1 = c.c_e1; - double e2 = c.c_e2;// + c.del_Seattle; - double e3 = c.c_e3;// + c.del_Seattle; - - // if (basin.equals(Basin.PACNW_NONE)) { - // e2 = c.c_e2 + c.del_None; - // e3 = c.c_e3 + c.del_None; - // } else if (basin.equals(Basin.SEATTLE)) { - // e2 = c.c_e2 + c.del_Seattle; - // e3 = c.c_e3 + c.del_Seattle; - // } else if (basin.equals(Basin.PACNW)) { - // e2 = c.c_e2; - // e3 = c.c_e3; - // } else { - // throw new IllegalArgumentException("Basin [" + basin + "] not - // supported"); - // } - - // Eq. 4.13/5.7 double lnμz2p5 = log(10) * (θ1 * (1 + Maths.erf((log10(vs30) - log10(νμ)) / νσ / sqrt(2))) + θ0); double δz2p5 = log(z2p5m) - lnμz2p5; double Fb; - if (δz2p5 <= (e1 / e3)) { - Fb = e1; - } else if (δz2p5 >= (e2 / e3)) { - Fb = e2; + if (δz2p5 <= (c.c_e1 / c.c_e3)) { + Fb = c.c_e1; + } else if (δz2p5 >= (c.c_e2 / c.c_e3)) { + Fb = c.c_e2; } else { - Fb = e3 * δz2p5; + Fb = c.c_e3 * δz2p5; } return Fb; } - /* Eqs. 6.3, 6.4, 6.5 */ - private static double getϕTotal(Coefficients c, double rRup, - double vs30) { + /* Eqsuation 18, 19, 20 */ + private static double phiTotal(Coefficients c, double rRup, double vs30) { - /* Eq. 6.5 */ double Δvar; - double rPrime = max(R1, min(R2, rRup)); - if (vs30 <= V1) { - Δvar = c.ϕ2v * (log(R2 / rPrime) / log(R2 / R1)); - } else if (vs30 >= V2) { + double rPrime = max(ΦR1, min(ΦR2, rRup)); + if (vs30 <= ΦV1) { + Δvar = c.ϕ2v * (log(ΦR2 / rPrime) / log(ΦR2 / ΦR1)); + } else if (vs30 >= ΦV2) { Δvar = 0; } else { - Δvar = c.ϕ2v * (log(V2 / vs30) / log(V2 / V1)) * (log(R2 / rPrime) / log(R2 / R1)); + Δvar = c.ϕ2v * + (log(ΦV2 / vs30) / log(ΦV2 / ΦV1)) * + (log(ΦR2 / rPrime) / log(ΦR2 / ΦR1)); } - /* Eq. 6.4 */ double ϕ2; - if (rRup <= R1) { + if (rRup <= ΦR1) { ϕ2 = c.ϕ21; - } else if (rRup >= R2) { + } else if (rRup >= ΦR2) { ϕ2 = c.ϕ22; } else { - ϕ2 = (c.ϕ22 - c.ϕ21) / log(R2 / R1) * log(rRup / R1) + c.ϕ21; + ϕ2 = (c.ϕ22 - c.ϕ21) / log(ΦR2 / ΦR1) * log(rRup / ΦR1) + c.ϕ21; } - /* Eq. 6.3 */ - double ϕTotal = sqrt(Δvar + ϕ2); - - // System.out.println("getϕTotal for T = " + c.imt.toString()); - // System.out.printf(" rRup = %.6f\n",rRup); - // System.out.printf(" vs30 = %.6f\n",vs30); - // System.out.printf(" rPrime = %.6f\n",rPrime); - // System.out.printf(" c.ϕ2v = %.6f\n",c.ϕ2v); - // System.out.printf(" c.ϕ21 = %.6f\n",c.ϕ21); - // System.out.printf(" c.ϕ22 = %.6f\n",c.ϕ22); - // System.out.printf(" Δvar = %.6f\n",Δvar); - // System.out.printf(" ϕ2 = %.6f\n",ϕ2); - // System.out.printf(" ϕTotal = %.6f\n",ϕTotal); - - return ϕTotal; + return sqrt(Δvar + ϕ2); } - /* Eqs. 6.6, 6.7 */ - private static double getϕS2S(Coefficients c, double rRup, double vs30) { - - /* Eq. 6.7 */ - double rPrime = max(R3, min(R4, rRup)); - double ΔvarS2S; - if (vs30 <= V3) { - ΔvarS2S = c.a1 * log(V3 / c.vm) * log(R4 / rPrime) / log(R4 / R3); - } else if (vs30 < c.vm) { - ΔvarS2S = c.a1 * log(vs30 / c.vm) * log(R4 / rPrime) / log(R4 / R3); - } else if (vs30 < V4) { - ΔvarS2S = c.a1 * log(vs30 / c.vm); - } else { // vs30 >= V4 - ΔvarS2S = c.a1 * log(V4 / c.vm); - } - - double ϕS2S = sqrt(c.ϕ2s2s0 + ΔvarS2S); - - // System.out.println("getϕS2S for T = " + c.imt.toString()); - // System.out.printf(" rRup = %.6f\n",rRup); - // System.out.printf(" vs30 = %.6f\n",vs30); - // System.out.printf(" rPrime = %.6f\n",rPrime); - // System.out.printf(" c.a1 = %.6f\n",c.a1); - // System.out.printf(" c.vm = %.6f\n",c.vm); - // System.out.printf(" ΔvarS2S = %.6f\n",ΔvarS2S); - // System.out.printf(" c.ϕ2s2s0 = %.6f\n",c.ϕ2s2s0); - // System.out.printf(" ϕS2S = %.6f\n",ϕS2S); - - return ϕS2S; - } - - /* Eqs. 6.8, 6.9, 6.10 */ - private double getϕSS(Coefficients c, double rRup, double vs30) { - /* Eq. 6.10 */ - double rPrime = max(R3, min(R4, rRup)); - double ΔvarSS; - if (vs30 <= V3) { - ΔvarSS = c.a2 * log(V3 / c.vm) * log(R4 / rPrime) / log(R4 / R3); - } else if (vs30 < c.vm) { - ΔvarSS = c.a2 * log(vs30 / c.vm) * log(R4 / rPrime) / log(R4 / R3); - } else if (vs30 < V4) { - ΔvarSS = c.a2 * log(vs30 / c.vm); - } else { // vs30 >= V4 - ΔvarSS = c.a2 * log(V4 / c.vm); - } - - /* Eq. 6.9 */ - double ϕSS2rRup; - if (rRup <= R5) { - ϕSS2rRup = c.ϕ2ss1; - } else if (rRup < R6) { - ϕSS2rRup = (c.ϕ2ss2 - c.ϕ2ss1) / log(R6 / R5) * log(rRup / R5) + c.ϕ2ss1; - } else { // rRup >= R6 - ϕSS2rRup = c.ϕ2ss2; - } - - double ϕSS = sqrt(ϕSS2rRup + ΔvarSS); - - // System.out.println("getϕSS for T = " + c.imt.toString()); - // System.out.printf(" rRup = %.6f\n",rRup); - // System.out.printf(" vs30 = %.6f\n",vs30); - // System.out.printf(" rPrime = %.6f\n",rPrime); - // System.out.printf(" c.a2 = %.6f\n",c.a1); - // System.out.printf(" c.vm = %.6f\n",c.vm); - // System.out.printf(" c.ϕ2ss1 = %.6f\n",c.ϕ2ss1); - // System.out.printf(" c.ϕ2ss2 = %.6f\n",c.ϕ2ss2); - // System.out.printf(" ΔvarSS = %.6f\n",ΔvarSS); - // System.out.printf(" ϕSS2rRup = %.6f\n",ϕSS2rRup); - // System.out.printf(" ϕS2S = %.6f\n",ϕSS); - - return ϕSS; - } - - /* Eq. 8.2 */ + /* Equation 27 */ private static double getσε(Coefficients c) { double σε; if (c.imt == Imt.PGA || c.imt == Imt.PGV || c.imt.period() < c.t1) { σε = c.σε1; } else if (c.imt.period() < c.t2) { σε = c.σε1 - (c.σε1 - c.σε2) * log(c.imt.period() / c.t1) / log(c.t2 / c.t1); - } else { // T > T2 + } else { σε = c.σε2; } return σε;