def DAFrelative(alp, bet, k2_k1):
    """Relative Dynamic Amplification Factor, based on ratio of beam
    stiffness to some reference beam stiffness.

    Parameters
    ----------
    alp : float
        velocity ratio. alp=v/vcr; vcr= SQRT(2*SQRT(k*EI)/(rho*A))
    bet : float
        damping ratio.  bet=c/ccr; ccr=2*SQRT(rho*A*k)
    k2_k1 : float
        stiffness ratio of beam to current beam

    Returns
    -------
    DAF : float
        Relative Dynamic Amplification Factor

    Notes
    -----
    DAFrelative is a terse implementation for an appendix of code in
    Walker, R.T.R. and Indraratna, B, (in press) "Moving loads on a
    viscoelastic foundation with special reference to railway
    transition zones". International Journal of Geomechanics.

    """
    import numpy as np
    alp /= k2_k1 ** 0.25
    bet *= (1 / k2_k1) ** 0.5
    gam = np.roots([1.0, 0.0, 4.0 * alp ** 2, -8 * alp * bet, 4])
    cmat = np.zeros((4, 4), dtype=complex)
    cmat[(0), :] = [1, 1, -1, -1]
    cmat[(1), :] = [gam[0], gam[1], -gam[2], -gam[3]]
    cmat[(2), :] = [gam[0] ** 2, gam[1] ** 2, -gam[2] ** 2, -gam[3] ** 2]
    cmat[(3), :] = [gam[0] ** 3, gam[1] ** 3, -gam[2] ** 3, -gam[3] ** 3]
    rhs = np.array([0, 0, 0, 8.0], dtype=complex)
    A = np.linalg.solve(cmat, rhs)
    s = np.linspace(-8, 8, 1001, dtype=complex)
    disp = np.empty(len(s), dtype=complex)
    s_neg = s[s < 0]
    s_pos = s[s >= 0]
    disp[s >= 0] = A[0] * np.exp(s_pos * gam[0]) + A[1] * np.exp(s_pos * gam[1]
        )
    disp[s < 0] = A[2] * np.exp(s_neg * gam[2]) + A[3] * np.exp(s_neg * gam[3])
    DAF = np.real(np.max(disp) * (1 / k2_k1) ** 0.75)
    return DAF
