Find all possible sigma values and corresponding rotation angles
        within a sigma value cutoff with known rotation axis in cubic system.
        The algorithm for this code is from reference, Acta Cryst, A40,108(1984)
        Args:
            cutoff (integer): the cutoff of sigma values.
            r_axis (list of three integers, e.g. u, v, w):
                    the rotation axis of the grain boundary, with the format of [u,v,w].
        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 angles of one grain respect to
                    the other grain.
                    When generate the microstructures of the grain boundary using these angles,
                    you need to analyze the symmetry of the structure. Different angles may
                    result in equivalent microstructures.

def enum_sigma_cubic(cutoff, r_axis):
        """
        Find all possible sigma values and corresponding rotation angles
        within a sigma value cutoff with known rotation axis in cubic system.
        The algorithm for this code is from reference, Acta Cryst, A40,108(1984)
        Args:
            cutoff (integer): the cutoff of sigma values.
            r_axis (list of three integers, e.g. u, v, w):
                    the rotation axis of the grain boundary, with the format of [u,v,w].
        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 angles of one grain respect to
                    the other grain.
                    When generate the microstructures of the grain boundary using these angles,
                    you need to analyze the symmetry of the structure. Different angles may
                    result in equivalent microstructures.

        """
        sigmas = {}
        # 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]

        # count the number of odds in r_axis
        odd_r = len(list(filter(lambda x: x % 2 == 1, r_axis)))
        # Compute the max n we need to enumerate.
        if odd_r == 3:
            a_max = 4
        elif odd_r == 0:
            a_max = 1
        else:
            a_max = 2
        n_max = int(np.sqrt(cutoff * a_max / sum(np.array(r_axis) ** 2)))
        # enumerate all possible n, m to give possible sigmas within the cutoff.
        for n_loop in range(1, n_max + 1):
            n = n_loop
            m_max = int(np.sqrt(cutoff * a_max - n ** 2 * sum(np.array(r_axis) ** 2)))
            for m in range(0, m_max + 1):
                if gcd(m, n) == 1 or m == 0:
                    if m == 0:
                        n = 1
                    else:
                        n = n_loop
                    # construct the quadruple [m, U,V,W], count the number of odds in
                    # quadruple to determine the parameter a, refer to the reference
                    quadruple = [m] + [x * n for x in r_axis]
                    odd_qua = len(list(filter(lambda x: x % 2 == 1, quadruple)))
                    if odd_qua == 4:
                        a = 4
                    elif odd_qua == 2:
                        a = 2
                    else:
                        a = 1
                    sigma = int(round((m ** 2 + n ** 2 * sum(np.array(r_axis) ** 2)) / a))
                    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 * np.sqrt(sum(np.array(r_axis) ** 2)) / m) \
                                        / np.pi * 180
                            sigmas[sigma] = [angle]
                        else:
                            if m == 0:
                                angle = 180.0
                            else:
                                angle = 2 * np.arctan(n * np.sqrt(sum(np.array(r_axis) ** 2)) / m) \
                                        / np.pi * 180
                            if angle not in sigmas[sigma]:
                                sigmas[sigma].append(angle)
        return sigmas