def rad_dis(theta, alpha, beta, order=2):
    """Radial distortion due to ellipticity or higher order distortion.

    Relative distortion, to be multiplied with radial distance.

    Parameters
    ----------
        theta : np.ndarray
            Angles at which to evaluate. Must be float.
        alpha : float
            Orientation of major axis.
        beta : float
            Strength of distortion (beta = (1-r_min/r_max)/(1+r_min/r_max).
        order : int
            Order of distortion.

    Returns
    -------
        : np.ndarray
            Distortion factor.
    """