def enum_sigma_rho(cutoff, r_axis, ratio_alpha):
        """
        Find all possible sigma values and corresponding rotation angles
        within a sigma value cutoff with known rotation axis in rhombohedral system.
        The algorithm for this code is from reference, Acta Cryst, A45,505(1989).

        Args:
            cutoff (integer): the cutoff of sigma values.
            r_axis (list of three integers, e.g. u, v, w
                    or four integers, e.g. u, v, t, w):
                    the rotation axis of the grain boundary, with the format of [u,v,w]
                    or Weber indices [u, v, t, w].
            ratio_alpha (list of two integers, e.g. mu, mv):
                    mu/mv is the ratio of (1+2*cos(alpha))/cos(alpha) with rational number.
                    If irrational, set ratio_alpha = None.
        Returns:
            sigmas (dict):
                    dictionary with keys as the possible integer sigma values
                    and values as list of the possible rotation angles to the
                    corresponding sigma values.
                    e.g. the format as
                    {sigma1: [angle11,angle12,...], sigma2: [angle21, angle22,...],...}
                    Note: the angles are the rotation angle of one grain respect to the
                    other grain.
                    When generate the microstructure of the grain boundary using these
                    angles, you need to analyze the symmetry of the structure. Different
                    angles may result in equivalent microstructures.

        """
        sigmas = {}
        # transform four index notation to three index notation
        if len(r_axis) == 4:
            u1 = r_axis[0]
            v1 = r_axis[1]
            w1 = r_axis[3]
            u = 2 * u1 + v1 + w1
            v = v1 + w1 - u1
            w = w1 - 2 * v1 - u1
            r_axis = [u, v, w]
        # make sure gcd(r_axis)==1
        if reduce(gcd, r_axis) != 1:
            r_axis = [int(round(x / reduce(gcd, r_axis))) for x in r_axis]
        u, v, w = r_axis
        # make sure mu, mv are coprime integers.
        if ratio_alpha is None:
            mu, mv = [1, 1]
            if u + v + w != 0:
                if u != v or u != w:
                    raise RuntimeError('For irrational ratio_alpha, CSL only exist for [1,1,1]'
                                       'or [u, v, -(u+v)] and m =0')
        else:
            mu, mv = ratio_alpha
            if gcd(mu, mv) != 1:
                temp = gcd(mu, mv)
                mu = int(round(mu / temp))
                mv = int(round(mv / temp))

        # refer to the meaning of d in reference
        d = (u ** 2 + v ** 2 + w ** 2) * (mu - 2 * mv) + \
            2 * mv * (v * w + w * u + u * v)
        # Compute the max n we need to enumerate.
        n_max = int(np.sqrt((cutoff * abs(4 * mu * (mu - 3 * mv))) / abs(d)))

        # Enumerate all possible n, m to give possible sigmas within the cutoff.
        for n in range(1, n_max + 1):
            if ratio_alpha is None and u + v + w == 0:
                m_max = 0
            else:
                m_max = int(np.sqrt((cutoff * abs(4 * mu * (mu - 3 * mv)) - n ** 2 * d) / (mu)))
            for m in range(0, m_max + 1):
                if gcd(m, n) == 1 or m == 0:
                    # construct the rotation matrix, refer to the reference
                    R_list = [(mu - 2 * mv) * (u ** 2 - v ** 2 - w ** 2) * n ** 2 +
                              2 * mv * (v - w) * m * n - 2 * mv * v * w * n ** 2 +
                              mu * m ** 2,
                              2 * (mv * u * n * (w * n + u * n - m) - (mu - mv) *
                                   m * w * n + (mu - 2 * mv) * u * v * n ** 2),
                              2 * (mv * u * n * (v * n + u * n + m) + (mu - mv) *
                                   m * v * n + (mu - 2 * mv) * w * u * n ** 2),
                              2 * (mv * v * n * (w * n + v * n + m) + (mu - mv) *
                                   m * w * n + (mu - 2 * mv) * u * v * n ** 2),
                              (mu - 2 * mv) * (v ** 2 - w ** 2 - u ** 2) * n ** 2 +
                              2 * mv * (w - u) * m * n - 2 * mv * u * w * n ** 2 +
                              mu * m ** 2,
                              2 * (mv * v * n * (v * n + u * n - m) - (mu - mv) *
                                   m * u * n + (mu - 2 * mv) * w * v * n ** 2),
                              2 * (mv * w * n * (w * n + v * n - m) - (mu - mv) *
                                   m * v * n + (mu - 2 * mv) * w * u * n ** 2),
                              2 * (mv * w * n * (w * n + u * n + m) + (mu - mv) *
                                   m * u * n + (mu - 2 * mv) * w * v * n ** 2),
                              (mu - 2 * mv) * (w ** 2 - u ** 2 - v ** 2) * n ** 2 +
                              2 * mv * (u - v) * m * n - 2 * mv * u * v * n ** 2 +
                              mu * m ** 2]
                    m = -1 * m
                    # inverse of the rotation matrix
                    R_list_inv = [(mu - 2 * mv) * (u ** 2 - v ** 2 - w ** 2) * n ** 2 +
                                  2 * mv * (v - w) * m * n - 2 * mv * v * w * n ** 2 +
                                  mu * m ** 2,
                                  2 * (mv * u * n * (w * n + u * n - m) - (mu - mv) *
                                       m * w * n + (mu - 2 * mv) * u * v * n ** 2),
                                  2 * (mv * u * n * (v * n + u * n + m) + (mu - mv) *
                                       m * v * n + (mu - 2 * mv) * w * u * n ** 2),
                                  2 * (mv * v * n * (w * n + v * n + m) + (mu - mv) *
                                       m * w * n + (mu - 2 * mv) * u * v * n ** 2),
                                  (mu - 2 * mv) * (v ** 2 - w ** 2 - u ** 2) * n ** 2 +
                                  2 * mv * (w - u) * m * n - 2 * mv * u * w * n ** 2 +
                                  mu * m ** 2,
                                  2 * (mv * v * n * (v * n + u * n - m) - (mu - mv) *
                                       m * u * n + (mu - 2 * mv) * w * v * n ** 2),
                                  2 * (mv * w * n * (w * n + v * n - m) - (mu - mv) *
                                       m * v * n + (mu - 2 * mv) * w * u * n ** 2),
                                  2 * (mv * w * n * (w * n + u * n + m) + (mu - mv) *
                                       m * u * n + (mu - 2 * mv) * w * v * n ** 2),
                                  (mu - 2 * mv) * (w ** 2 - u ** 2 - v ** 2) * n ** 2 +
                                  2 * mv * (u - v) * m * n - 2 * mv * u * v * n ** 2 +
                                  mu * m ** 2]
                    m = -1 * m
                    F = mu * m ** 2 + d * n ** 2
                    all_list = R_list_inv + R_list + [F]
                    # Compute the max common factors for the elements of the rotation matrix
                    #  and its inverse.
                    com_fac = reduce(gcd, all_list)
                    sigma = int(round(abs(F / com_fac)))
                    if (sigma <= cutoff) and (sigma > 1):
                        if sigma not in list(sigmas.keys()):
                            if m == 0:
                                angle = 180.0
                            else:
                                angle = 2 * np.arctan(n / m * np.sqrt(d / mu)) \
                                        / np.pi * 180
                            sigmas[sigma] = [angle]
                        else:
                            if m == 0:
                                angle = 180
                            else:
                                angle = 2 * np.arctan(n / m * np.sqrt(d / mu)) \
                                        / np.pi * 180.0
                            if angle not in sigmas[sigma]:
                                sigmas[sigma].append(angle)
            if m_max == 0:
                break
        return sigmas