Find all possible sigma values and corresponding rotation angles
        within a sigma value cutoff with known rotation axis in tetragonal system.
        The algorithm for this code is from reference, Acta Cryst, B46,117(1990)

        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].
            c2_a2_ratio (list of two integers, e.g. mu, mv):
                    mu/mv is the square of the tetragonal axial ratio with rational number.
                    if irrational, set c2_a2_ratio = 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.

def enum_sigma_tet(cutoff, r_axis, c2_a2_ratio):
        """
        Find all possible sigma values and corresponding rotation angles
        within a sigma value cutoff with known rotation axis in tetragonal system.
        The algorithm for this code is from reference, Acta Cryst, B46,117(1990)

        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].
            c2_a2_ratio (list of two integers, e.g. mu, mv):
                    mu/mv is the square of the tetragonal axial ratio with rational number.
                    if irrational, set c2_a2_ratio = 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 = {}
        # 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 c2_a2_ratio is None:
            mu, mv = [1, 1]
            if w != 0:
                if u != 0 or (v != 0):
                    raise RuntimeError('For irrational c2/a2, CSL only exist for [0,0,1] '
                                       'or [u,v,0] and m = 0')
        else:
            mu, mv = c2_a2_ratio
            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) * mv + w ** 2 * mu

        # Compute the max n we need to enumerate.
        n_max = int(np.sqrt((cutoff * 4 * mu * mv) / d))

        # Enumerate all possible n, m to give possible sigmas within the cutoff.
        for n in range(1, n_max + 1):
            if c2_a2_ratio is None and w == 0:
                m_max = 0
            else:
                m_max = int(np.sqrt((cutoff * 4 * mu * 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 = [(u ** 2 * mv - v ** 2 * mv - w ** 2 * mu) * n ** 2 +
                              mu * m ** 2,
                              2 * v * u * mv * n ** 2 - 2 * w * mu * m * n,
                              2 * u * w * mu * n ** 2 + 2 * v * mu * m * n,
                              2 * u * v * mv * n ** 2 + 2 * w * mu * m * n,
                              (v ** 2 * mv - u ** 2 * mv - w ** 2 * mu) * n ** 2 +
                              mu * m ** 2,
                              2 * v * w * mu * n ** 2 - 2 * u * mu * m * n,
                              2 * u * w * mv * n ** 2 - 2 * v * mv * m * n,
                              2 * v * w * mv * n ** 2 + 2 * u * mv * m * n,
                              (w ** 2 * mu - u ** 2 * mv - v ** 2 * mv) * n ** 2 +
                              mu * m ** 2]
                    m = -1 * m
                    # inverse of rotation matrix
                    R_list_inv = [(u ** 2 * mv - v ** 2 * mv - w ** 2 * mu) * n ** 2 +
                                  mu * m ** 2,
                                  2 * v * u * mv * n ** 2 - 2 * w * mu * m * n,
                                  2 * u * w * mu * n ** 2 + 2 * v * mu * m * n,
                                  2 * u * v * mv * n ** 2 + 2 * w * mu * m * n,
                                  (v ** 2 * mv - u ** 2 * mv - w ** 2 * mu) * n ** 2 +
                                  mu * m ** 2,
                                  2 * v * w * mu * n ** 2 - 2 * u * mu * m * n,
                                  2 * u * w * mv * n ** 2 - 2 * v * mv * m * n,
                                  2 * v * w * mv * n ** 2 + 2 * u * mv * m * n,
                                  (w ** 2 * mu - u ** 2 * mv - v ** 2 * mv) * n ** 2 +
                                  mu * m ** 2]
                    m = -1 * m
                    F = mu * m ** 2 + d * n ** 2
                    all_list = R_list + R_list_inv + [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((mu * m ** 2 + d * n ** 2) / 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.0
                            else:
                                angle = 2 * np.arctan(n / m * np.sqrt(d / mu)) \
                                        / np.pi * 180
                            if angle not in sigmas[sigma]:
                                sigmas[sigma].append(angle)
            if m_max == 0:
                break

        return sigmas