diff --git a/4_grains/gtGetReflections.m b/4_grains/gtGetReflections.m
new file mode 100644
index 0000000000000000000000000000000000000000..5171c16e2b37bae72e77ace3571b3412bc3dc4ed
--- /dev/null
+++ b/4_grains/gtGetReflections.m
@@ -0,0 +1,39 @@
+function hkl_reflections = gtGetReflections(hkl, spacegroup)
+% GTGETREFLECTIONS  Return matrix(n x 3) of all symmetry related reflections
+%                   from input hkl or hkil
+%
+%     hkl_reflections = gtGetReflections(hkl, spacegroup)
+%     ---------------------------------------------------------
+%
+%     INPUT:
+%       hkl        = <double> hkl or hkil vector, depending on crystal system
+%       spacegroup = <int>    Crystal spacegroup
+%
+%     OUTPUT:
+%       hkl_reflections = <double> hkl or hkil reflections, depending on 
+%                                  crystal system
+%     from input hkil [h k i l], return matrix(n x 4) of all symmetry related
+%     reflections.
+%     
+%     Note:
+%       Now, uses gtCrystGetSymmetryOperators.
+%       Before, was using gtGetHexagonalSymOp_sab / GetCubicSymOp modified
+%       from Risoe group
+
+% Get symmetry operators
+symm = gtCrystGetSymmetryOperators([], spacegroup);
+if size(hkl,2) == 4
+    symm = {symm.g};
+elseif size(hkl,2) == 3
+    symm = {symm.g3};
+end
+
+% Apply to input reflection
+for ii=1:length(symm)
+    hkl_reflections(ii,:) = hkl*symm{ii};
+end
+
+% Remove duplicates
+hkl_reflections = unique(hkl_reflections, 'rows');
+
+end % end of function