diff --git a/conda.recipe/meta.yaml b/conda.recipe/meta.yaml
index 57f0945fb..cb5ea1b74 100644
--- a/conda.recipe/meta.yaml
+++ b/conda.recipe/meta.yaml
@@ -59,6 +59,9 @@ requirements:
     # Installing TBB means numba will use it instead of OpenMP.
     - tbb
     - tqdm
+    # needed for the phase transformation variants plotting.
+    # could also be used for pole density plotting in the future
+    - mplstereonet
 
 test:         # [build_platform == target_platform]
   imports:
diff --git a/hexrd/core/material/spacegroup.py b/hexrd/core/material/spacegroup.py
index 647f3777b..db261e4b5 100644
--- a/hexrd/core/material/spacegroup.py
+++ b/hexrd/core/material/spacegroup.py
@@ -87,7 +87,6 @@
 
 
 class SpaceGroup:
-
     def __init__(self, sgnum):
         """Constructor for SpaceGroup
 
@@ -373,6 +372,60 @@ def _sgrange(min, max):
 }
 
 
+def get_symmetry_directions(mat):
+    """
+    helper function to get a list of primary,
+    secondary and tertiary directions of the
+    space group of mat. e.g. cubic systems have
+        primary: [001]
+        secondary: [111]
+        tertiary: [110]
+
+    For some symmetries like monoclinic and trigonal,
+    some of the symmetry directions are not present. In
+    that case, we have choden them to be the same as the
+    crystal axis
+
+    For trigonal systems, it is ALWAYS assumed that they
+    are represented in the hexagonal basis. so the directions
+    are the same for the two
+    """
+    ltype = mat.latticeType
+
+    if ltype == "triclinic":
+        primary = np.array([0, 0, 1])
+        secondary = np.array([0, 1, 0])
+        tertiary = np.array([1, 0, 0])
+    elif ltype == "monoclinic":
+        primary = np.array([0, 1, 0])
+        secondary = np.array([1, 0, 0])
+        tertiary = np.array([0, 0, 1])
+    elif ltype == "orthorhombic":
+        primary = np.array([1, 0, 0])
+        secondary = np.array([0, 1, 0])
+        tertiary = np.array([1, 0, 1])
+    elif ltype == "tetragonal":
+        primary = np.array([0, 0, 1])
+        secondary = np.array([1, 0, 0])
+        tertiary = np.array([1, 1, 0])
+    elif ltype == "trigonal":
+        # it is assumed that trigonal crystals are
+        # represented in the hexagonal basis
+        primary = np.array([0, 0, 1])
+        secondary = np.array([1, 0, 0])
+        tertiary = np.array([1, 2, 0])
+    elif ltype == "hexagonal":
+        primary = np.array([0, 0, 1])
+        secondary = np.array([1, 0, 0])
+        tertiary = np.array([1, 2, 0])
+    elif ltype == "cubic":
+        primary = np.array([0, 0, 1])
+        secondary = np.array([1, 1, 1])
+        tertiary = np.array([1, 1, 0])
+
+    return (primary, secondary, tertiary)
+
+
 def Allowed_HKLs(sgnum, hkllist):
     """
     this function checks if a particular g vector is allowed
diff --git a/hexrd/core/valunits.py b/hexrd/core/valunits.py
index e0d00dbe7..d9e766910 100644
--- a/hexrd/core/valunits.py
+++ b/hexrd/core/valunits.py
@@ -313,6 +313,18 @@ def valWithDflt(val, dflt, toUnit=None):
     return retval
 
 
+def _nm(x):
+    return valWUnit("lp", "length", x, "nm")
+
+
+def _angstrom(x):
+    return valWUnit("lp", "length", x, "angstrom")
+
+
+def _kev(x):
+    return valWUnit("kev", "energy", x, "keV")
+
+
 if __name__ == '__main__':
     #
     #  doc testing
diff --git a/hexrd/singlextal/phasetransformation/variants.py b/hexrd/singlextal/phasetransformation/variants.py
new file mode 100644
index 000000000..da192d916
--- /dev/null
+++ b/hexrd/singlextal/phasetransformation/variants.py
@@ -0,0 +1,245 @@
+import numpy as np
+from hexrd.rotations import misorientation
+from hexrd.material import Material
+from matplotlib import pyplot as plt
+from hexrd.valunits import _angstrom, _kev
+import mplstereonet
+from numba import njit
+from hexrd.material.spacegroup import get_symmetry_directions
+
+"""
+define some helper functions
+"""
+
+"""
+    remove all hkl, -hkl pairs from list of symmetrically
+equivalent hkls"""
+
+
+def removeinversion(ksym):
+    klist = []
+    for i in range(ksym.shape[0]):
+        k = ksym[i, :]
+        kk = list(k)
+        nkk = list(-k)
+        if klist == []:
+            if np.sum(k) > np.sum(-k):
+                klist.append(kk)
+            else:
+                klist.append(nkk)
+
+        else:
+            if (kk in klist) or (nkk in klist):
+                pass
+            else:
+                klist.append(kk)
+    klist = np.array(klist)
+    return klist
+
+
+"""
+    get expected number of phase transformation variants
+    from group theoretic calculations
+"""
+
+
+def expected_num_variants(mat1, mat2, R1, R2):
+    """
+    calculate expected number of orientational
+    variants using the formula given in
+    C. Cayron, J. Appl. Cryst. (2007). 40, 1179–1182
+    page 2
+    N_alpha = |G_beta|/|H_beta|
+    """
+    T = np.dot(R1, R2.T)
+    sym1 = mat1.unitcell.SYM_PG_c
+    sym2 = mat2.unitcell.SYM_PG_c
+    ctr = 0
+    for s2 in sym2:
+        ss = np.dot(T, np.dot(s2, T.T))
+        for s1 in sym1:
+            diff = np.sum(np.abs(ss - s1))
+            if diff < 1e-6:
+                ctr += 1
+    return int(len(sym1) / ctr)
+
+
+"""
+    get rotation matrix
+"""
+
+
+def get_rmats(p1, d1, mat, toprint=False):
+    sym = mat.unitcell.SYM_PG_c
+    z = p1
+    x = d1
+    y = np.cross(z, x)
+    return np.vstack((x, y, z)).T
+
+
+"""
+    check if rotation matrix generated is a new one
+"""
+
+
+def isnew(mat, rmat, sym):
+    res = True
+    for r in rmat:
+        for s in sym:
+            rr = np.dot(s, r)
+            diff = np.sum(np.abs(rr - mat))
+            if diff < 1e-6:
+                res = False
+                break
+    return res
+
+
+def prepare_data(mat1, mat2, parallel_planes, parallel_directions):
+    """
+    prepare the planes and directions by:
+        1. converting to cartesian space
+        2. normalizing magnitude"""
+    p1, p2 = parallel_planes
+    d1, d2 = parallel_directions
+
+    p1 = mat1.unitcell.TransSpace(p1, "r", "c")
+    p1 = mat1.unitcell.NormVec(p1, "c")
+    d1 = mat1.unitcell.TransSpace(d1, "d", "c")
+    d1 = mat1.unitcell.NormVec(d1, "c")
+
+    p2 = mat2.unitcell.TransSpace(p2, "r", "c")
+    p2 = mat2.unitcell.NormVec(p2, "c")
+    d2 = mat2.unitcell.TransSpace(d2, "d", "c")
+    d2 = mat2.unitcell.NormVec(d2, "c")
+
+    return (p1, p2), (d1, d2)
+
+
+# main function
+# get the variants
+def getOR(R1, R2, mat1, mat2):
+    """
+    R1 ---> mat1
+    R2 ---> mat2
+    """
+    rmat = []
+    sym1 = mat1.unitcell.SYM_PG_c
+    sym2 = mat2.unitcell.SYM_PG_c
+    for sa in sym1:
+        ma = np.dot(sa, R1)
+        # for sb in sym2:
+        #     mb = np.dot(sb, R2)
+        mb = R2
+        m = np.dot(mb, ma.T)
+        if isnew(m, rmat, sym2):
+            rmat.append(m)
+
+    rmat_t = [r.T for r in rmat]
+    return rmat_t
+
+
+def plot_OR_mat(rmat, mat, fig, ax, title):
+    font = {
+        "family": "serif",
+        "weight": "bold",
+        "size": 12,
+    }
+
+    Gs = get_symmetry_directions(mat)
+    markers = ["ok", "^g", "sb"]
+
+    leg = []
+    for G in Gs:
+        leg.append(rf"$[{G[0]},{G[1]},{G[2]}]$")
+
+    for m, g in zip(markers, Gs):
+        hkl = mat.unitcell.CalcStar(g, "d")
+        dip = []
+        strike = []
+
+        for v in hkl:
+            vv = mat.unitcell.TransSpace(v, "d", "c")
+            vec = mat.unitcell.NormVec(vv, "c")
+            for r in rmat:
+                xyz = np.dot(r, vec)
+                if xyz[2] < 0.0:
+                    xyz = -xyz
+                dip.append(np.degrees(np.arccos(xyz[2])))
+                strike.append(180.0 + np.degrees(np.arctan2(xyz[1], xyz[0])))
+
+        ax.pole(strike, dip, m, markersize=12, markeredgecolor="red")
+    ax.tick_params(axis="both", size=24)
+    ax.grid(True, which="both", ls=":", lw=1)
+    ax.set_longitude_grid(10)
+    # ax.set_latitude_grid(10)
+    ax.set_azimuth_ticks([])
+    ax.legend(leg, bbox_to_anchor=(1.15, 1.15), loc="upper right", prop=font)
+    ax.set_title(title, **font)
+
+
+def plot_OR(rmat_parent, rmat_variants, mat1, mat2, fig=None, ax=None):
+    if fig is None or ax is None:
+        fig, ax = mplstereonet.subplots(
+            ncols=2, figsize=(10, 5), projection='equal_angle_stereonet'
+        )
+
+    plot_OR_mat(rmat_parent, mat1, fig, ax[0], title="Parent")
+    plot_OR_mat(rmat_variants, mat2, fig, ax[1], title="Variants")
+    # cosmetic
+    fig.subplots_adjust(
+        hspace=0, wspace=0.05, left=0.01, bottom=0.1, right=0.99
+    )
+    fig.show()
+
+
+def getPhaseTransformationVariants(
+    mat1,
+    mat2,
+    parallel_planes,
+    parallel_directions,
+    rmat_parent=[np.eye(3)],
+    plot=False,
+    verbose=False,
+):
+    """
+    main function to get the phase transformation variants between two materials
+    give a set of parallel planes and parallel directions
+
+    Parameters
+    ----------
+    mat1 : hexrd.material.Material
+        Material class for parent phase
+    mat2 : hexrd.material.Material
+        Material class for child phase
+    parallel_planes : list/tuple
+        list of numpy arrays with parallel planes, length=2
+    parallel_planes : list/tuple
+        list of numpy arrays with parallel directions, length=2
+    plot : boolean
+        plot the data in stereographic projection if true
+
+
+    """
+    (p1, p2), (d1, d2) = prepare_data(
+        mat1, mat2, parallel_planes, parallel_directions
+    )
+
+    R1 = get_rmats(p1, d1, mat1)
+    R2 = get_rmats(p2, d2, mat2)
+
+    rmat_variants = getOR(R1, R2, mat1, mat2)
+
+    num_var = expected_num_variants(mat1, mat2, R1, R2)
+
+    if verbose:
+        print("Expected # of orientational variants = ", num_var)
+        print(
+            "number of orientation variants in phase transition = ",
+            len(rmat_variants),
+            "\n",
+        )
+
+    if plot:
+        plot_OR(rmat_parent, rmat_variants, mat1, mat2, fig=None, ax=None)
+
+    return rmat_variants, num_var
diff --git a/setup.py b/setup.py
index dd567034c..844a7f9eb 100644
--- a/setup.py
+++ b/setup.py
@@ -31,6 +31,7 @@
     'scipy',
     'tqdm',
     'xxhash',
+    'mplstereonet',
 ]
 
 if platform.machine() == 'x86_64':
diff --git a/tests/data/materials/materials_variants.h5 b/tests/data/materials/materials_variants.h5
new file mode 100644
index 000000000..41578477c
Binary files /dev/null and b/tests/data/materials/materials_variants.h5 differ
diff --git a/tests/test_variants.py b/tests/test_variants.py
new file mode 100644
index 000000000..22002cb48
--- /dev/null
+++ b/tests/test_variants.py
@@ -0,0 +1,45 @@
+from hexrd.core.material import Material, load_materials_hdf5
+import h5py
+import numpy as np
+from hexrd.singlextal.phasetransformation import variants
+from hexrd.valunits import _kev, _nm
+import pytest
+
+
+@pytest.fixture
+def test_variants_material_file(test_data_dir):
+    return f'{test_data_dir}/materials/materials_variants.h5'
+
+
+def test_variants(test_variants_material_file):
+    beamenergy = _kev(10)
+    dmin = _nm(0.075)
+    print(test_variants_material_file)
+    mats = load_materials_hdf5(
+        test_variants_material_file, dmin=dmin, kev=beamenergy
+    )
+
+    # 4H-SiC phase crystal structure
+    fourH = mats["4H_SiC_HP"]
+
+    # B1 phase crystal structure
+    B1 = mats["SiC_B1"]
+
+    p1 = np.array([0.0, 0.0, 1.0])
+    p2 = np.array([0.0, 0.0, 1.0])
+    d1 = np.array([2.0, 1.0, 0.0])
+    d2 = np.array([1.0, 1.0, 0.0])
+    parallel_planes = (p1, p2)
+    parallel_directions = (d1, d2)
+
+    rmat_variants, num_variants = variants.getPhaseTransformationVariants(
+        fourH,
+        B1,
+        parallel_planes,
+        parallel_directions,
+        plot=False,
+        verbose=False,
+    )
+
+    assert len(rmat_variants) == num_variants
+    assert len(rmat_variants) == 3