def circ_hz_pu(a, b):
    """Calculate the circularised horizontal PU from the semi-major and
    semi-minor axes of an error ellipse

    :param a: semi-major axis
    :param b: semi-minor axis
    :return: r the radius of the circularised error
    """
    q0 = 1.96079
    q1 = 0.004071
    q2 = 0.114276
    q3 = 0.371625
    c = b / a
    k = q0 + q1 * c + q2 * c ** 2 + q3 * c ** 3
    r = a * k
    return r
