Add feature for multiple lattice parameters/phases in polycrystal_from_odf

docs/source/examples/example_polycrystal_from_odf_2phases.py

@@ -0,0 +1,31 @@
+import numpy as np
+from xrd_simulator.templates import polycrystal_from_odf
+
+
+# uniform orientation distribution function.
+def ODF(x, q): return 1. / (np.pi**2)
+
+
+number_of_crystals = 500
+bounding_height = 50.0
+bounding_radius = 25.0
+unit_cell = [[3.579, 3.579, 3.579, 90., 90., 90.0],
+            [5.46745, 5.46745, 5.46745, 90., 90., 90.0]]
+sgname = ['F432', 'F432']
+max_bin = np.radians(10.0)
+path_to_cif_file = None
+
+def strain_tensor(x): return np.array([[0, 0, 0.02 * x[2] / bounding_height],
+                                       [0, 0, 0],
+                                       [0, 0, 0]])  # Linear strain gradient along rotation axis.
+
+
+polycrystal = polycrystal_from_odf(ODF,
+                                   number_of_crystals,
+                                   bounding_height,
+                                   bounding_radius,
+                                   unit_cell,
+                                   sgname,
+                                   path_to_cif_file,
+                                   max_bin,
+                                   strain_tensor)

tests/test_templates.py

@@ -175,6 +175,116 @@ def strain_tensor(x): return np.array(
             nosequences,
             10,
             msg="Few or no rings appeared from diffraction.")
+
+    def test_polycrystal_from_odf_2phases(self):
+
+        unit_cell = [[3.579, 3.579, 3.579, 90., 90., 90.0],
+                     [5.46745, 5.46745, 5.46745, 90., 90., 90.0]]
+        sgname = ['F432', 'F432']
+
+        def orientation_density_function(
+            x, q): return 1. / (np.pi**2)  # uniform ODF
+        number_of_crystals = 500
+        sample_bounding_cylinder_height = 50
+        sample_bounding_cylinder_radius = 25
+        maximum_sampling_bin_seperation = np.radians(10.0)
+        # Linear strain gradient along rotation axis.
+        def strain_tensor(x): return np.array(
+            [[0, 0, 0.02 * x[2] / sample_bounding_cylinder_height], [0, 0, 0], [0, 0, 0]])
+
+        polycrystal = templates.polycrystal_from_odf(
+            orientation_density_function,
+            number_of_crystals,
+            sample_bounding_cylinder_height,
+            sample_bounding_cylinder_radius,
+            unit_cell,
+            sgname,
+            path_to_cif_file=None,
+            maximum_sampling_bin_seperation=maximum_sampling_bin_seperation,
+            strain_tensor=strain_tensor)
+
+        # Compare Euler angle distributions to scipy random uniform orientation
+        # sampler
+        euler1 = np.array([Rotation.from_matrix(U).as_euler(
+            'xyz', degrees=True) for U in polycrystal.orientation_lab])
+        euler2 = Rotation.random(10 * euler1.shape[0]).as_euler('xyz')
+
+        for i in range(3):
+            hist1, bins = np.histogram(euler1[:, i])
+            hist2, bins = np.histogram(euler2[:, i])
+            hist2 = hist2 / 10.
+            # These histograms should look roughly the same
+            self.assertLessEqual(
+                np.max(
+                    np.abs(
+                        hist1 -
+                        hist2)),
+                number_of_crystals *
+                0.05,
+                "ODF not sampled correctly.")
+
+        parameters = {
+            "detector_distance": 191023.9164,
+            "detector_center_pixel_z": 256.2345,
+            "detector_center_pixel_y": 255.1129,
+            "pixel_side_length_z": 181.4234,
+            "pixel_side_length_y": 180.2343,
+            "number_of_detector_pixels_z": 512,
+            "number_of_detector_pixels_y": 512,
+            "wavelength": 0.285227,
+            "beam_side_length_z": 512 * 200.,
+            "beam_side_length_y": 512 * 200.,
+            "rotation_step": np.radians(20.0),
+            "rotation_axis": np.array([0., 0., 1.0])
+        }
+
+        beam, detector, motion = templates.s3dxrd(parameters)
+
+        number_of_crystals = 100
+        sample_bounding_cylinder_height = 256 * 180 / 128.
+        sample_bounding_cylinder_radius = 256 * 180 / 128.
+
+        polycrystal = templates.polycrystal_from_odf(
+            orientation_density_function,
+            number_of_crystals,
+            sample_bounding_cylinder_height,
+            sample_bounding_cylinder_radius,
+            unit_cell,
+            sgname,
+            path_to_cif_file=None,
+            maximum_sampling_bin_seperation=maximum_sampling_bin_seperation,
+            strain_tensor=strain_tensor)
+
+        polycrystal.transform(motion, time=0.134)
+        polycrystal.diffract(
+            beam,
+            detector,
+            motion,
+            min_bragg_angle=0,
+            max_bragg_angle=None,
+            verbose=True)
+
+        diffraction_pattern = detector.render(
+            frames_to_render=0,
+            lorentz=False,
+            polarization=False,
+            structure_factor=False,
+            method="centroid",
+            verbose=True)
+        bins, histogram = utils._diffractogram(
+            diffraction_pattern, parameters['detector_center_pixel_z'], parameters['detector_center_pixel_y'])
+        histogram[histogram < 0.5 * np.median(histogram)] = 0
+        csequence, nosequences = 0, 0
+        for i in range(histogram.shape[0]):
+            if histogram[i] > 0:
+                csequence += 1
+            elif csequence >= 1:
+                nosequences += 1
+                csequence = 0
+        self.assertGreaterEqual(
+            nosequences,
+            10,
+            msg="Few or no rings appeared from diffraction.")
 
     def test_get_uniform_powder_sample(self):
         sample_bounding_radius = 256 * 180 / 128.

xrd_simulator/templates.py

@@ -94,6 +94,7 @@ def polycrystal_from_odf(
     path_to_cif_file=None,
     maximum_sampling_bin_seperation=np.radians(5.0),
     strain_tensor=lambda x: np.zeros((3, 3)),
+    phase_fractions=None
 ):
     """Fill a cylinder with crystals from a given orientation density function.
 
@@ -117,12 +118,19 @@ def polycrystal_from_odf(
              microns.
         sample_bounding_cylinder_radius (:obj:`float`): Radius of sample cylinder in units of
             microns.
-        unit_cell (:obj:`list` of :obj:`float`): Crystal unit cell representation of the form
-            [a,b,c,alpha,beta,gamma], where alpha,beta and gamma are in units of degrees while
-            a,b and c are in units of anstrom.
-        sgname (:obj:`string`): Name of space group , e.g 'P3221' for quartz, SiO2, for instance
-        path_to_cif_file (:obj:`string`): Path to CIF file. Defaults to None, in which case no structure
-            factors are computed.
+        unit_cell (:obj:`list of lists` of :obj:`float`): Crystal unit cell representation of the form
+            [a,b,c,alpha,beta,gamma] or [[a,b,c,alpha,beta,gamma],], where alpha,beta and gamma are 
+            in units of degrees while a,b and c are in units of anstrom. When the unit_cell is just a 
+            list (first example), the input represents single-phase material. When the unit_cell is an 
+            iterable list (second example), the input represents multi-phase material.
+        sgname (:obj:`list of strings`): Name of space group , e.g 'P3221' or ['P3221',] for quartz, 
+            SiO2, for instance. When the input is just a string, it represents space group for 
+            single-phase material. When the input is a list of string (second example), it represents 
+            space groups for multi-phase material.
+        path_to_cif_file (:obj:`list of strings`): Path to CIF file or [Path to CIF file,]. 
+            Defaults to None, in which case no structure factors are computed. When the input is a 
+            single file path (first example), the input is for single-phase material. When the 
+            input is a list of file paths (second example), the input is for multi-phase material.
         maximum_sampling_bin_seperation (:obj:`float`): Discretization steplength of orientation
             space using spherical coordinates over the unit quarternions in units of radians.
             A smaller steplength gives more accurate sampling of the input
@@ -131,6 +139,8 @@ def polycrystal_from_odf(
             strain_tensor(x) -> :obj:`numpy array` of shape ``(3,3)`` where input variable ``x`` is
             a :obj:`numpy array` of shape ``(3,)`` representing a spatial coordinate in the
             cylinder (x,y,z).
+        phase_fraction (:obj: `list` of :obj:`float`): Phase fraction represented as probability, summing up to 1.
+            Default None; will take even distribution of crystals.
 
     Returns:
         (:obj:`xrd_simulator.polycrystal.Polycrystal`)
@@ -168,10 +178,29 @@ def polycrystal_from_odf(
     )
 
     mesh = TetraMesh._build_tetramesh(cylinder)
-
-    # Sample is uniformly single phase
-    phases = [Phase(unit_cell, sgname, path_to_cif_file)]
-    element_phase_map = np.zeros((mesh.number_of_elements,)).astype(int)
+
+    if all(isinstance(x, list) for x in unit_cell) \
+        and isinstance(sgname, list):
+        # Multi Phase Material
+        if not isinstance(path_to_cif_file, list):
+            if path_to_cif_file is None:
+                path_to_cif_file = [None]* len(unit_cell)
+            else:
+                print("Single cif file input. Please enter list of corresponding cif files.")
+                exit
+        phases = [Phase(uc, sg, cf) for uc, sg, cf in zip(unit_cell, sgname, path_to_cif_file)]
+        if phase_fractions is None:
+            # Sample is uniformly distributed phases
+            element_phase_map = np.random.randint(len(phases), size=(mesh.number_of_elements,))
+        else :
+            # Sample is distributed by phase fraction
+            probabilities = np.array(phase_fractions)
+            element_phase_map = np.random.choice(len(phases), size=(mesh.number_of_elements,), p=probabilities)
+    else:
+        # Single Phase Material
+        phases = [Phase(unit_cell, sgname, path_to_cif_file)]
+        # Sample is uniformly single phase
+        element_phase_map = np.zeros((mesh.number_of_elements,)).astype(int)
 
     # Sample spatial texture
     orientation = _sample_ODF(