diff --git a/docs/source/examples/example_polycrystal_from_odf_2phases.py b/docs/source/examples/example_polycrystal_from_odf_2phases.py new file mode 100644 index 0000000..c77b704 --- /dev/null +++ b/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) diff --git a/tests/test_templates.py b/tests/test_templates.py index 3911e38..3e2d68f 100644 --- a/tests/test_templates.py +++ b/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. diff --git a/xrd_simulator/templates.py b/xrd_simulator/templates.py index 37268dd..6e7750c 100644 --- a/xrd_simulator/templates.py +++ b/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(