Parker et al. cleanup

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 σε;