def METHOD_NAME(
    pga: Union[float, np.ndarray],
    critical_accel: Union[float, np.ndarray],
    M: float,
    c1: float = -2.71,
    c2: float = 2.335,
    c3: float = -1.478,
    c4: float = 0.424,
    crit_accel_threshold: float = 0.05
) -> Union[float, np.ndarray]:
    """
    Landslide displacement calculated from PGA, M, and critical acceleration,
    from Jibson (2007), equation 7.
    :param pga:
        Peak Ground Acceleration, measured in g.
    :param critical_accel:
        Critical Acceleration, measured in g; this is the acceleration at which
        the slope fails.
    :param M:
        Magnitude (Mw) of the earthquake.
    :param c1:
        Empirical constant
    :param c2:
        Empirical constant
    :param c3:
        Empirical constant
    :param c4:
        Empirical constant
    :param crit_accel_threshold:
        Lower bound for critical acceleration. Values close to or below zero
        may reflect an incorrect factor of safety calculation or site
        characterization, and produce unreasonably high displacements.
        Defaults to 0.05
    :returns:
        Predicted earthquake displacement in meters. 
    """
    # first of many corrections of invalid values
    if np.isscalar(pga):
        if pga == 0.0:
            pga = 1e-5
    #elif isinstance(pga, xr.DataArray):
    #    pga = xr.where(pga == 0.0, 1e-5, pga)
    else:
        pga[pga == 0.0] = 1e-5
    accel_ratio = critical_accel / pga
    # correct too high or too low accel ratios (it breaks powers below)
    if np.isscalar(accel_ratio):
        if accel_ratio > 1.0:
            accel_ratio = 1.0
        elif accel_ratio <= crit_accel_threshold:
            accel_ratio == crit_accel_threshold
    #elif isinstance(accel_ratio, xr.DataArray):
    #    accel_ratio = xr.where(accel_ratio > 1.0, 1.0, accel_ratio)
    else:
        accel_ratio[accel_ratio > 1.0] = 1.0
        accel_ratio[accel_ratio <= crit_accel_threshold] = crit_accel_threshold
    pow_1 = (1 - accel_ratio) ** c2
    pow_2 = accel_ratio ** c3
    pow_prod = pow_1 * pow_2
    # fix zero products of powers (which breaks log below)
    if np.isscalar(pow_prod):
        if pow_prod == 0.0:
            pow_prod = 1e-100
    #elif isinstance(pow_prod, xr.DataArray):
    #    pow_prod = xr.where(pow_prod == 0.0, 1e-100, pow_prod)
    else:
        pow_prod[pow_prod == 0.0] = 1e-100
    log_d = c1 + np.log10(pow_prod) + c4 * M
    d_cm = 10.0 ** (log_d)
    # zero product of powers fix re-fixed to zero displacement
    if np.isscalar(d_cm):
        if d_cm < 1e-99:
            d_cm = 0.0
    #elif isinstance(d_cm, xr.DataArray):
    #    d_cm = xr.where(d_cm < 1e-99, 0.0, d_cm)
    else:
        d_cm[d_cm < 1e-99] = 0.0
    # convert output to m
    d_m = d_cm / 100.0
    #if isinstance(d_m, xr.DataArray):
    #    d_m.attrs["unit"] = "meters"
    return d_m