PRVI legacy gmms

Updated and added GMMs for legacy 2003 PRVI NSHM

src/main/java/gov/usgs/earthquake/nshmp/gmm/MotazedianAtkinson_2003.java

+package gov.usgs.earthquake.nshmp.gmm;
+
+import static gov.usgs.earthquake.nshmp.gmm.GmmInput.Field.MW;
+import static gov.usgs.earthquake.nshmp.gmm.GmmInput.Field.RRUP;
+import static gov.usgs.earthquake.nshmp.gmm.GmmInput.Field.VS30;
+import static gov.usgs.earthquake.nshmp.gmm.GmmUtils.BASE_10_TO_E;
+import static gov.usgs.earthquake.nshmp.gmm.GmmUtils.LN_G_CM_TO_M;
+import static java.lang.Math.log10;
+import static java.lang.Math.sqrt;
+
+import java.util.Map;
+
+import com.google.common.collect.Range;
+
+import gov.usgs.earthquake.nshmp.gmm.GmmInput.Constraints;
+import gov.usgs.earthquake.nshmp.tree.LogicTree;
+
+/**
+ * Implementation of the preliminary MA05 model that was used for PRVI.
+ * Functional form is the same as MA05 with diffeent coefficients and constants.
+ * Coefficients were gleaned from forttran codes and only three IMTs are
+ * supported (PGA, 0.2 s, and 1 s).
+ *
+ * <p><b>Component:</b> not specified (avg horizontal implied)
+ *
+ * @author U.S. Geological Survey
+ * @see MotazedianAtkinson_2005
+ * @see Gmm#MA_03_PRVI
+ */
+public class MotazedianAtkinson_2003 implements GroundMotionModel {
+
+  static final String NAME = "Motazedian & Atkinson (2003) PRVI Preliminary";
+
+  static final Constraints CONSTRAINTS = Constraints.builder()
+      .set(MW, Range.closed(4.0, 9.5))
+      .set(RRUP, Range.closed(0.0, 1000.0))
+      .set(VS30, Range.closedOpen(150.0, 1500.0))
+      .build();
+
+  static final CoefficientContainer COEFFS =
+      new CoefficientContainer("MotazedianAtkinson03.csv");
+
+  // distance centering coefficients
+  private static final double C5 = -7.333;
+  private static final double C6 = 2.333;
+  // hinge function coefficients
+  private static final double C7 = -1.88;
+  private static final double C8 = 0.14;
+  private static final double DIST_HINGE_LOWER = 75.0;
+  private static final double DIST_HINGE_UPPER = 100.0;
+
+  // convert log10 sigma to ln
+  private static final double SIGMA = 0.28 * BASE_10_TO_E;
+
+  private static final class Coefficients {
+
+    final Imt imt;
+    final double c1, c2, c3, c4;
+
+    Coefficients(Imt imt, CoefficientContainer cc) {
+      this.imt = imt;
+      Map<String, Double> coeffs = cc.get(imt);
+      c1 = coeffs.get("c1");
+      c2 = coeffs.get("c2");
+      c3 = coeffs.get("c3");
+      c4 = coeffs.get("c4");
+    }
+  }
+
+  private final Coefficients coeffs;
+
+  MotazedianAtkinson_2003(Imt imt) {
+    coeffs = new Coefficients(imt, COEFFS);
+  }
+
+  @Override
+  public Imt imt() {
+    return coeffs.imt;
+  }
+
+  @Override
+  public LogicTree<GroundMotion> calc(GmmInput in) {
+    return calc(coeffs, in);
+  }
+
+  private static LogicTree<GroundMotion> calc(
+      Coefficients c,
+      GmmInput in) {
+
+    double μ = calcMean(c, in.rRup, in.Mw);
+    return GroundMotions.createTree(μ, SIGMA);
+  }
+
+  private static final double calcMean(
+      Coefficients c,
+      double rRup,
+      double Mw) {
+
+    /* Eq. 6, log10(PSA) in cm/s/s */
+
+    double r = getR(rRup, Mw);
+
+    double hingeFunction = getHingeFunction(r, Mw);
+
+    double μ = c.c1 + c.c2 * (Mw - 6) + c.c3 * (Mw - 6) * (Mw - 6) + hingeFunction + c.c4 * r;
+    // convert from base 10 to base e and from cm/s/s to g
+    return μ * BASE_10_TO_E - LN_G_CM_TO_M;
+  }
+
+  private static final double getR(double rRup, double Mw) {
+    // Equation and coefficients defined in Table 2 caption
+    double delta = C5 + C6 * Mw;
+    return sqrt(rRup * rRup + delta * delta);
+  }
+
+  private static final double getHingeFunction(double dist, double Mw) {
+    // Equation and coefficients defined in Table 2 caption
+    double mFac = C7 + C8 * Mw;
+    if (dist <= DIST_HINGE_LOWER) {
+      // R <= 75 km
+      return mFac * log10(dist);
+    } else if (dist <= DIST_HINGE_UPPER) {
+      // 75 km < R <= 100 km
+      return mFac * log10(DIST_HINGE_LOWER);
+    } else {
+      // R >= 100 km
+      return mFac * log10(DIST_HINGE_LOWER) - 0.5 * log10(dist / DIST_HINGE_UPPER);
+    }
+  }
+
+}