Compute the shear modulus reduction using Equation (1).
503280
def _calc_mod_reduc(self, strains, strain_ref, x_1, x_1_mean, x_2,
                        x_2_mean, x_3, x_3_mean):
        """"""

        ones = np.ones_like(strains)
        # Predictor
        x_4 = np.log(self._lab_consol_ratio) * ones
        x = np.c_[ones, x_1, x_2, x_3, x_4, (x_1 - x_1_mean) * (
            x_2 - x_2_mean), (x_1 - x_1_mean) * (x_3 - x_3_mean), (
                x_2 - x_2_mean) * (x_3 - x_3_mean), (x_1 - x_1_mean) * (
                    x_2 - x_2_mean) * (x_3 - x_3_mean)]
        # Coefficients
        denom = np.log(1 / strain_ref + strains / strain_ref)  # TODO: is this percent or decimal?
        b = np.c_[5.11 * ones, -0.729 * ones, (1 - 0.37 * x_3_mean * (1 + ((
            np.log(strain_ref) - x_1_mean) / denom))), -0.693 * ones, 0.8 - 0.4
                  * x_3, 0.37 * x_3_mean / denom, 0.0 * ones, -0.37 * (1 + (
                      np.log(strain_ref) - x_1_mean) / denom), 0.37 / denom, ]
        ln_shear_mod = (b * x).sum(axis=1)
        shear_mod = np.exp(ln_shear_mod)
        mod_reduc = shear_mod / shear_mod[0]
        return mod_reduc
628535
function () {
		var parts = this._parts,
		    tolerance = this.options.smoothFactor;

		for (var i = 0, len = parts.length; i < len; i++) {
			parts[i] = L.LineUtil.simplify(parts[i], tolerance);
		}
	}
941112
def FD(self):
    
    if self.Verbose:
      print("Finite Difference Solution Technique")
    # Used to check for coeff_matrix here, but now doing so in self.bc_check()
    # called by f1d and f2d at the start
    # 
    # Define a stress-based qs = q0
    # But only if the latter has not already been defined
    # (e.g., by the getters and setters)
    try:
      self.qs
    except:
      self.qs = self.q0.copy()
      # Remove self.q0 to avoid issues with multiply-defined inputs
      # q0 is the parsable input to either a qs grid or contains (x,(y),q)
      del self.q0
    # Give it x and y dimensions for help with plotting tools
    # (not implemented internally, but a help with external methods)
    self.x = np.arange(self.dx/2., self.dx * self.qs.shape[0], self.dx)
    if self.dimension == 2:
      self.y = np.arange(self.dy/2., self.dy * self.qs.shape[1], self.dy)
    # Is there a solver defined
    try:
      self.Solver # See if it exists already
    except:
      # Well, will fail if it doesn't see this, maybe not the most reasonable
      # error message.
      if self.filename:
        self.Solver = self.configGet("string", "numerical", "Solver")
      else:
        sys.exit("No solver defined!")
    # Check consistency of size if coeff array was loaded
    if self.filename:
      # In the case that it is iterative, find the convergence criterion
      self.iterative_ConvergenceTolerance = self.configGet("float", "numerical", "ConvergenceTolerance")    
      # Try to import Te grid or scalar for the finite difference solution
      try:
        self.Te = self.configGet("float", "input", "ElasticThickness", optional=False)
        if self.Te is None:
          Tepath = self.configGet("string", "input", "ElasticThickness", optional=False)
          self.Te = Tepath
        else:
          Tepath = None
      except:
        Tepath = self.configGet("string", "input", "ElasticThickness", optional=False)
        self.Te = Tepath
      if self.Te is None:
        if self.coeff_matrix is not None:
          pass
        else:
          # Have to bring this out here in case it was discovered in the 
          # try statement that there is no value given
          sys.exit("No input elastic thickness or coefficient matrix supplied.")
    # or if getter/setter
    if type(self.Te) == str: 
      # Try to import Te grid or scalar for the finite difference solution
      Tepath = self.Te
    else:
      Tepath = None # in case no self.filename present (like for GRASS GIS)
    # If there is a Tepath, import Te
    # Assume that even if a coeff_matrix is defined
    # That the user wants Te if they gave the path
    if Tepath:
      self.Te = self.loadFile(self.Te, close_on_fail = False)
      if self.Te is None:
        print("Requested Te file is provided but cannot be located.")
        print("No scalar elastic thickness is provided in configuration file")
        print("(Typo in path to input Te grid?)")
        if self.coeff_matrix is not None:
          print("But a coefficient matrix has been found.")
          print("Calculations will be carried forward using it.")
        else:
          print("Exiting.")
          sys.exit()

      # Check that Te is the proper size if it was loaded
      # Will be array if it was loaded
      if self.Te.any():
        self.TeArraySizeCheck()
1202362
function cubicBezierAtTime(t,p1x,p1y,p2x,p2y,duration) {
      var ax = 0,bx = 0,cx = 0,ay = 0,by = 0,cy = 0;
      function sampleCurveX(t) {return ((ax * t + bx) * t + cx) * t;}
      function sampleCurveY(t) {return ((ay * t + by) * t + cy) * t;}
      function sampleCurveDerivativeX(t) {return (3.0 * ax * t + 2.0 * bx) * t + cx;}
      function solveEpsilon(duration) {return 1.0 / (200.0 * duration);}
      function solve(x,epsilon) {return sampleCurveY(solveCurveX(x,epsilon));}
      function fabs(n) {if (n >= 0) {return n;}else {return 0 - n;}}
      function solveCurveX(x,epsilon) {
        var t0,t1,t2,x2,d2,i;
        for (t2 = x, i = 0; i < 8; i++) {x2 = sampleCurveX(t2) - x; if (fabs(x2) < epsilon) {return t2;} d2 = sampleCurveDerivativeX(t2); if (fabs(d2) < 1e-6) {break;} t2 = t2 - x2 / d2;}
        t0 = 0.0; t1 = 1.0; t2 = x; if (t2 < t0) {return t0;} if (t2 > t1) {return t1;}
        while (t0 < t1) {x2 = sampleCurveX(t2); if (fabs(x2 - x) < epsilon) {return t2;} if (x > x2) {t0 = t2;}else {t1 = t2;} t2 = (t1 - t0) * 0.5 + t0;}
        return t2; // Failure.
      }
      cx = 3.0 * p1x; bx = 3.0 * (p2x - p1x) - cx; ax = 1.0 - cx - bx; cy = 3.0 * p1y; by = 3.0 * (p2y - p1y) - cy; ay = 1.0 - cy - by;
      return solve(t, solveEpsilon(duration));
    }
1222093
function calculateBezier(t, p0, p1, p2, p3) {
  var mt = 1 - t;
  return mt * mt * mt * p0 + 3 * mt * mt * t * p1 + 3 * mt * t * t * p2 + t * t * t * p3;
}
754983
def rhochange(self):
        """"""

        if self.opt['HighMemSolve'] and self.cri.Cd == 1:
            self.c = sl.solvedbd_sm_c(
                self.Df, np.conj(self.Df),
                (self.mu / self.rho) * self.GHGf + 1.0, self.cri.axisM)
457325
def ccmod_xstep(k):
    

    YU = mp_D_Y - mp_D_U[k]
    b = mp_ZSf[k] + mp_drho * sl.rfftn(YU, None, mp_cri.axisN)
    Xf = sl.solvedbi_sm(mp_Zf[k], mp_drho, b, axis=mp_cri.axisM)
    mp_D_X[k] = sl.irfftn(Xf, mp_cri.Nv, mp_cri.axisN)
367460
function spring(string, duration) {

  var params = parseEasingParameters(string);
  var mass = minMax(is.und(params[0]) ? 1 : params[0], .1, 100);
  var stiffness = minMax(is.und(params[1]) ? 100 : params[1], .1, 100);
  var damping = minMax(is.und(params[2]) ? 10 : params[2], .1, 100);
  var velocity =  minMax(is.und(params[3]) ? 0 : params[3], .1, 100);
  var w0 = Math.sqrt(stiffness / mass);
  var zeta = damping / (2 * Math.sqrt(stiffness * mass));
  var wd = zeta < 1 ? w0 * Math.sqrt(1 - zeta * zeta) : 0;
  var a = 1;
  var b = zeta < 1 ? (zeta * w0 + -velocity) / wd : -velocity + w0;

  function solver(t) {
    var progress = duration ? (duration * t) / 1000 : t;
    if (zeta < 1) {
      progress = Math.exp(-progress * zeta * w0) * (a * Math.cos(wd * progress) + b * Math.sin(wd * progress));
    } else {
      progress = (a + b * progress) * Math.exp(-progress * w0);
    }
    if (t === 0 || t === 1) { return t; }
    return 1 - progress;
  }

  function getDuration() {
    var cached = cache.springs[string];
    if (cached) { return cached; }
    var frame = 1/6;
    var elapsed = 0;
    var rest = 0;
    while(true) {
      elapsed += frame;
      if (solver(elapsed) === 1) {
        rest++;
        if (rest >= 16) { break; }
      } else {
        rest = 0;
      }
    }
    var duration = elapsed * frame * 1000;
    cache.springs[string] = duration;
    return duration;
  }

  return duration ? solver : getDuration;

}
498498
def get_shear_vel(self, saturated):
        
        try:
            if saturated:
                return np.sqrt(self.g_mod / self.unit_sat_mass)
            else:
                return np.sqrt(self.g_mod / self.unit_dry_mass)
        except TypeError:
            return None
888741
def newton_update(x, a, c, step=1.0):
    
    return x - step*np.linalg.inv(hessian(x, a)).dot(gradient(x, a, c))
1369181
def _get_distance_term(self, C, rjb, mag):
        
        c_3 = self._get_anelastic_coeff(C)
        rval = np.sqrt(rjb ** 2. + C["h"] ** 2.)
        return (C["c1"] + C["c2"] * (mag - self.CONSTS["Mref"])) *\
            np.log(rval / self.CONSTS["Rref"]) +\
            c_3 * (rval - self.CONSTS["Rref"])
641132
def hessian(self, x, y, theta_E, r_trunc, center_x=0, center_y=0):
        
        x_shift = x - center_x
        y_shift = y - center_y
        dphi_dr = self._dphi_dr(x_shift, y_shift, theta_E, r_trunc)
        d2phi_dr2 = self._d2phi_dr2(x_shift, y_shift, theta_E, r_trunc)
        dr_dx, dr_dy = self._dr_dx(x, y)
        d2r_dx2, d2r_dy2, d2r_dxy = self._d2r_dx2(x_shift, y_shift)
        f_xx = d2r_dx2*dphi_dr + dr_dx**2*d2phi_dr2
        f_yy = d2r_dy2*dphi_dr + dr_dy**2*d2phi_dr2
        f_xy = d2r_dxy*dphi_dr + dr_dx*dr_dy*d2phi_dr2
        return f_xx, f_yy, f_xy
362568
def polynomial_reduce_mod( poly, polymod, p ):
  """"""

  # This module has been tested only by extensive use
  # in calculating modular square roots.

  # Just to make this easy, require a monic polynomial:
  assert polymod[-1] == 1

  assert len( polymod ) > 1

  while len( poly ) >= len( polymod ):
    if poly[-1] != 0:
      for i in range( 2, len( polymod ) + 1 ):
        poly[-i] = ( poly[-i] - poly[-1] * polymod[-i] ) % p
    poly = poly[0:-1]

  return poly
234502
def _make_complex(self):
        """"""
        rcomplex_coeffs = _shtools.SHrtoc(self.coeffs,
                                          convention=1, switchcs=0)

        # These coefficients are using real floats, and need to be
        # converted to complex form.
        complex_coeffs = _np.zeros((2, self.lmax+1, self.lmax+1),
                                   dtype='complex')
        complex_coeffs[0, :, :] = (rcomplex_coeffs[0, :, :] + 1j *
                                   rcomplex_coeffs[1, :, :])
        complex_coeffs[1, :, :] = complex_coeffs[0, :, :].conjugate()
        for m in self.degrees():
            if m % 2 == 1:
                complex_coeffs[1, :, m] = - complex_coeffs[1, :, m]

        # complex_coeffs is initialized in this function and can be
        # passed as reference
        return SHCoeffs.from_array(complex_coeffs,
                                   normalization=self.normalization,
                                   csphase=self.csphase, copy=False)
1240381
def q12d(grid, spacing=None, E=1e5, nu=0.3, dirichlet_boundary=True,
         format=None):
    
    X, Y = tuple(grid)

    if X < 1 or Y < 1:
        raise ValueError('invalid grid shape')

    if dirichlet_boundary:
        X += 1
        Y += 1

    pts = np.mgrid[0:X+1, 0:Y+1]
    pts = np.hstack((pts[0].T.reshape(-1, 1) - X / 2.0,
                     pts[1].T.reshape(-1, 1) - Y / 2.0))

    if spacing is None:
        DX, DY = 1, 1
    else:
        DX, DY = tuple(spacing)
        pts *= [DX, DY]

    # compute local stiffness matrix
    lame = E * nu / ((1 + nu) * (1 - 2*nu))  # Lame's first parameter
    mu = E / (2 + 2*nu)                   # shear modulus

    vertices = np.array([[0, 0], [DX, 0], [DX, DY], [0, DY]])
    K = q12d_local(vertices, lame, mu)

    nodes = np.arange((X+1)*(Y+1)).reshape(X+1, Y+1)
    LL = nodes[:-1, :-1]
    Id = (2*LL).repeat(K.size).reshape(-1, 8, 8)
    J = Id.copy()
    Id += np.tile([0, 1, 2, 3, 2*X + 4, 2*X + 5, 2*X + 2, 2*X + 3], (8, 1))
    J += np.tile([0, 1, 2, 3, 2*X + 4, 2*X + 5, 2*X + 2, 2*X + 3], (8, 1)).T
    V = np.tile(K, (X*Y, 1))

    Id = np.ravel(Id)
    J = np.ravel(J)
    V = np.ravel(V)

    # sum duplicates
    A = coo_matrix((V, (Id, J)), shape=(pts.size, pts.size)).tocsr()
    A = A.tobsr(blocksize=(2, 2))

    del Id, J, V, LL, nodes

    B = np.zeros((2 * (X+1)*(Y+1), 3))
    B[0::2, 0] = 1
    B[1::2, 1] = 1
    B[0::2, 2] = -pts[:, 1]
    B[1::2, 2] = pts[:, 0]

    if dirichlet_boundary:
        mask = np.zeros((X+1, Y+1), dtype='bool')
        mask[1:-1, 1:-1] = True
        mask = np.ravel(mask)
        data = np.zeros(((X-1)*(Y-1), 2, 2))
        data[:, 0, 0] = 1
        data[:, 1, 1] = 1
        indices = np.arange((X-1)*(Y-1))
        indptr = np.concatenate((np.array([0]), np.cumsum(mask)))
        P = bsr_matrix((data, indices, indptr),
                       shape=(2*(X+1)*(Y+1), 2*(X-1)*(Y-1)))
        Pt = P.T
        A = P.T * A * P

        B = Pt * B

    return A.asformat(format), B
254440
def _get_anelastic_attenuation_term(self, C, rrup):
        
        f_atn = np.zeros(len(rrup))
        idx = rrup >= 80.0
        f_atn[idx] = (C["c20"] + C["Dc20"]) * (rrup[idx] - 80.0)
        return f_atn
941615
def apply(self, f_del, h):
        

        fd_rule = self.rule

        ne = h.shape[0]
        nr = fd_rule.size - 1
        _assert(nr < ne, 'num_steps ({0:d}) must  be larger than '
                '({1:d}) n + order - 1 = {2:d} + {3:d} -1'
                ' ({4:s})'.format(ne, nr+1, self.n, self.order, self.method))
        f_diff = convolve(f_del, fd_rule[::-1], axis=0, origin=nr // 2)

        der_init = f_diff / (h ** self.n)
        ne = max(ne - nr, 1)
        return der_init[:ne], h[:ne]