def get_subcritical_coeffs(b, q0, hL, g):
    """2nd- & 0th-order term coefficients of subcritical Bernoulli relation.

    Args:
    -----
        b: a scalar or 1 1D numpy.ndarray; the topo elevation at target locations.
        q0: a scalar; the conservative quantity hu at x = 0 (the left) boundary.
        hL: a scalar; the depth h at x = L (the right) boundary.
        g: gravity.

    Returns:
    --------
        C0: the coefficient of the constant term of the polynomial.
        C2: the coefficient of the quadratic term of the polynomial.
    """

    C0 = q0 * q0 / (2. * g)
    C2 = b - C0 / (hL * hL) - hL

    return C0, C2