diff --git a/hexrd/constants.py b/hexrd/constants.py index aa813e52..932814ae 100644 --- a/hexrd/constants.py +++ b/hexrd/constants.py @@ -82,7 +82,8 @@ zeros_3x1 = np.zeros((3, 1)) zeros_6x1 = np.zeros((6, 1)) -# reference beam direction and eta=0 ref in LAB FRAME for standard geometry +'''reference beam direction and +eta=0 ref in LAB FRAME for standard geometry''' beam_vec = -lab_z eta_vec = lab_x @@ -162,7 +163,9 @@ 3628800.]).astype(np.complex128) """ ->> @AUTHOR: Saransh Singh, Lawrence Livermore National Lab, saransh1@llnl.gov +>> @AUTHOR: Saransh Singh, + Lawrence Livermore National Lab, + saransh1@llnl.gov >> @DATE: 11/28/2022 SS 1.0 original >> @DETAILS: constants for rodrigues FZ """ @@ -177,7 +180,9 @@ ) ''' ->> @AUTHOR: Saransh Singh, Lawrence Livermore National Lab, saransh1@llnl.gov +>> @AUTHOR: Saransh Singh, + Lawrence Livermore National Lab, + saransh1@llnl.gov >> @DATE: 10/28/2020 SS 1.0 original >> @DETAILS: constants for sphere sectors used for IPF coloring ''' @@ -263,7 +268,8 @@ def set_numba_cache(): 4. The NUMBA_CACHE_DIR environment variable is not defined If all of these are true, then numba will try to write to a - directory where it doesn't have permission, and cause the application to + directory where it doesn't have permission, + and cause the application to freeze. Avoid that by setting the cache dir ourselves. """ @@ -316,8 +322,10 @@ def is_writable_file(path): cRestmass = 9.1093837090E-31 ''' -adding another parametrization of the scattering factors. these are more recent -and more accurate. also used in Vesta (copied from there). see: +adding another parametrization of the +scattering factors. these are more recent +and more accurate. also used in Vesta +(copied from there). see: New Analytical coherent Scattering-Factor Functions for Free Atoms and Ions BY D. WAASMAIER AND A. KIRFEL @@ -638,7 +646,9 @@ def is_writable_file(path): 'Cf': ['0'] } ''' -this dictionary tabulates the small nuclear Thomson term fNT for all elements up to Z=92 +this dictionary tabulates the small +nuclear Thomson term fNT for all +elements up to Z=92 ''' fNT = { 'H':-0.00054423,'He':-0.00054817,'Li':-0.00071131,'Be':-0.00097394,'B':-0.0012687,'C':-0.0016442,'N':-0.0019191,'O':-0.0021944, @@ -737,8 +747,10 @@ def is_writable_file(path): 229 ]) -''' this variable encodes all the generators (including translations) for all 230 space groups - will be used to compute the full space group symmetry operators +''' this variable encodes all the generators +(including translations) for all 230 space groups +will be used to compute the full space group symmetry +operators ''' SYM_GL = [ "000 ", "100 ", "01cOOO0 ", @@ -996,8 +1008,10 @@ def is_writable_file(path): } ''' -this dictionary contains the generators encoded in each letter of the generator string -the full symmetry is generated by the repeated action of the generator matrix +this dictionary contains the generators encoded +in each letter of the generator string +the full symmetry is generated by the repeated +action of the generator matrix ''' ''' rotational, inversions, mirrors etc. components @@ -1100,17 +1114,21 @@ def is_writable_file(path): SYM_GENERATORS['Y'] = -1./4. SYM_GENERATORS['Z'] = -1./8. - ''' - @AUTHOR Saransh Singh, Lawrence Livermore National Lab, saransh1@llnl.gov + @AUTHOR Saransh Singh, + Lawrence Livermore National Lab, + saransh1@llnl.gov @DATE 11/23/2020 SS 1.0 original - @DETAIL. this list of symbols will help us to genrate the point group symmetries - in the cartesian space for any point group. this is needed for the - supergroup symmetry usd in the coloring scheme used in the package. this - needs to be a separate set of routines because the supergroup can be a - point group which is not the laue group of the crystal (e.g. m-3 --> m-3m) - the notation used will be the same as the one used for the space group - without any translations. + @DETAIL. this list of symbols will help us to genrate + the point group symmetries in the cartesian + space for any point group. this is needed for + the supergroup symmetry usd in the coloring + scheme used in the package. this needs to be a + separate set of routines because the supergroup + can be a point group which is not the laue group + of the crystal (e.g. m-3 --> m-3m) the notation + used will be the same as the one used for the + space group without any translations. ''' SYM_GL_PG = { 'c1': '1a', # only identity rotation diff --git a/hexrd/findorientations.py b/hexrd/findorientations.py index b04fe4fd..aa24aeb8 100755 --- a/hexrd/findorientations.py +++ b/hexrd/findorientations.py @@ -88,7 +88,7 @@ def generate_orientation_fibers(cfg, eta_ome): pd = eta_ome.planeData tTh = pd.getTTh() bMat = pd.latVecOps['B'] - csym = pd.getLaueGroup() + csym = pd.laueGroup # !!! changed recently where iHKLList are now master hklIDs pd_hkl_ids = eta_ome.iHKLList[seed_hkl_ids] @@ -885,8 +885,8 @@ def find_orientations(cfg, logger.info("\tmean reflections per grain: %d", mean_rpg) logger.info("\tneighborhood size: %d", min_samples) - qbar, _ = run_cluster( - completeness, qfib, plane_data.getQSym(), cfg, + qbar, cl = run_cluster( + completeness, qfib, plane_data.q_sym, cfg, min_samples=min_samples, compl_thresh=compl_thresh, radius=cl_radius) diff --git a/hexrd/instrument/detector.py b/hexrd/instrument/detector.py index 9071a489..227f2037 100644 --- a/hexrd/instrument/detector.py +++ b/hexrd/instrument/detector.py @@ -1131,14 +1131,7 @@ def make_powder_rings( # in case you want to give it tth angles directly if isinstance(pd, PlaneData): - # Okay, we have a PlaneData object - try: - pd = pd.makeNew() # make a copy to munge - except TypeError: - # !!! have some other object here, - # likely a dummy plane data object of sorts - raise - + pd = PlaneData(None, pd) if delta_tth is not None: pd.tThWidth = np.radians(delta_tth) else: diff --git a/hexrd/material/crystallography.py b/hexrd/material/crystallography.py index 76d74fd7..13148d29 100644 --- a/hexrd/material/crystallography.py +++ b/hexrd/material/crystallography.py @@ -28,16 +28,21 @@ import re import copy from math import pi +from typing import Optional, Union, Dict, List, Tuple import numpy as np -import csv -import os +from hexrd.material.unitcell import unitcell +from hexrd.deprecation import deprecated from hexrd import constants from hexrd.matrixutil import unitVector -from hexrd.rotations import \ - rotMatOfExpMap, mapAngle, applySym, \ - ltypeOfLaueGroup, quatOfLaueGroup +from hexrd.rotations import ( + rotMatOfExpMap, + mapAngle, + applySym, + ltypeOfLaueGroup, + quatOfLaueGroup, +) from hexrd.transforms import xfcapi from hexrd import valunits from hexrd.valunits import toFloat @@ -51,24 +56,60 @@ outputDegrees_bak = outputDegrees -def hklToStr(x): - return re.sub(r'[\[\]\(\)\{\},]', '', str(x)) +def hklToStr(hkl: np.ndarray) -> str: + """ + Converts hkl representation to a string. + + Parameters + ---------- + hkl : np.ndarray + 3 element list of h, k, and l values (Miller indices). + + Returns + ------- + str + Space-separated string representation of h, k, and l values. + + """ + return re.sub(r'[\[\]\(\)\{\},]', '', str(hkl)) + + +def tempSetOutputDegrees(val: bool) -> None: + """ + Set the global outputDegrees flag temporarily. Can be reverted with + revertOutputDegrees(). + + Parameters + ---------- + val : bool + True to output angles in degrees, False to output angles in radians. + Returns + ------- + None -def tempSetOutputDegrees(val): + """ global outputDegrees, outputDegrees_bak outputDegrees_bak = outputDegrees outputDegrees = val - return -def revertOutputDegrees(): +def revertOutputDegrees() -> None: + """ + Revert the effect of tempSetOutputDegrees(), resetting the outputDegrees + flag to its previous value (True to output in degrees, False for radians). + + Returns + ------- + None + """ global outputDegrees, outputDegrees_bak outputDegrees = outputDegrees_bak - return -def cosineXform(a, b, c): +def cosineXform( + a: np.ndarray, b: np.ndarray, c: np.ndarray +) -> tuple[np.ndarray, np.ndarray]: """ Spherical trig transform to take alpha, beta, gamma to expressions for cos(alpha*). See ref below. @@ -77,67 +118,59 @@ def cosineXform(a, b, c): the relations between direct and reciprocal lattice quantities''. Acta Cryst. (1968), A24, 247--248 + Parameters + ---------- + a : np.ndarray + List of alpha angle values (radians). + b : np.ndarray + List of beta angle values (radians). + c : np.ndarray + List of gamma angle values (radians). + + Returns + ------- + np.ndarray + List of cos(alpha*) values. + np.ndarray + List of sin(alpha*) values. + """ - cosar = (np.cos(b)*np.cos(c) - np.cos(a)) / (np.sin(b)*np.sin(c)) + cosar = (np.cos(b) * np.cos(c) - np.cos(a)) / (np.sin(b) * np.sin(c)) sinar = np.sqrt(1 - cosar**2) return cosar, sinar -def processWavelength(arg): +def processWavelength(arg: Union[valunits.valWUnit, float]) -> float: """ Convert an energy value to a wavelength. If argument has units of length or energy, will convert to globally specified unit type for wavelength (dUnit). If argument is a scalar, assumed input units are keV. """ - if hasattr(arg, 'getVal'): + if isinstance(arg, valunits.valWUnit): + # arg is a valunits.valWUnit object if arg.isLength(): - retval = arg.getVal(dUnit) + return arg.getVal(dUnit) elif arg.isEnergy(): e = arg.getVal('keV') - retval = valunits.valWUnit( + return valunits.valWUnit( 'wavelength', 'length', constants.keVToAngstrom(e), 'angstrom' ).getVal(dUnit) else: - raise RuntimeError('do not know what to do with '+str(arg)) + raise RuntimeError('do not know what to do with ' + str(arg)) else: # !!! assuming arg is in keV - retval = valunits.valWUnit( + return valunits.valWUnit( 'wavelength', 'length', constants.keVToAngstrom(arg), 'angstrom' ).getVal(dUnit) - return retval - - -def latticeParameters(lvec): - """ - Generates direct and reciprocal lattice vector components in a - crystal-relative RHON basis, X. The convention for fixing X to the - lattice is such that a || x1 and c* || x3, where a and c* are - direct and reciprocal lattice vectors, respectively. - """ - lnorm = np.sqrt(np.sum(lvec**2, 0)) - - a = lnorm[0] - b = lnorm[1] - c = lnorm[2] - - ahat = lvec[:, 0]/a - bhat = lvec[:, 1]/b - chat = lvec[:, 2]/c - gama = np.arccos(np.dot(ahat, bhat)) - beta = np.arccos(np.dot(ahat, chat)) - alfa = np.arccos(np.dot(bhat, chat)) - if outputDegrees: - gama = r2d*gama - beta = r2d*beta - alfa = r2d*alfa - - return [a, b, c, alfa, beta, gama] - - -def latticePlanes(hkls, lparms, - ltype='cubic', wavelength=1.54059292, strainMag=None): +def latticePlanes( + hkls: np.ndarray, + lparms: np.ndarray, + ltype: Optional[str] = 'cubic', + wavelength: Optional[float] = 1.54059292, + strainMag: Optional[float] = None, +) -> Dict[str, np.ndarray]: """ Generates lattice plane data in the direct lattice for a given set of Miller indices. Vector components are written in the @@ -158,9 +191,9 @@ def latticePlanes(hkls, lparms, 2) lparms (1 x m float list) is the array of lattice parameters, where m depends on the symmetry group (see below). - 3) The following optional keyword arguments are recognized: + The following optional arguments are recognized: - *) ltype=(string) is a string representing the symmetry type of + 3) ltype=(string) is a string representing the symmetry type of the implied Laue group. The 11 available choices are shown below. The default value is 'cubic'. Note that each group expects a lattice parameter array of the indicated length @@ -177,11 +210,11 @@ def latticePlanes(hkls, lparms, 'monoclinic' a, b, c, beta (in degrees) 'triclinic' a, b, c, alpha, beta, gamma (in degrees) - *) wavelength= is a value represented the wavelength in + 4) wavelength= is a value represented the wavelength in Angstroms to calculate bragg angles for. The default value is for Cu K-alpha radiation (1.54059292 Angstrom) - *) strainMag=None + 5) strainMag=None OUTPUTS: @@ -194,7 +227,7 @@ def latticePlanes(hkls, lparms, dspacings (n, ) double array array of the d-spacings for each {hkl} - 2thetas (n, ) double array array of the Bragg angles for + tThetas (n, ) double array array of the Bragg angles for each {hkl} relative to the specified wavelength @@ -216,8 +249,9 @@ def latticePlanes(hkls, lparms, """ location = 'latticePlanes' - assert hkls.shape[0] == 3, \ - "hkls aren't column vectors in call to '%s'!" % location + assert ( + hkls.shape[0] == 3 + ), f"hkls aren't column vectors in call to '{location}'!" tag = ltype wlen = wavelength @@ -231,43 +265,46 @@ def latticePlanes(hkls, lparms, # magnitudes d = 1 / np.sqrt(np.sum(G**2, 0)) - aconv = 1. + aconv = 1.0 if outputDegrees: aconv = r2d # two thetas - sth = wlen / 2. / d - mask = (np.abs(sth) < 1.) + sth = wlen / 2.0 / d + mask = np.abs(sth) < 1.0 tth = np.zeros(sth.shape) tth[~mask] = np.nan - tth[mask] = aconv * 2. * np.arcsin(sth[mask]) + tth[mask] = aconv * 2.0 * np.arcsin(sth[mask]) - p = dict(normals=unitVector(G), - dspacings=d, - tThetas=tth) + p = dict(normals=unitVector(G), dspacings=d, tThetas=tth) if strainMag is not None: p['tThetasLo'] = np.zeros(sth.shape) p['tThetasHi'] = np.zeros(sth.shape) - mask = ( - (np.abs(wlen / 2. / (d * (1. + strainMag))) < 1.) & - (np.abs(wlen / 2. / (d * (1. - strainMag))) < 1.) + mask = (np.abs(wlen / 2.0 / (d * (1.0 + strainMag))) < 1.0) & ( + np.abs(wlen / 2.0 / (d * (1.0 - strainMag))) < 1.0 ) p['tThetasLo'][~mask] = np.nan p['tThetasHi'][~mask] = np.nan - p['tThetasLo'][mask] = \ - aconv * 2 * np.arcsin(wlen/2./(d[mask]*(1. + strainMag))) - p['tThetasHi'][mask] = \ - aconv * 2 * np.arcsin(wlen/2./(d[mask]*(1. - strainMag))) + p['tThetasLo'][mask] = ( + aconv * 2 * np.arcsin(wlen / 2.0 / (d[mask] * (1.0 + strainMag))) + ) + p['tThetasHi'][mask] = ( + aconv * 2 * np.arcsin(wlen / 2.0 / (d[mask] * (1.0 - strainMag))) + ) return p -def latticeVectors(lparms, tag='cubic', radians=False, debug=False): +def latticeVectors( + lparms: np.ndarray, + tag: Optional[str] = 'cubic', + radians: Optional[bool] = False, +) -> Dict[str, Union[np.ndarray, float]]: """ Generates direct and reciprocal lattice vector components in a crystal-relative RHON basis, X. The convention for fixing X to the @@ -283,7 +320,7 @@ def latticeVectors(lparms, tag='cubic', radians=False, debug=False): 1) lparms (1 x n float list) is the array of lattice parameters, where n depends on the symmetry group (see below). - 2) symTag (string) is a case-insensitive string representing the + 2) tag (string) is a case-insensitive string representing the symmetry type of the implied Laue group. The 11 available choices are shown below. The default value is 'cubic'. Note that each group expects a lattice parameter array of the indicated length @@ -300,6 +337,11 @@ def latticeVectors(lparms, tag='cubic', radians=False, debug=False): 'monoclinic' a, b, c, beta (in degrees) 'triclinic' a, b, c, alpha, beta, gamma (in degrees) + The following optional arguments are recognized: + + 3) radians= is a boolean flag indicating usage of radians rather + than degrees, defaults to false. + OUTPUTS: 1) lattice is a dictionary containing the following keys/items: @@ -380,60 +422,64 @@ def latticeVectors(lparms, tag='cubic', radians=False, debug=False): 'tetragonal', 'orthorhombic', 'monoclinic', - 'triclinic' - ] + 'triclinic', + ] if radians: - aconv = 1. + aconv = 1.0 else: - aconv = pi/180. # degToRad - deg90 = pi/2. - deg120 = 2.*pi/3. + aconv = pi / 180.0 # degToRad + deg90 = pi / 2.0 + deg120 = 2.0 * pi / 3.0 # if tag == lattStrings[0]: # cubic - cellparms = np.r_[np.tile(lparms[0], (3,)), deg90*np.ones((3,))] + cellparms = np.r_[np.tile(lparms[0], (3,)), deg90 * np.ones((3,))] elif tag == lattStrings[1] or tag == lattStrings[2]: # hexagonal | trigonal (hex indices) - cellparms = np.r_[lparms[0], lparms[0], lparms[1], - deg90, deg90, deg120] + cellparms = np.r_[ + lparms[0], lparms[0], lparms[1], deg90, deg90, deg120 + ] elif tag == lattStrings[3]: # rhombohedral - cellparms = np.r_[np.tile(lparms[0], (3,)), - np.tile(aconv*lparms[1], (3,))] + cellparms = np.r_[ + np.tile(lparms[0], (3,)), np.tile(aconv * lparms[1], (3,)) + ] elif tag == lattStrings[4]: # tetragonal - cellparms = np.r_[lparms[0], lparms[0], lparms[1], - deg90, deg90, deg90] + cellparms = np.r_[lparms[0], lparms[0], lparms[1], deg90, deg90, deg90] elif tag == lattStrings[5]: # orthorhombic - cellparms = np.r_[lparms[0], lparms[1], lparms[2], - deg90, deg90, deg90] + cellparms = np.r_[lparms[0], lparms[1], lparms[2], deg90, deg90, deg90] elif tag == lattStrings[6]: # monoclinic - cellparms = np.r_[lparms[0], lparms[1], lparms[2], - deg90, aconv*lparms[3], deg90] + cellparms = np.r_[ + lparms[0], lparms[1], lparms[2], deg90, aconv * lparms[3], deg90 + ] elif tag == lattStrings[7]: # triclinic - # FIXME: fixed DP 2/24/16 - cellparms = np.r_[lparms[0], lparms[1], lparms[2], - aconv*lparms[3], aconv*lparms[4], aconv*lparms[5]] + cellparms = np.r_[ + lparms[0], + lparms[1], + lparms[2], + aconv * lparms[3], + aconv * lparms[4], + aconv * lparms[5], + ] else: - raise RuntimeError('lattice tag \'%s\' is not recognized' % (tag)) - - if debug: - print((str(cellparms[0:3]) + ' ' + str(r2d*cellparms[3:6]))) - alfa = cellparms[3] - beta = cellparms[4] - gama = cellparms[5] + raise RuntimeError(f'lattice tag "{tag}" is not recognized') - cosalfar, sinalfar = cosineXform(alfa, beta, gama) + alpha, beta, gamma = cellparms[3:6] + cosalfar, sinalfar = cosineXform(alpha, beta, gamma) - a = cellparms[0]*np.r_[1, 0, 0] - b = cellparms[1]*np.r_[np.cos(gama), np.sin(gama), 0] - c = cellparms[2]*np.r_[np.cos(beta), - -cosalfar*np.sin(beta), - sinalfar*np.sin(beta)] + a = cellparms[0] * np.r_[1, 0, 0] + b = cellparms[1] * np.r_[np.cos(gamma), np.sin(gamma), 0] + c = ( + cellparms[2] + * np.r_[ + np.cos(beta), -cosalfar * np.sin(beta), sinalfar * np.sin(beta) + ] + ) ad = np.sqrt(np.sum(a**2)) bd = np.sqrt(np.sum(b**2)) @@ -446,200 +492,62 @@ def latticeVectors(lparms, tag='cubic', radians=False, debug=False): F = np.c_[a, b, c] # Reciprocal lattice vectors - astar = np.cross(b, c)/V - bstar = np.cross(c, a)/V - cstar = np.cross(a, b)/V + astar = np.cross(b, c) / V + bstar = np.cross(c, a) / V + cstar = np.cross(a, b) / V # and parameters ar = np.sqrt(np.sum(astar**2)) br = np.sqrt(np.sum(bstar**2)) cr = np.sqrt(np.sum(cstar**2)) - alfar = np.arccos(np.dot(bstar, cstar)/br/cr) - betar = np.arccos(np.dot(cstar, astar)/cr/ar) - gamar = np.arccos(np.dot(astar, bstar)/ar/br) + alfar = np.arccos(np.dot(bstar, cstar) / br / cr) + betar = np.arccos(np.dot(cstar, astar) / cr / ar) + gamar = np.arccos(np.dot(astar, bstar) / ar / br) # B takes components in the reciprocal lattice to X B = np.c_[astar, bstar, cstar] cosalfar2, sinalfar2 = cosineXform(alfar, betar, gamar) - afable = ar*np.r_[1, 0, 0] - bfable = br*np.r_[np.cos(gamar), np.sin(gamar), 0] - cfable = cr*np.r_[np.cos(betar), - -cosalfar2*np.sin(betar), - sinalfar2*np.sin(betar)] + afable = ar * np.r_[1, 0, 0] + bfable = br * np.r_[np.cos(gamar), np.sin(gamar), 0] + cfable = ( + cr + * np.r_[ + np.cos(betar), + -cosalfar2 * np.sin(betar), + sinalfar2 * np.sin(betar), + ] + ) BR = np.c_[afable, bfable, cfable] U0 = np.dot(B, np.linalg.inv(BR)) if outputDegrees: - dparms = np.r_[ad, bd, cd, r2d*np.r_[alfa, beta, gama]] - rparms = np.r_[ar, br, cr, r2d*np.r_[alfar, betar, gamar]] + dparms = np.r_[ad, bd, cd, r2d * np.r_[alpha, beta, gamma]] + rparms = np.r_[ar, br, cr, r2d * np.r_[alfar, betar, gamar]] else: - dparms = np.r_[ad, bd, cd, np.r_[alfa, beta, gama]] + dparms = np.r_[ad, bd, cd, np.r_[alpha, beta, gamma]] rparms = np.r_[ar, br, cr, np.r_[alfar, betar, gamar]] - L = {'F': F, - 'B': B, - 'BR': BR, - 'U0': U0, - 'vol': V, - 'dparms': dparms, - 'rparms': rparms} - - return L - - -def hexagonalIndicesFromRhombohedral(hkl): - """ - converts rhombohedral hkl to hexagonal indices - """ - HKL = np.zeros((3, hkl.shape[1]), dtype='int') - - HKL[0, :] = hkl[0, :] - hkl[1, :] - HKL[1, :] = hkl[1, :] - hkl[2, :] - HKL[2, :] = hkl[0, :] + hkl[1, :] + hkl[2, :] - - return HKL - - -def rhombohedralIndicesFromHexagonal(HKL): - """ - converts hexagonal hkl to rhombohedral indices - """ - hkl = np.zeros((3, HKL.shape[1]), dtype='int') - - hkl[0, :] = 2 * HKL[0, :] + HKL[1, :] + HKL[2, :] - hkl[1, :] = -HKL[0, :] + HKL[1, :] + HKL[2, :] - hkl[2, :] = -HKL[0, :] - 2 * HKL[1, :] + HKL[2, :] - - hkl = hkl / 3. - return hkl - - -def rhombohedralParametersFromHexagonal(a_h, c_h): - """ - converts hexagonal lattice parameters (a, c) to rhombohedral - lattice parameters (a, alpha) - """ - a_r = np.sqrt(3 * a_h**2 + c_h**2) / 3. - alfa_r = 2 * np.arcsin(3. / (2 * np.sqrt(3 + (c_h / a_h)**2))) - if outputDegrees: - alfa_r = r2d * alfa_r - return a_r, alfa_r - - -def convert_Miller_direction_to_cartesian(uvw, a=1., c=1., normalize=False): - """ - Converts 3-index hexagonal Miller direction indices to components in the - crystal reference frame. - - Parameters - ---------- - uvw : array_like - The (n, 3) array of 3-index hexagonal indices to convert. - a : scalar, optional - The `a` lattice parameter. The default value is 1. - c : scalar, optional - The `c` lattice parameter. The default value is 1. - normalize : bool, optional - Flag for whether or not to normalize output vectors - - Returns - ------- - numpy.ndarray - The (n, 3) array of cartesian components associated with the input - direction indices. - - Notes - ----- - 1) The [uv.w] the Miller-Bravais convention is in the hexagonal basis - {a1, a2, a3, c}. The basis for the output, {o1, o2, o3}, is - chosen such that - - o1 || a1 - o3 || c - o2 = o3 ^ o1 - - """ - u, v, w = np.atleast_2d(uvw).T - retval = np.vstack([1.5*u*a, sqrt3by2*a*(2*v + u), w*c]) - if normalize: - return unitVector(retval).T - else: - return retval.T - - -def convert_Miller_direction_to_MillerBravias(uvw, suppress_redundant=True): - """ - Converts 3-index hexagonal Miller direction indices to 4-index - Miller-Bravais direction indices. - - Parameters - ---------- - uvw : array_like - The (n, 3) array of 3-index hexagonal Miller indices to convert. - suppress_redundant : bool, optional - Flag to suppress the redundant 3rd index. The default is True. - - Returns - ------- - numpy.ndarray - The (n, 3) or (n, 4) array -- depending on kwarg -- of Miller-Bravis - components associated with the input Miller direction indices. - - Notes - ----- - * NOT for plane normals!!! - - """ - u, v, w = np.atleast_2d(uvw).T - retval = np.vstack([(2*u - v)/3, (2*v - u)/3, w]).T - rem = np.vstack([np.mod(np.tile(i[0], 2), i[1:]) for i in retval]) - rem[abs(rem) < epsf] = np.nan - lcm = np.nanmin(rem, axis=1) - lcm[np.isnan(lcm)] = 1 - retval = retval / np.tile(lcm, (3, 1)).T - if suppress_redundant: - return retval - else: - t = np.atleast_2d(1 - np.sum(retval[:2], axis=1)).T - return np.hstack([retval[:, :2], t, np.atleast_2d(retval[:, 2]).T]) - - -def convert_MillerBravias_direction_to_Miller(UVW): - """ - Converts 4-index hexagonal Miller-Bravais direction indices to - 3-index Miller direction indices. - - Parameters - ---------- - UVW : array_like - The (n, 3) array of **non-redundant** Miller-Bravais direction indices - to convert. - - Returns - ------- - numpy.ndarray - The (n, 3) array of Miller direction indices associated with the - input Miller-Bravais indices. - - Notes - ----- - * NOT for plane normals!!! - - """ - U, V, W = np.atleast_2d(UVW).T - return np.vstack([2*U + V, 2*V + U, W]) + return { + 'F': F, + 'B': B, + 'BR': BR, + 'U0': U0, + 'vol': V, + 'dparms': dparms, + 'rparms': rparms, + } class PlaneData(object): """ Careful with ordering: Outputs are ordered by the 2-theta for the - hkl unless you get self.__hkls directly, and this order can change + hkl unless you get self._hkls directly, and this order can change with changes in lattice parameters (lparms); setting and getting exclusions works on the current hkl ordering, not the original - ordering (in self.__hkls), but exclusions are stored in the + ordering (in self._hkls), but exclusions are stored in the original ordering in case the hkl ordering does change with lattice parameters @@ -648,77 +556,94 @@ class PlaneData(object): tThWidth """ - def __init__(self, - hkls, - *args, - **kwargs): + def __init__(self, hkls: Optional[np.ndarray], *args, **kwargs) -> None: + """ + Constructor for PlaneData + + Parameters + ---------- + hkls : np.ndarray + Miller indices to be used in the plane data. Can be None if + args is another PlaneData object + + *args + Unnamed arguments. Could be in the format of `lparms, laueGroup, + wavelength, strainMag`, or just a `PlaneData` object. + + **kwargs + Valid keyword arguments include: + - doTThSort + - exclusions + - tThMax + - tThWidth + """ + self._doTThSort = True + self._exclusions = None + self._tThMax = None - self.phaseID = None - self.__doTThSort = True - self.__exclusions = None - self.__tThMax = None - # if len(args) == 4: lparms, laueGroup, wavelength, strainMag = args tThWidth = None - self.__wavelength = processWavelength(wavelength) - self.__lparms = self.__parseLParms(lparms) - elif len(args) == 1 and hasattr(args[0], 'getParams'): + self._wavelength = processWavelength(wavelength) + self._lparms = self._parseLParms(lparms) + elif len(args) == 1 and isinstance(args[0], PlaneData): other = args[0] - lparms, laueGroup, wavelength, strainMag, tThWidth = \ + lparms, laueGroup, wavelength, strainMag, tThWidth = ( other.getParams() - self.__wavelength = wavelength - self.__lparms = lparms - self.phaseID = other.phaseID - self.__doTThSort = other.__doTThSort - self.__exclusions = other.__exclusions - self.__tThMax = other.__tThMax + ) + self._wavelength = wavelength + self._lparms = lparms + self._doTThSort = other._doTThSort + self._exclusions = other._exclusions + self._tThMax = other._tThMax if hkls is None: - hkls = other.__hkls + hkls = other._hkls else: - raise NotImplementedError('args : '+str(args)) + raise NotImplementedError(f'args : {args}') - self.__laueGroup = laueGroup - self.__qsym = quatOfLaueGroup(self.__laueGroup) - self.__hkls = copy.deepcopy(hkls) - self.__strainMag = strainMag - self.__structFact = np.ones(self.__hkls.shape[1]) + self._laueGroup = laueGroup + self._hkls = copy.deepcopy(hkls) + self._strainMag = strainMag + self._structFact = np.ones(self._hkls.shape[1]) self.tThWidth = tThWidth # ... need to implement tThMin too - if 'phaseID' in kwargs: - self.phaseID = kwargs.pop('phaseID') if 'doTThSort' in kwargs: - self.__doTThSort = kwargs.pop('doTThSort') + self._doTThSort = kwargs.pop('doTThSort') if 'exclusions' in kwargs: - self.__exclusions = kwargs.pop('exclusions') + self._exclusions = kwargs.pop('exclusions') if 'tThMax' in kwargs: - self.__tThMax = toFloat(kwargs.pop('tThMax'), 'radians') + self._tThMax = toFloat(kwargs.pop('tThMax'), 'radians') if 'tThWidth' in kwargs: self.tThWidth = kwargs.pop('tThWidth') if len(kwargs) > 0: - raise RuntimeError('have unparsed keyword arguments with keys: ' - + str(list(kwargs.keys()))) + raise RuntimeError( + f'have unparsed keyword arguments with keys: {kwargs.keys()}' + ) # This is only used to calculate the structure factor if invalidated - self.__unitcell = None - - self.__calc() - - return - - def __calc(self): - symmGroup = ltypeOfLaueGroup(self.__laueGroup) - latPlaneData, latVecOps, hklDataList = PlaneData.makePlaneData( - self.__hkls, self.__lparms, self.__qsym, - symmGroup, self.__strainMag, self.wavelength) + self._unitcell: unitcell = None + + self._calc() + + def _calc(self): + symmGroup = ltypeOfLaueGroup(self._laueGroup) + self._q_sym = quatOfLaueGroup(self._laueGroup) + _, latVecOps, hklDataList = PlaneData.makePlaneData( + self._hkls, + self._lparms, + self._q_sym, + symmGroup, + self._strainMag, + self.wavelength, + ) 'sort by tTheta' tThs = np.array( [hklDataList[iHKL]['tTheta'] for iHKL in range(len(hklDataList))] ) - if self.__doTThSort: - # sorted hkl -> __hkl - # __hkl -> sorted hkl + if self._doTThSort: + # sorted hkl -> _hkl + # _hkl -> sorted hkl self.tThSort = np.argsort(tThs) self.tThSortInv = np.empty(len(hklDataList), dtype=int) self.tThSortInv[self.tThSort] = np.arange(len(hklDataList)) @@ -729,101 +654,131 @@ def __calc(self): self.hklDataList = hklDataList self._latVecOps = latVecOps self.nHKLs = len(self.getHKLs()) - return def __str__(self): s = '========== plane data ==========\n' s += 'lattice parameters:\n ' + str(self.lparms) + '\n' - s += 'two theta width: (%s)\n' % str(self.tThWidth) - s += 'strain magnitude: (%s)\n' % str(self.strainMag) - s += 'beam energy (%s)\n' % str(self.wavelength) + s += f'two theta width: ({str(self.tThWidth)})\n' + s += f'strain magnitude: ({str(self.strainMag)})\n' + s += f'beam energy ({str(self.wavelength)})\n' s += 'hkls: (%d)\n' % self.nHKLs s += str(self.getHKLs()) return s - def getNHKLs(self): - return self.nHKLs + def getParams(self): + """ + Getter for the parameters of the plane data. + + Returns + ------- + tuple + The parameters of the plane data. In the order of + _lparams, _laueGroup, _wavelength, _strainMag, tThWidth - def getPhaseID(self): - 'may return None if not set' - return self.phaseID + """ + return ( + self._lparms, + self._laueGroup, + self._wavelength, + self._strainMag, + self.tThWidth, + ) - def getParams(self): - return (self.__lparms, self.__laueGroup, self.__wavelength, - self.__strainMag, self.tThWidth) + def getNhklRef(self) -> int: + """ + Get the total number of hkl's in the plane data, not ignoring + ones that are excluded in exclusions. - def getNhklRef(self): - 'does not use exclusions or the like' - retval = len(self.hklDataList) - return retval + Returns + ------- + int + The total number of hkl's in the plane data. + """ + return len(self.hklDataList) - def get_hkls(self): + @property + def hkls(self) -> np.ndarray: """ - do not do return self.__hkls, as everywhere else hkls are returned - in 2-theta order; transpose is to comply with lparm convention + hStacked Hkls of the plane data (Miller indices). """ return self.getHKLs().T - def set_hkls(self, hkls): - raise RuntimeError('for now, not allowing hkls to be reset') - # self.__exclusions = None - # self.__hkls = hkls - # self.__calc() - return - hkls = property(get_hkls, set_hkls, None) + @hkls.setter + def hkls(self, hkls): + raise NotImplementedError('for now, not allowing hkls to be reset') + + @property + def tThMax(self) -> Optional[float]: + """ + Maximum 2-theta value of the plane data. - def get_tThMax(self): - return self.__tThMax + float or None + """ + return self._tThMax - def set_tThMax(self, tThMax): - self.__tThMax = toFloat(tThMax, 'radians') - # self.__calc() # no need to redo calc for tThMax - return - tThMax = property(get_tThMax, set_tThMax, None) + @tThMax.setter + def tThMax(self, t_th_max: Union[float, valunits.valWUnit]) -> None: + self._tThMax = toFloat(t_th_max, 'radians') - def get_exclusions(self): + @property + def exclusions(self) -> np.ndarray: + """ + Excluded HKL's the plane data. + + Set as type np.ndarray, as a mask of length getNhklRef(), a list of + indices to be excluded, or a list of ranges of indices. + + Read as a mask of length getNhklRef(). + """ retval = np.zeros(self.getNhklRef(), dtype=bool) - if self.__exclusions is not None: + if self._exclusions is not None: # report in current hkl ordering - retval[:] = self.__exclusions[self.tThSortInv] - if self.__tThMax is not None: + retval[:] = self._exclusions[self.tThSortInv] + if self._tThMax is not None: for iHKLr, hklData in enumerate(self.hklDataList): - if hklData['tTheta'] > self.__tThMax: + if hklData['tTheta'] > self._tThMax: retval[iHKLr] = True return retval - def set_exclusions(self, exclusions): + @exclusions.setter + def exclusions(self, new_exclusions: Optional[np.ndarray]) -> None: excl = np.zeros(len(self.hklDataList), dtype=bool) - if exclusions is not None: - exclusions = np.atleast_1d(exclusions) + if new_exclusions is not None: + exclusions = np.atleast_1d(new_exclusions) if len(exclusions) == len(self.hklDataList): - assert exclusions.dtype == 'bool', \ - 'exclusions should be bool if full length' - # convert from current hkl ordering to __hkl ordering + assert ( + exclusions.dtype == 'bool' + ), 'Exclusions should be bool if full length' + # convert from current hkl ordering to _hkl ordering excl[:] = exclusions[self.tThSort] else: if len(exclusions.shape) == 1: # treat exclusions as indices excl[self.tThSort[exclusions]] = True elif len(exclusions.shape) == 2: - raise NotImplementedError( - 'have not yet coded treating exclusions as ranges' - ) + # treat exclusions as ranges of indices + for r in exclusions: + excl[self.tThSort[r[0]:r[1]]] = True else: raise RuntimeError( - 'do not now what to do with exclusions with shape ' - + str(exclusions.shape) + f'Unclear behavior for shape {exclusions.shape}' ) - self.__exclusions = excl - self.nHKLs = np.sum(np.logical_not(self.__exclusions)) - return - exclusions = property(get_exclusions, set_exclusions, None) + self._exclusions = excl + self.nHKLs = np.sum(np.logical_not(self._exclusions)) def exclude( - self, dmin=None, dmax=None, tthmin=None, tthmax=None, - sfacmin=None, sfacmax=None, pintmin=None, pintmax=None - ): - """Set exclusions according to various parameters + self, + dmin: Optional[float] = None, + dmax: Optional[float] = None, + tthmin: Optional[float] = None, + tthmax: Optional[float] = None, + sfacmin: Optional[float] = None, + sfacmax: Optional[float] = None, + pintmin: Optional[float] = None, + pintmax: Optional[float] = None, + ) -> None: + """ + Set exclusions according to various parameters Any hkl with a value below any min or above any max will be excluded. So to be included, an hkl needs to have values between the min and max @@ -856,43 +811,42 @@ def exclude( if (dmin is not None) or (dmax is not None): d = np.array(self.getPlaneSpacings()) - if (dmin is not None): + if dmin is not None: excl[d < dmin] = True - if (dmax is not None): + if dmax is not None: excl[d > dmax] = True if (tthmin is not None) or (tthmax is not None): tth = self.getTTh() - if (tthmin is not None): + if tthmin is not None: excl[tth < tthmin] = True - if (tthmax is not None): + if tthmax is not None: excl[tth > tthmax] = True if (sfacmin is not None) or (sfacmax is not None): sfac = self.structFact - sfac = sfac/sfac.max() - if (sfacmin is not None): + sfac = sfac / sfac.max() + if sfacmin is not None: excl[sfac < sfacmin] = True - if (sfacmax is not None): + if sfacmax is not None: excl[sfac > sfacmax] = True if (pintmin is not None) or (pintmax is not None): pint = self.powder_intensity - pint = pint/pint.max() - if (pintmin is not None): + pint = pint / pint.max() + if pintmin is not None: excl[pint < pintmin] = True - if (pintmax is not None): + if pintmax is not None: excl[pint > pintmax] = True self.exclusions = excl - def get_lparms(self): - return self.__lparms - - def __parseLParms(self, lparms): + def _parseLParms( + self, lparms: List[Union[valunits.valWUnit, float]] + ) -> List[float]: lparmsDUnit = [] for lparmThis in lparms: - if hasattr(lparmThis, 'getVal'): + if isinstance(lparmThis, valunits.valWUnit): if lparmThis.isLength(): lparmsDUnit.append(lparmThis.getVal(dUnit)) elif lparmThis.isAngle(): @@ -901,70 +855,104 @@ def __parseLParms(self, lparms): lparmsDUnit.append(lparmThis.getVal('degrees')) else: raise RuntimeError( - 'do not know what to do with ' - + str(lparmThis) + f'Do not know what to do with {lparmThis}' ) else: lparmsDUnit.append(lparmThis) return lparmsDUnit - def set_lparms(self, lparms): - self.__lparms = self.__parseLParms(lparms) - self.__calc() - return - lparms = property(get_lparms, set_lparms, None) + @property + def lparms(self) -> List[float]: + """ + Lattice parameters of the plane data. + + Can be set as a List[float | valWUnit], but will be converted to + List[float]. + """ + return self._lparms - def get_strainMag(self): - return self.__strainMag + @lparms.setter + def lparms(self, lparms: List[Union[valunits.valWUnit, float]]) -> None: + self._lparms = self._parseLParms(lparms) + self._calc() + + @property + def strainMag(self) -> Optional[float]: + """ + Strain magnitude of the plane data. + + float or None + """ + return self._strainMag - def set_strainMag(self, strainMag): - self.__strainMag = strainMag + @strainMag.setter + def strainMag(self, strain_mag: float) -> None: + self._strainMag = strain_mag self.tThWidth = None - self.__calc() - return - strainMag = property(get_strainMag, set_strainMag, None) + self._calc() - def get_wavelength(self): - return self.__wavelength + @property + def wavelength(self) -> float: + """ + Wavelength of the plane data. + + Set as float or valWUnit. - def set_wavelength(self, wavelength): + Read as float + """ + return self._wavelength + + @wavelength.setter + def wavelength(self, wavelength: Union[float, valunits.valWUnit]) -> None: wavelength = processWavelength(wavelength) - if np.isclose(self.__wavelength, wavelength): - # Do not re-compute if it is almost the same + # Do not re-compute if it is almost the same + if np.isclose(self._wavelength, wavelength): return - self.__wavelength = wavelength - self.__calc() - - wavelength = property(get_wavelength, set_wavelength, None) + self._wavelength = wavelength + self._calc() - def invalidate_structure_factor(self, unitcell): - # It can be expensive to compute the structure factor, so provide the - # option to just invalidate it, while providing a unit cell, so that - # it can be lazily computed from the unit cell. - self.__structFact = None + def invalidate_structure_factor(self, ucell: unitcell) -> None: + """ + It can be expensive to compute the structure factor + This method just invalidates it, providing a unit cell, + so that it can be lazily computed from the unit cell. + + Parameters: + ----------- + unitcell : unitcell + The unit cell to be used to compute the structure factor + """ + self._structFact = None self._hedm_intensity = None self._powder_intensity = None - self.__unitcell = unitcell + self._unitcell = ucell def _compute_sf_if_needed(self): any_invalid = ( - self.__structFact is None or - self._hedm_intensity is None or - self._powder_intensity is None + self._structFact is None + or self._hedm_intensity is None + or self._powder_intensity is None ) - if any_invalid and self.__unitcell is not None: + if any_invalid and self._unitcell is not None: # Compute the structure factor first. # This can be expensive to do, so we lazily compute it when needed. hkls = self.getHKLs(allHKLs=True) - self.set_structFact(self.__unitcell.CalcXRSF(hkls)) + self.structFact = self._unitcell.CalcXRSF(hkls) + + @property + def structFact(self) -> np.ndarray: + """ + Structure factors for each hkl. - def get_structFact(self): + np.ndarray + """ self._compute_sf_if_needed() - return self.__structFact[~self.exclusions] + return self._structFact[~self.exclusions] - def set_structFact(self, structFact): - self.__structFact = structFact + @structFact.setter + def structFact(self, structFact: np.ndarray) -> None: + self._structFact = structFact multiplicity = self.getMultiplicity(allHKLs=True) tth = self.getTTh(allHKLs=True) @@ -981,33 +969,70 @@ def set_structFact(self, structFact): self._hedm_intensity = hedm_intensity self._powder_intensity = powderI - structFact = property(get_structFact, set_structFact, None) - @property - def powder_intensity(self): + def powder_intensity(self) -> np.ndarray: + """ + Powder intensity for each hkl. + """ self._compute_sf_if_needed() return self._powder_intensity[~self.exclusions] @property - def hedm_intensity(self): + def hedm_intensity(self) -> np.ndarray: + """ + HEDM (high energy x-ray diffraction microscopy) intensity for each hkl. + """ self._compute_sf_if_needed() return self._hedm_intensity[~self.exclusions] @staticmethod - def makePlaneData(hkls, lparms, qsym, symmGroup, strainMag, wavelength): + def makePlaneData( + hkls: np.ndarray, + lparms: np.ndarray, + qsym: np.ndarray, + symmGroup, + strainMag, + wavelength, + ) -> Tuple[ + Dict[str, np.ndarray], Dict[str, Union[np.ndarray, float]], List[Dict] + ]: """ - hkls : need to work with crystallography.latticePlanes - lparms : need to work with crystallography.latticePlanes - laueGroup : see symmetry module - wavelength : wavelength - strainMag : swag of strian magnitudes + Generate lattice plane data from inputs. + + Parameters: + ----------- + hkls: np.ndarray + Miller indices, as in crystallography.latticePlanes + lparms: np.ndarray + Lattice parameters, as in crystallography.latticePlanes + qsym: np.ndarray + (4, n) containing quaternions of symmetry + symmGroup: str + Tag for the symmetry (Laue) group of the lattice. Can generate from + ltypeOfLaueGroup + strainMag: float + Swag of strain magnitudes + wavelength: float + Wavelength + + Returns: + ------- + dict: + Dictionary containing lattice plane data + dict: + Dictionary containing lattice vector operators + list: + List of dictionaries, each containing the data for one hkl """ tempSetOutputDegrees(False) - latPlaneData = latticePlanes(hkls, lparms, - ltype=symmGroup, - strainMag=strainMag, - wavelength=wavelength) + latPlaneData = latticePlanes( + hkls, + lparms, + ltype=symmGroup, + strainMag=strainMag, + wavelength=wavelength, + ) latVecOps = latticeVectors(lparms, symmGroup) @@ -1027,113 +1052,128 @@ def makePlaneData(hkls, lparms, qsym, symmGroup, strainMag, wavelength): np.dot(latVecOps['B'], hkls[:, iHKL].reshape(3, 1)), qsym, csFlag=True, - cullPM=False) + cullPM=False, + ) # check for +/- in symmetry group latPlnNrmlsM = applySym( np.dot(latVecOps['B'], hkls[:, iHKL].reshape(3, 1)), qsym, csFlag=False, - cullPM=False) + cullPM=False, + ) csRefl = latPlnNrmls.shape[1] == latPlnNrmlsM.shape[1] # added this so that I retain the actual symmetric # integer hkls as well symHKLs = np.array( - np.round( - np.dot(latVecOps['F'].T, latPlnNrmls) - ), dtype='int' + np.round(np.dot(latVecOps['F'].T, latPlnNrmls)), dtype='int' ) hklDataList.append( - dict(hklID=iHKL, - hkl=hkls[:, iHKL], - tTheta=latPlaneData['tThetas'][iHKL], - dSpacings=latPlaneData['dspacings'][iHKL], - tThetaLo=latPlaneData['tThetasLo'][iHKL], - tThetaHi=latPlaneData['tThetasHi'][iHKL], - latPlnNrmls=unitVector(latPlnNrmls), - symHKLs=symHKLs, - centrosym=csRefl - ) + dict( + hklID=iHKL, + hkl=hkls[:, iHKL], + tTheta=latPlaneData['tThetas'][iHKL], + dSpacings=latPlaneData['dspacings'][iHKL], + tThetaLo=latPlaneData['tThetasLo'][iHKL], + tThetaHi=latPlaneData['tThetasHi'][iHKL], + latPlnNrmls=unitVector(latPlnNrmls), + symHKLs=symHKLs, + centrosym=csRefl, ) + ) revertOutputDegrees() return latPlaneData, latVecOps, hklDataList - def getLatticeType(self): - """This is the lattice type""" - return ltypeOfLaueGroup(self.__laueGroup) - - def getLaueGroup(self): - """This is the Schoenflies tag""" - return self.__laueGroup + @property + def laueGroup(self) -> str: + """ + This is the Schoenflies tag, describing symmetry group of the lattice. + Note that setting this with incompatible lattice parameters will + cause an error. If changing both, use set_laue_and_lparms. - def setLaueGroup(self, laueGroup): - self.__laueGroup = laueGroup - self.__calc() + str + """ + return self._laueGroup - laueGroup = property(getLaueGroup, setLaueGroup, None) + @laueGroup.setter + def laueGroup(self, laueGroup: str) -> None: + self._laueGroup = laueGroup + self._calc() - def set_laue_and_lparms(self, laueGroup, lparms): - """Set the Laue group and lattice parameters simultaneously + def set_laue_and_lparms( + self, laueGroup: str, lparms: List[Union[valunits.valWUnit, float]] + ) -> None: + """ + Set the Laue group and lattice parameters simultaneously When the Laue group changes, the lattice parameters may be - incompatible, and cause an error in self.__calc(). This function + incompatible, and cause an error in self._calc(). This function allows us to update both the Laue group and lattice parameters simultaneously to avoid this issue. + + Parameters: + ----------- + laueGroup : str + The symmetry (Laue) group to be set + lparms : List[valunits.valWUnit | float] + Lattice parameters to be set """ - self.__laueGroup = laueGroup - self.__lparms = self.__parseLParms(lparms) - self.__calc() + self._laueGroup = laueGroup + self._lparms = self._parseLParms(lparms) + self._calc() - def getQSym(self): - return self.__qsym # rotations.quatOfLaueGroup(self.__laueGroup) + @property + def q_sym(self) -> np.ndarray: + """ + Quaternions of symmetry for each hkl, generated from the Laue group + + np.ndarray((4, n)) + """ + return self._q_sym # rotations.quatOfLaueGroup(self._laueGroup) - def getPlaneSpacings(self): + def getPlaneSpacings(self) -> List[float]: """ - gets plane spacings + Plane spacings for each hkl. + + Returns: + ------- + List[float] + List of plane spacings for each hkl """ dspacings = [] for iHKLr, hklData in enumerate(self.hklDataList): - if not self.__thisHKL(iHKLr): + if not self._thisHKL(iHKLr): continue dspacings.append(hklData['dSpacings']) return dspacings - def getPlaneNormals(self): - """ - gets both +(hkl) and -(hkl) normals - """ - plnNrmls = [] - for iHKLr, hklData in enumerate(self.hklDataList): - if not self.__thisHKL(iHKLr): - continue - plnNrmls.append(hklData['latPlnNrmls']) - return plnNrmls - @property - def latVecOps(self): + def latVecOps(self) -> Dict[str, Union[np.ndarray, float]]: """ gets lattice vector operators as a new (deepcopy) + + Returns: + ------- + Dict[str, np.ndarray | float] + Dictionary containing lattice vector operators """ return copy.deepcopy(self._latVecOps) - def __thisHKL(self, iHKLr): - retval = True + def _thisHKL(self, iHKLr: int) -> bool: hklData = self.hklDataList[iHKLr] - if self.__exclusions is not None: - if self.__exclusions[self.tThSortInv[iHKLr]]: - retval = False - if self.__tThMax is not None: - # FIXME: check for nans here??? - if hklData['tTheta'] > self.__tThMax \ - or np.isnan(hklData['tTheta']): - retval = False - return retval - - def __getTThRange(self, iHKLr): + if self._exclusions is not None: + if self._exclusions[self.tThSortInv[iHKLr]]: + return False + if self._tThMax is not None: + if hklData['tTheta'] > self._tThMax or np.isnan(hklData['tTheta']): + return False + return True + + def _getTThRange(self, iHKLr: int) -> Tuple[float, float]: hklData = self.hklDataList[iHKLr] if self.tThWidth is not None: # tThHi-tThLo < self.tThWidth tTh = hklData['tTheta'] @@ -1144,40 +1184,58 @@ def __getTThRange(self, iHKLr): tThLo = hklData['tThetaLo'] return (tThLo, tThHi) - def getTThRanges(self, strainMag=None, lparms=None): + def getTThRanges(self, strainMag: Optional[float] = None) -> np.ndarray: """ - Return 2-theta ranges for included hkls + Get the 2-theta ranges for included hkls + + Parameters: + ----------- + strainMag : Optional[float] + Optional swag of strain magnitude - return array is n x 2 + Returns: + ------- + np.ndarray: + hstacked array of hstacked tThLo and tThHi for each hkl (n x 2) """ - if lparms is None: - tThRanges = [] - for iHKLr, hklData in enumerate(self.hklDataList): - if not self.__thisHKL(iHKLr): - continue - # tThRanges.append([hklData['tThetaLo'], hklData['tThetaHi']]) - if strainMag is None: - tThRanges.append(self.__getTThRange(iHKLr)) - else: - hklData = self.hklDataList[iHKLr] - d = hklData['dSpacings'] - tThLo = 2.0 * np.arcsin( - self.__wavelength / 2. / (d*(1.+strainMag)) - ) - tThHi = 2.0 * np.arcsin( - self.__wavelength / 2. / (d*(1.-strainMag)) - ) - tThRanges.append((tThLo, tThHi)) - else: - new = self.__class__(self.__hkls, self) - new.lparms = lparms - tThRanges = new.getTThRanges(strainMag=strainMag) + tThRanges = [] + for iHKLr, hklData in enumerate(self.hklDataList): + if not self._thisHKL(iHKLr): + continue + if strainMag is None: + tThRanges.append(self._getTThRange(iHKLr)) + else: + hklData = self.hklDataList[iHKLr] + d = hklData['dSpacings'] + tThLo = 2.0 * np.arcsin( + self._wavelength / 2.0 / (d * (1.0 + strainMag)) + ) + tThHi = 2.0 * np.arcsin( + self._wavelength / 2.0 / (d * (1.0 - strainMag)) + ) + tThRanges.append((tThLo, tThHi)) return np.array(tThRanges) - def getMergedRanges(self, cullDupl=False): + def getMergedRanges( + self, cullDupl: Optional[bool] = False + ) -> Tuple[List[List[int]], List[List[float]]]: """ - return indices and ranges for specified planeData, merging where + Return indices and ranges for specified planeData, merging where there is overlap based on the tThWidth and line positions + + Parameters: + ----------- + cullDupl : (optional) bool + If True, cull duplicate 2-theta values (within sqrt_epsf). Defaults + to False. + + Returns: + -------- + List[List[int]] + List of indices for each merged range + + List[List[float]] + List of merged ranges, (n x 2) """ tThs = self.getTTh() tThRanges = self.getTThRanges() @@ -1191,11 +1249,11 @@ def getMergedRanges(self, cullDupl=False): mergedRanges = [] hklsCur = [] tThLoIdx = 0 - tThHiCur = 0. + tThHiCur = 0.0 for iHKL, nonoverlapNext in enumerate(nonoverlapNexts): tThHi = tThRanges[iHKL, -1] if not nonoverlapNext: - if cullDupl and abs(tThs[iHKL] - tThs[iHKL+1]) < sqrt_epsf: + if cullDupl and abs(tThs[iHKL] - tThs[iHKL + 1]) < sqrt_epsf: continue else: hklsCur.append(iHKL) @@ -1209,71 +1267,55 @@ def getMergedRanges(self, cullDupl=False): hklsCur = [] return iHKLLists, mergedRanges - def makeNew(self): - new = self.__class__(None, self) - return new - - def getTTh(self, lparms=None, allHKLs=False): - if lparms is None: - tTh = [] - for iHKLr, hklData in enumerate(self.hklDataList): - if not allHKLs: - if not self.__thisHKL(iHKLr): - continue - tTh.append(hklData['tTheta']) - else: - tTh.append(hklData['tTheta']) - else: - new = self.makeNew() - new.lparms = lparms - tTh = new.getTTh() - return np.array(tTh) - - def getDD_tThs_lparms(self): - """ - derivatives of tThs with respect to lattice parameters; - have not yet done coding for analytic derivatives, just wimp out - and finite difference + def getTTh(self, allHKLs: Optional[bool] = False) -> np.ndarray: """ - pert = 1.0e-5 # assume they are all around unity - pertInv = 1.0/pert + Get the 2-theta values for each hkl. - lparmsRef = copy.deepcopy(self.__lparms) - tThRef = self.getTTh() - ddtTh = np.empty((len(tThRef), len(lparmsRef))) + Parameters: + ----------- + allHKLs : (optional) bool + If True, return all 2-theta values, even if they are excluded in + the current planeData. Default is False. - for iLparm in range(len(lparmsRef)): - self.__lparms = copy.deepcopy(lparmsRef) - self.__lparms[iLparm] += pert - self.__calc() - - iTTh = 0 - for iHKLr, hklData in enumerate(self.hklDataList): - if not self.__thisHKL(iHKLr): - continue - ddtTh[iTTh, iLparm] = \ - np.r_[hklData['tTheta'] - tThRef[iTTh]]*pertInv - iTTh += 1 + Returns: + ------- + np.ndarray + Array of 2-theta values for each hkl + """ + tTh = [] + for iHKLr, hklData in enumerate(self.hklDataList): + if not allHKLs and not self._thisHKL(iHKLr): + continue + tTh.append(hklData['tTheta']) + return np.array(tTh) - 'restore' - self.__lparms = lparmsRef - self.__calc() + def getMultiplicity(self, allHKLs: Optional[bool] = False) -> np.ndarray: + """ + Get the multiplicity for each hkl (number of symHKLs). - return ddtTh + Paramters: + ---------- + allHKLs : (optional) bool + If True, return all multiplicities, even if they are excluded in + the current planeData. Defaults to false. - def getMultiplicity(self, allHKLs=False): + Returns + ------- + np.ndarray + Array of multiplicities for each hkl + """ # ... JVB: is this incorrect? multip = [] for iHKLr, hklData in enumerate(self.hklDataList): - if not allHKLs: - if not self.__thisHKL(iHKLr): - continue - multip.append(hklData['symHKLs'].shape[1]) - else: + if allHKLs or self._thisHKL(iHKLr): multip.append(hklData['symHKLs'].shape[1]) return np.array(multip) - def getHKLID(self, hkl, master=False): + def getHKLID( + self, + hkl: Union[int, Tuple[int, int, int], np.ndarray], + master: Optional[bool] = False, + ) -> Union[List[int], int]: """ Return the unique ID of a list of hkls. @@ -1293,7 +1335,7 @@ def getHKLID(self, hkl, master=False): Returns ------- - retval : list + hkl_ids : list The list of requested hklID values associate with the input. Notes @@ -1305,23 +1347,24 @@ def getHKLID(self, hkl, master=False): 2020-05-21 (JVB) -- modified to handle all symmetric equavlent reprs. """ if hasattr(hkl, '__setitem__'): # tuple does not have __setitem__ - if hasattr(hkl, 'shape'): + if isinstance(hkl, np.ndarray): # if is ndarray, assume is 3xN - # retval = list(map(self.__getHKLID, hkl.T)) - retval = [self.__getHKLID(x, master=master) for x in hkl.T] + return [self._getHKLID(x, master=master) for x in hkl.T] else: - # retval = list(map(self.__getHKLID, hkl)) - retval = [self.__getHKLID(x, master=master) for x in hkl] + return [self._getHKLID(x, master=master) for x in hkl] else: - retval = self.__getHKLID(hkl, master=master) - return retval + return self._getHKLID(hkl, master=master) - def __getHKLID(self, hkl, master=False): + def _getHKLID( + self, + hkl: Union[int, Tuple[int, int, int], np.ndarray], + master: Optional[bool] = False, + ) -> int: """ for hkl that is a tuple, return externally visible hkl index """ if isinstance(hkl, int): - retval = hkl + return hkl else: hklList = self.getSymHKLs() # !!! list, reduced by exclusions intl_hklIDs = np.asarray([i['hklID'] for i in self.hklDataList]) @@ -1332,14 +1375,13 @@ def __getHKLID(self, hkl, master=False): for thisHKL in symHKLs.T: dHKLInv[tuple(thisHKL)] = idx try: - retval = dHKLInv[tuple(hkl)] - except(KeyError): + return dHKLInv[tuple(hkl)] + except KeyError: raise RuntimeError( f"hkl '{tuple(hkl)}' is not present in this material!" ) - return retval - def getHKLs(self, *hkl_ids, **kwargs): + def getHKLs(self, *hkl_ids: int, **kwargs) -> Union[List[str], np.ndarray]: """ Returns the powder HKLs subject to specified options. @@ -1366,7 +1408,7 @@ def getHKLs(self, *hkl_ids, **kwargs): Returns ------- - retval : list | numpy.ndarray + hkls : list | numpy.ndarray Either a list of hkls as strings (if asStr=True) or a vstacked array of hkls. @@ -1396,10 +1438,10 @@ def getHKLs(self, *hkl_ids, **kwargs): if len(hkl_ids) == 0: for iHKLr, hklData in enumerate(self.hklDataList): if not opts['allHKLs']: - if not self.__thisHKL(iHKLr): + if not self._thisHKL(iHKLr): continue if opts['thisTTh'] is not None: - tThLo, tThHi = self.__getTThRange(iHKLr) + tThLo, tThHi = self._getTThRange(iHKLr) if opts['thisTTh'] < tThHi and opts['thisTTh'] > tThLo: hkls.append(hklData['hkl']) else: @@ -1415,7 +1457,7 @@ def getHKLs(self, *hkl_ids, **kwargs): # find ordinal index of current hklID try: idx[i] = int(np.where(all_hkl_ids == hkl_id)[0]) - except(TypeError): + except TypeError: raise RuntimeError( f"Requested hklID '{hkl_id}'is invalid!" ) @@ -1427,88 +1469,104 @@ def getHKLs(self, *hkl_ids, **kwargs): # handle output kwarg if opts['asStr']: - retval = list(map(hklToStr, np.array(hkls))) + return list(map(hklToStr, np.array(hkls))) else: - retval = np.array(hkls) - return retval - - def getSymHKLs(self, asStr=False, withID=False, indices=None): + return np.array(hkls) + + def getSymHKLs( + self, + asStr: Optional[bool] = False, + withID: Optional[bool] = False, + indices: Optional[List[int]] = None, + ) -> Union[List[List[str]], List[np.ndarray]]: """ - new function that returns all symmetric hkls + Return all symmetry HKLs. + + Parameters + ---------- + asStr : bool, optional + If True, return the symmetry HKLs as strings. The default is False. + withID : bool, optional + If True, return the symmetry HKLs with the hklID. The default is + False. Does nothing if asStr is True. + indices : list[inr], optional + Optional list of indices of hkls to include. + + Returns + ------- + sym_hkls : list list of strings, or list of numpy.ndarray + List of symmetry HKLs for each HKL, either as strings or as a + vstacked array. """ - retval = [] - iRetval = 0 + sym_hkls = [] + hkl_index = 0 if indices is not None: indB = np.zeros(self.nHKLs, dtype=bool) indB[np.array(indices)] = True else: indB = np.ones(self.nHKLs, dtype=bool) for iHKLr, hklData in enumerate(self.hklDataList): - if not self.__thisHKL(iHKLr): + if not self._thisHKL(iHKLr): continue - if indB[iRetval]: + if indB[hkl_index]: hkls = hklData['symHKLs'] if asStr: - retval.append(list(map(hklToStr, np.array(hkls).T))) + sym_hkls.append(list(map(hklToStr, np.array(hkls).T))) elif withID: - retval.append( + sym_hkls.append( np.vstack( - [np.tile(hklData['hklID'], (1, hkls.shape[1])), - hkls] + [ + np.tile(hklData['hklID'], (1, hkls.shape[1])), + hkls, + ] ) ) else: - retval.append(np.array(hkls)) - iRetval += 1 - return retval - - def getCentroSymHKLs(self): - retval = [] - for iHKLr, hklData in enumerate(self.hklDataList): - if not self.__thisHKL(iHKLr): - continue - retval.append(hklData['centrosym']) - return retval - - def makeTheseScatteringVectors(self, hklList, rMat, - bMat=None, wavelength=None, chiTilt=None): - iHKLList = np.atleast_1d(self.getHKLID(hklList)) - fHKLs = np.hstack(self.getSymHKLs(indices=iHKLList)) - if bMat is None: - bMat = self._latVecOps['B'] - if wavelength is None: - wavelength = self.__wavelength - retval = PlaneData.makeScatteringVectors( - fHKLs, rMat, bMat, wavelength, chiTilt=chiTilt - ) - return retval - - def makeAllScatteringVectors(self, rMat, - bMat=None, wavelength=None, chiTilt=None): - fHKLs = np.hstack(self.getSymHKLs()) - if bMat is None: - bMat = self._latVecOps['B'] - if wavelength is None: - wavelength = self.__wavelength - retval = PlaneData.makeScatteringVectors( - fHKLs, rMat, bMat, wavelength, chiTilt=chiTilt - ) - return retval + sym_hkls.append(np.array(hkls)) + hkl_index += 1 + return sym_hkls @staticmethod - def makeScatteringVectors(hkls, rMat_c, bMat, wavelength, chiTilt=None): + def makeScatteringVectors( + hkls: np.ndarray, + rMat_c: np.ndarray, + bMat: np.ndarray, + wavelength: float, + chiTilt: Optional[float] = None, + ) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: """ Static method for calculating g-vectors and scattering vector angles for specified hkls, subject to the bragg conditions specified by lattice vectors, orientation matrix, and wavelength + Parameters + ---------- + hkls : np.ndarray + (3, n) array of hkls. + rMat_c : np.ndarray + (3, 3) rotation matrix from the crystal to the sample frame. + bMat : np.ndarray, optional + (3, 3) COB from reciprocal lattice frame to the crystal frame. + wavelength : float + xray wavelength in Angstroms. + chiTilt : float, optional + 0 <= chiTilt <= 90 degrees, defaults to 0 + + Returns + ------- + gVec_s : np.ndarray + (3, n) array of g-vectors (reciprocal lattice) in the sample frame. + oangs0 : np.ndarray + (3, n) array containing the feasible (2-theta, eta, ome) triplets + for each input hkl (first solution) + oangs1 : np.ndarray + (3, n) array containing the feasible (2-theta, eta, ome) triplets + for each input hkl (second solution) + FIXME: must do testing on strained bMat """ # arg munging - if chiTilt is None: - chi = 0. - else: - chi = float(chiTilt) + chi = float(chiTilt) if chiTilt is not None else 0.0 rMat_c = rMat_c.squeeze() # these are the reciprocal lattice vectors in the SAMPLE FRAME @@ -1517,9 +1575,9 @@ def makeScatteringVectors(hkls, rMat_c, bMat, wavelength, chiTilt=None): # strained [a, b, c] in the CRYSTAL FRAME gVec_s = np.dot(rMat_c, np.dot(bMat, hkls)) - dim0, nRefl = gVec_s.shape - assert dim0 == 3, \ - "Looks like something is wrong with your lattice plane normals" + dim0 = gVec_s.shape[0] + if dim0 != 3: + raise ValueError(f'Number of lattice plane normal dims is {dim0}') # call model from transforms now oangs0, oangs1 = xfcapi.oscill_angles_of_hkls( @@ -1528,7 +1586,12 @@ def makeScatteringVectors(hkls, rMat_c, bMat, wavelength, chiTilt=None): return gVec_s, oangs0.T, oangs1.T - def __makeScatteringVectors(self, rMat, bMat=None, chiTilt=None): + def _makeScatteringVectors( + self, + rMat: np.ndarray, + bMat: Optional[np.ndarray] = None, + chiTilt: Optional[float] = None, + ) -> Tuple[List[np.ndarray], List[np.ndarray], List[np.ndarray]]: """ modeled after QFromU.m """ @@ -1540,11 +1603,14 @@ def __makeScatteringVectors(self, rMat, bMat=None, chiTilt=None): Qs_ang0 = [] Qs_ang1 = [] for iHKLr, hklData in enumerate(self.hklDataList): - if not self.__thisHKL(iHKLr): + if not self._thisHKL(iHKLr): continue thisQs, thisAng0, thisAng1 = PlaneData.makeScatteringVectors( - hklData['symHKLs'], rMat, bMat, - self.__wavelength, chiTilt=chiTilt + hklData['symHKLs'], + rMat, + bMat, + self._wavelength, + chiTilt=chiTilt, ) Qs_vec.append(thisQs) Qs_ang0.append(thisAng0) @@ -1552,70 +1618,30 @@ def __makeScatteringVectors(self, rMat, bMat=None, chiTilt=None): return Qs_vec, Qs_ang0, Qs_ang1 - def calcStructFactor(self, atominfo): - """ - Calculates unit cell structure factors as a function of hkl - - USAGE: - - FSquared = calcStructFactor(atominfo,hkls,B) - - INPUTS: - - 1) atominfo (m x 1 float ndarray) the first threee columns of the - matrix contain fractional atom positions [uvw] of atoms in the unit - cell. The last column contains the number of electrons for a given atom - - 2) hkls (3 x n float ndarray) is the array of Miller indices for - the planes of interest. The vectors are assumed to be - concatenated along the 1-axis (horizontal) - - 3) B (3 x 3 float ndarray) is a matrix of reciprocal lattice basis - vectors,where each column contains a reciprocal lattice basis vector - ({g}=[B]*{hkl}) - + # OLD DEPRECATED PLANE_DATA STUFF ==================================== + @deprecated(new_func="len(self.hkls.T)", removal_date="2025-08-01") + def getNHKLs(self): + return len(self.getHKLs()) - OUTPUTS: + @deprecated(new_func="self.exclusions", removal_date="2025-08-01") + def get_exclusions(self): + return self.exclusions - 1) FSquared (n x 1 float ndarray) array of structure factors, - one for each hkl passed into the function - """ - r = atominfo[:, 0:3] - elecNum = atominfo[:, 3] - hkls = self.hkls - B = self.latVecOps['B'] - sinThOverLamdaList, ffDataList = LoadFormFactorData() - FSquared = np.zeros(hkls.shape[1]) - - for jj in np.arange(0, hkls.shape[1]): - # ???: probably have other functions for this - # Calculate G for each hkl - # Calculate magnitude of G for each hkl - G = hkls[0, jj]*B[:, 0] + hkls[1, jj]*B[:, 1] + hkls[2, jj]*B[:, 2] - magG = np.sqrt(G[0]**2 + G[1]**2 + G[2]**2) - - # Begin calculating form factor - F = 0 - for ii in np.arange(0, r.shape[0]): - ff = RetrieveAtomicFormFactor( - elecNum[ii], magG, sinThOverLamdaList, ffDataList - ) - exparg = complex( - 0., - 2.*np.pi*(hkls[0, jj]*r[ii, 0] - + hkls[1, jj]*r[ii, 1] - + hkls[2, jj]*r[ii, 2]) - ) - F = F + ff*np.exp(exparg) + @deprecated(new_func="self.exclusions=...", removal_date="2025-08-01") + def set_exclusions(self, exclusions): + self.exclusions = exclusions - """ - F = sum_atoms(ff(Q)*e^(2*pi*i(hu+kv+lw))) - """ - FSquared[jj] = np.real(F*np.conj(F)) + @deprecated(new_func="rotations.ltypeOfLaueGroup(self.laueGroup)", + removal_date="2025-08-01") + def getLatticeType(self): + return ltypeOfLaueGroup(self.laueGroup) - return FSquared + @deprecated(new_func="self.q_sym", removal_date="2025-08-01") + def getQSym(self): + return self.q_sym +@deprecated(removal_date='2025-01-01') def getFriedelPair(tth0, eta0, *ome0, **kwargs): """ Get the diffractometer angular coordinates in degrees for @@ -1692,22 +1718,12 @@ def getFriedelPair(tth0, eta0, *ome0, **kwargs): dispFlag = False fableFlag = False chi = None - c1 = 1. - c2 = pi/180. - - # cast to arrays (in case they aren't) - if np.isscalar(eta0): - eta0 = [eta0] - - if np.isscalar(tth0): - tth0 = [tth0] - - if np.isscalar(ome0): - ome0 = [ome0] + c1 = 1.0 + c2 = pi / 180.0 - eta0 = np.asarray(eta0) - tth0 = np.asarray(tth0) - ome0 = np.asarray(ome0) + eta0 = np.atleast_1d(eta0) + tth0 = np.atleast_1d(tth0) + ome0 = np.atleast_1d(ome0) if eta0.ndim != 1: raise RuntimeError('azimuthal input must be 1-D') @@ -1719,25 +1735,24 @@ def getFriedelPair(tth0, eta0, *ome0, **kwargs): else: if len(tth0) != npts: if len(tth0) == 1: - tth0 = tth0*np.ones(npts) + tth0 *= np.ones(npts) elif npts == 1: npts = len(tth0) - eta0 = eta0*np.ones(npts) + eta0 *= np.ones(npts) else: raise RuntimeError( 'the azimuthal and Bragg angle inputs are inconsistent' ) if len(ome0) == 0: - ome0 = np.zeros(npts) # dummy ome0 + ome0 = np.zeros(npts) # dummy ome0 elif len(ome0) == 1 and npts > 1: - ome0 = ome0*np.ones(npts) + ome0 *= np.ones(npts) else: if len(ome0) != npts: raise RuntimeError( 'your oscialltion angle input is inconsistent; ' - + 'it has length %d while it should be %d' - % (len(ome0), npts) + + f'it has length {len(ome0)} while it should be {npts}' ) # keyword args processing @@ -1752,8 +1767,8 @@ def getFriedelPair(tth0, eta0, *ome0, **kwargs): fableFlag = True elif argkeys[i] == 'units': if kwargs[argkeys[i]] == 'radians': - c1 = 180./pi - c2 = 1. + c1 = 180.0 / pi + c2 = 1.0 elif argkeys[i] == 'chiTilt': if kwargs[argkeys[i]] is not None: chi = kwargs[argkeys[i]] @@ -1767,17 +1782,17 @@ def getFriedelPair(tth0, eta0, *ome0, **kwargs): # mapped eta input # - in DEGREES, thanks to c1 - eta0 = mapAngle(c1*eta0, [-180, 180], units='degrees') + eta0 = mapAngle(c1 * eta0, [-180, 180], units='degrees') if fableFlag: eta0 = 90 - eta0 # must put args into RADIANS # - eta0 is in DEGREES, # - the others are in whatever was entered, hence c2 - eta0 = d2r*eta0 - tht0 = c2*tth0/2 + eta0 = d2r * eta0 + tht0 = c2 * tth0 / 2 if chi is not None: - chi = c2*chi + chi = c2 * chi else: chi = 0 @@ -1789,15 +1804,15 @@ def getFriedelPair(tth0, eta0, *ome0, **kwargs): = sin(theta) - sin(chi)sin(eta)cos(theta) - Identity: a sin x + b cos x = sqrt(a**2 + b**2) sin (x + alfa) + Identity: a sin x + b cos x = sqrt(a**2 + b**2) sin (x + alpha) / | atan(b/a) for a > 0 - alfa < + alpha < | pi + atan(b/a) for a < 0 \ - => sin (x + alfa) = c / sqrt(a**2 + b**2) + => sin (x + alpha) = c / sqrt(a**2 + b**2) must use both branches for sin(x) = n: x = u (+ 2k*pi) | x = pi - u (+ 2k*pi) @@ -1809,23 +1824,20 @@ def getFriedelPair(tth0, eta0, *ome0, **kwargs): ctht = np.cos(tht0) stht = np.sin(tht0) - nchi = np.c_[0., cchi, schi].T + nchi = np.c_[0.0, cchi, schi].T - gHat0_l = -np.vstack([ceta * ctht, - seta * ctht, - stht]) + gHat0_l = -np.vstack([ceta * ctht, seta * ctht, stht]) - a = cchi*ceta*ctht - b = -cchi*stht - c = stht + schi*seta*ctht + a = cchi * ceta * ctht + b = -cchi * stht + c = stht + schi * seta * ctht # form solution - abMag = np.sqrt(a*a + b*b) - assert np.all(abMag > 0), \ - "Beam vector specification is infeasible!" + abMag = np.sqrt(a * a + b * b) + assert np.all(abMag > 0), "Beam vector specification is infeasible!" phaseAng = np.arctan2(b, a) rhs = c / abMag - rhs[abs(rhs) > 1.] = np.nan + rhs[abs(rhs) > 1.0] = np.nan rhsAng = np.arcsin(rhs) # write ome angle output arrays (NaNs persist here) @@ -1849,9 +1861,7 @@ def getFriedelPair(tth0, eta0, *ome0, **kwargs): tmp_eta = np.empty(numGood) tmp_gvec = gHat0_l[:, goodOnes] for i in range(numGood): - rchi = rotMatOfExpMap( - np.tile(ome_min[goodOnes][i], (3, 1)) * nchi - ) + rchi = rotMatOfExpMap(np.tile(ome_min[goodOnes][i], (3, 1)) * nchi) gHat_l = np.dot(rchi, tmp_gvec[:, i].reshape(3, 1)) tmp_eta[i] = np.arctan2(gHat_l[1], gHat_l[0]) eta_min[goodOnes] = tmp_eta @@ -1860,129 +1870,78 @@ def getFriedelPair(tth0, eta0, *ome0, **kwargs): # - ome1 is in RADIANS here # - convert and put into [-180, 180] ome1 = mapAngle( - mapAngle( - r2d*ome_min, [-180, 180], units='degrees' - ) + c1*ome0, [-180, 180], units='degrees' + mapAngle(r2d * ome_min, [-180, 180], units='degrees') + c1 * ome0, + [-180, 180], + units='degrees', ) # put eta1 in [-180, 180] - eta1 = mapAngle(r2d*eta_min, [-180, 180], units='degrees') + eta1 = mapAngle(r2d * eta_min, [-180, 180], units='degrees') if not outputDegrees: - ome1 = d2r * ome1 - eta1 = d2r * eta1 + ome1 *= d2r + eta1 *= d2r return ome1, eta1 -def getDparms(lp, lpTag, radians=True): +def getDparms( + lp: np.ndarray, lpTag: str, radians: Optional[bool] = True +) -> np.ndarray: """ Utility routine for getting dparms, that is the lattice parameters without symmetry -- 'triclinic' - """ - latVecOps = latticeVectors(lp, tag=lpTag, radians=radians) - return latVecOps['dparms'] - -def LoadFormFactorData(): - """ - Script to read in a csv file containing information relating the - magnitude of Q (sin(th)/lambda) to atomic form factor - - - Notes: - Atomic form factor data gathered from the International Tables of - Crystallography: - - P. J. Brown, A. G. Fox, E. N. Maslen, M. A. O'Keefec and B. T. M. Willis, - "Chapter 6.1. Intensity of diffracted intensities", International Tables - for Crystallography (2006). Vol. C, ch. 6.1, pp. 554-595 - """ - - dir1 = os.path.split(valunits.__file__) - dataloc = os.path.join(dir1[0], 'data', 'FormFactorVsQ.csv') - - data = np.zeros((62, 99), float) - - # FIXME: marked broken by DP - jj = 0 - with open(dataloc, 'rU') as csvfile: - datareader = csv.reader(csvfile, dialect=csv.excel) - for row in datareader: - ii = 0 - for val in row: - data[jj, ii] = float(val) - ii += 1 - jj += 1 - - sinThOverLamdaList = data[:, 0] - ffDataList = data[:, 1:] - - return sinThOverLamdaList, ffDataList - - -def RetrieveAtomicFormFactor(elecNum, magG, sinThOverLamdaList, ffDataList): - """ Interpolates between tabulated data to find the atomic form factor - for an atom with elecNum electrons for a given magnitude of Q - - USAGE: - - ff = RetrieveAtomicFormFactor(elecNum,magG,sinThOverLamdaList,ffDataList) - - INPUTS: - - 1) elecNum, (1 x 1 float) number of electrons for atom of interest - - 2) magG (1 x 1 float) magnitude of G - - 3) sinThOverLamdaList (n x 1 float ndarray) form factor data is tabulated - in terms of sin(theta)/lambda (A^-1). - - 3) ffDataList (n x m float ndarray) form factor data is tabulated in terms - of sin(theta)/lambda (A^-1). Each column corresponds to a different - number of electrons - - OUTPUTS: - - 1) ff (n x 1 float) atomic form factor for atom and hkl of interest - - NOTES: - Data should be calculated in terms of G at some point + Parameters + ---------- + lp : np.ndarray + Parsed lattice parameters + lpTag : str + Tag for the symmetry group of the lattice (from Laue group) + radians : bool, optional + Whether or not to use radians for angles, default is True + Returns + ------- + np.ndarray + The lattice parameters without symmetry. """ - sinThOverLambda = 0.5*magG - # lambda=2*d*sin(th) - # lambda=2*sin(th)/G - # 1/2*G=sin(th)/lambda - - ff = np.interp( - sinThOverLambda, sinThOverLamdaList, ffDataList[:, (elecNum - 1)] - ) - - return ff + latVecOps = latticeVectors(lp, tag=lpTag, radians=radians) + return latVecOps['dparms'] -def lorentz_factor(tth): +def lorentz_factor(tth: np.ndarray) -> np.ndarray: """ 05/26/2022 SS adding lorentz factor computation to the detector so that it can be compenstated for in the intensity correction - parameters: tth two theta of every pixel in radians + Parameters + ---------- + tth: np.ndarray + 2-theta of every pixel in radians + + Returns + ------- + np.ndarray + Lorentz factor for each pixel """ - theta = 0.5*tth + theta = 0.5 * tth cth = np.cos(theta) - sth2 = np.sin(theta)**2 + sth2 = np.sin(theta) ** 2 - L = 1./(4.0*cth*sth2) + return 1.0 / (4.0 * cth * sth2) - return L - -def polarization_factor(tth, unpolarized=True, eta=None, - f_hor=None, f_vert=None): +def polarization_factor( + tth: np.ndarray, + unpolarized: Optional[bool] = True, + eta: Optional[np.ndarray] = None, + f_hor: Optional[float] = None, + f_vert: Optional[float] = None, +) -> np.ndarray: """ 06/14/2021 SS adding lorentz polarization factor computation to the detector so that it can be compenstated for in the @@ -1998,13 +1957,16 @@ def polarization_factor(tth, unpolarized=True, eta=None, f_vert fraction of vertical polarization (~0 for XFELs) notice f_hor + f_vert = 1 + + FIXME, called without parameters like eta, f_hor, f_vert, but they default + to none in the current implementation, which will throw an error. """ - ctth2 = np.cos(tth)**2 + ctth2 = np.cos(tth) ** 2 if unpolarized: return (1 + ctth2) / 2 - seta2 = np.sin(eta)**2 - ceta2 = np.cos(eta)**2 - return f_hor*(seta2 + ceta2*ctth2) + f_vert*(ceta2 + seta2*ctth2) + seta2 = np.sin(eta) ** 2 + ceta2 = np.cos(eta) ** 2 + return f_hor * (seta2 + ceta2 * ctth2) + f_vert * (ceta2 + seta2 * ctth2) diff --git a/hexrd/material/material.py b/hexrd/material/material.py index 30b89e9b..b447c6fd 100644 --- a/hexrd/material/material.py +++ b/hexrd/material/material.py @@ -36,8 +36,11 @@ from hexrd.material.crystallography import PlaneData as PData from hexrd.material import symmetry, unitcell +from hexrd.material.symbols import two_origin_choice from hexrd.valunits import valWUnit -from hexrd.constants import ptable, ptableinverse, chargestate +from hexrd.constants import (ptable, + ptableinverse, + chargestate) from os import path from pathlib import Path @@ -77,6 +80,12 @@ def _key(x): return x.name +def get_default_sgsetting(sgnum): + if sgnum in two_origin_choice: + return 1 + else: + return 0 + # # ---------------------------------------------------CLASS: Material # @@ -148,7 +157,7 @@ def __init__( material_file=None, dmin=None, kev=None, - sgsetting=DFLT_SGSETTING, + sgsetting=None, ): """Constructor for Material @@ -165,8 +174,6 @@ def __init__( self.description = '' - self.sgsetting = sgsetting - if material_file: # >> @ date 08/20/2020 SS removing dependence on hklmax if isinstance(material_file, (Path, str)): @@ -180,6 +187,9 @@ def __init__( self._readCif(material_file) elif isinstance(material_file, h5py.Group) or form in h5_suffixes: self._readHDFxtal(fhdf=material_file, xtal=name) + if sgsetting is not None: + if sgsetting in [0, 1]: + self._sgsetting = sgsetting else: # default name self._name = Material.DFLT_XTAL @@ -189,8 +199,8 @@ def __init__( # self.description = '' # - self.sgnum = Material.DFLT_SGNUM self._sgsetting = Material.DFLT_SGSETTING + self.sgnum = Material.DFLT_SGNUM # self._atominfo = Material.DFLT_ATOMINFO # @@ -645,6 +655,7 @@ def _readCif(self, fcif=DFLT_NAME + '.cif'): lparms[i] = _degrees(lparms[i]) self._lparms = lparms + self._sgsetting = get_default_sgsetting(sgnum) self.sgnum = sgnum # fractional atomic site, occ and vibration amplitude @@ -795,7 +806,6 @@ def _readCif(self, fcif=DFLT_NAME + '.cif'): self._atomtype = np.asarray(atomtype).astype(np.int32) self._charge = charge - self._sgsetting = 0 self._dmin = Material.DFLT_DMIN self._beamEnergy = Material.DFLT_KEV @@ -864,6 +874,9 @@ def _readHDFxtal(self, fhdf=DFLT_NAME, xtal=DFLT_NAME): self._lparms = lparms # fill space group and lattice parameters + self._sgsetting = np.array( + gid.get('SpaceGroupSetting'), dtype=np.int32 + ).item() self.sgnum = sgnum # the U factors are related to B by the relation B = 8pi^2 U @@ -880,10 +893,6 @@ def _readHDFxtal(self, fhdf=DFLT_NAME, xtal=DFLT_NAME): self._charge = ['0'] * self._atomtype.shape[0] self._atom_ntype = self._atomtype.shape[0] - self._sgsetting = np.array( - gid.get('SpaceGroupSetting'), dtype=np.int32 - ).item() - if 'stiffness' in gid: # we're assuming the stiffness is in units of GPa self.stiffness = np.array(gid.get('stiffness')) @@ -1428,8 +1437,7 @@ def loadMaterialList(cfgFile): def load_materials_hdf5( - f, dmin=None, kev=None, sgsetting=Material.DFLT_SGSETTING -): + f, dmin=None, kev=None): """Load materials from an HDF5 file The file uses the HDF5 file format. @@ -1452,7 +1460,6 @@ def load_materials_hdf5( 'material_file': f, 'dmin': dmin, 'kev': kev, - 'sgsetting': sgsetting, } return {name: Material(name, **kwargs) for name in names} diff --git a/hexrd/material/mksupport.py b/hexrd/material/mksupport.py index 08327262..73f71870 100644 --- a/hexrd/material/mksupport.py +++ b/hexrd/material/mksupport.py @@ -1,6 +1,9 @@ -from hexrd.material.symbols import pstr_Elements, sitesym, \ - tworig, PrintPossibleSG, TRIG, pstr_spacegroup,\ - pstr_mkxtal +from hexrd.material.symbols import (pstr_Elements, + two_origin_choice, + PrintPossibleSG, + TRIG, + pstr_spacegroup, + pstr_mkxtal) import h5py import os import numpy as np @@ -8,6 +11,7 @@ import getpass from hexrd.material.unitcell import _StiffnessDict, _pgDict + def mk(filename, xtalname): # print some initial information for the user @@ -23,7 +27,7 @@ def mk(filename, xtalname): AtomInfo = GetAtomInfo() AtomInfo.update({'file': filename, 'xtalname': xtalname, - 'xtal_sys': xtal_sys, 'SG': space_group,\ + 'xtal_sys': xtal_sys, 'SG': space_group, 'SGsetting': iset}) Write2H5File(AtomInfo, lat_param) @@ -275,14 +279,15 @@ def GetSpaceGroup(xtal_sys, btrigonal, bhexset): def SpaceGroupSetting(sgnum): iset = 1 - if(sgnum in tworig): - idx = tworig.index(sgnum) + if(sgnum in two_origin_choice): + sitesym1 = two_origin_choice[sgnum][0] + sitesym2 = two_origin_choice[sgnum][1] print(' ---------------------------------------------') print(' This space group has two origin settings.') print(' The first setting has site symmetry : ' + - sitesym[2*idx - 2]) + sitesym1) print(' the second setting has site symmetry : ' + - sitesym[2*idx - 1]) + sitesym2) iset = input(' Which setting do you wish to use (1/2) : ') if(not iset.isdigit()): raise ValueError("Invalid integer value for atomic number.") @@ -294,7 +299,6 @@ def SpaceGroupSetting(sgnum): return iset-1 - def GetAtomInfo(): print(pstr_Elements) diff --git a/hexrd/material/symbols.py b/hexrd/material/symbols.py index ccc65ef3..637a5136 100644 --- a/hexrd/material/symbols.py +++ b/hexrd/material/symbols.py @@ -1,6 +1,5 @@ pstr_mkxtal = "\n\n This is a program to create a HDF5 file for storing crystallographic information.\n " -pstr_mkxtal = pstr_mkxtal + " This format is the same format as used in the EMsoft (electron microscoy) suite.\n " pstr_mkxtal = pstr_mkxtal + " The following inputs are required:\n " pstr_mkxtal = pstr_mkxtal + " Crystal System:\n" pstr_mkxtal = pstr_mkxtal + " 1. Cubic\n" @@ -24,88 +23,100 @@ pstr_mkxtal = pstr_mkxtal + " trigonal symmetry is the hexagonal setting. When you select\n" pstr_mkxtal = pstr_mkxtal + " crystal system 5 above, you will be prompted for the setting. \n" -pstr_spacegroup = [ " P 1 " ," P -1 ", \ -# MONOCLINIC SPACE GROUPS -" P 2 " ," P 21 " ," C 2 " ," P m ", \ -" P c " ," C m " ," C c " ," P 2/m ", \ -" P 21/m " ," C 2/m " ," P 2/c " ," P 21/c ", \ -" C 2/c ", \ -# ORTHORHOMBIC SPACE GROUPS -" P 2 2 2 " ," P 2 2 21 " ," P 21 21 2 " ," P 21 21 21", \ -" C 2 2 21 " ," C 2 2 2 " ," F 2 2 2 " ," I 2 2 2 ", \ -" I 21 21 21" ," P m m 2 " ," P m c 21 " ," P c c 2 ", \ -" P m a 2 " ," P c a 21 " ," P n c 2 " ," P m n 21 ", \ -" P b a 2 " ," P n a 21 " ," P n n 2 " ," C m m 2 ", \ -" C m c 21 " ," C c c 2 " ," A m m 2 " ," A b m 2 ", \ -" A m a 2 " ," A b a 2 " ," F m m 2 " ," F d d 2 ", \ -" I m m 2 " ," I b a 2 " ," I m a 2 " ," P m m m ", \ -" P n n n " ," P c c m " ," P b a n " ," P m m a ", \ -" P n n a " ," P m n a " ," P c c a " ," P b a m ", \ -" P c c n " ," P b c m " ," P n n m " ," P m m n ", \ -" P b c n " ," P b c a " ," P n m a " ," C m c m ", \ -" C m c a " ," C m m m " ," C c c m " ," C m m a ", \ -" C c c a " ," F m m m " ," F d d d " ," I m m m ", \ -" I b a m " ," I b c a " ," I m m a ", \ -# TETRAGONAL SPACE GROUPS -" P 4 " ," P 41 " ," P 42 " ," P 43 ", \ -" I 4 " ," I 41 " ," P -4 " ," I -4 ", \ -" P 4/m " ," P 42/m " ," P 4/n " ," P 42/n ", \ -" I 4/m " ," I 41/a " ," P 4 2 2 " ," P 4 21 2 ", \ -" P 41 2 2 " ," P 41 21 2 " ," P 42 2 2 " ," P 42 21 2 ", \ -" P 43 2 2 " ," P 43 21 2 " ," I 4 2 2 " ," I 41 2 2 ", \ -" P 4 m m " ," P 4 b m " ," P 42 c m " ," P 42 n m ", \ -" P 4 c c " ," P 4 n c " ," P 42 m c " ," P 42 b c ", \ -" I 4 m m " ," I 4 c m " ," I 41 m d " ," I 41 c d ", \ -" P -4 2 m " ," P -4 2 c " ," P -4 21 m " ," P -4 21 c ", \ -" P -4 m 2 " ," P -4 c 2 " ," P -4 b 2 " ," P -4 n 2 ", \ -" I -4 m 2 " ," I -4 c 2 " ," I -4 2 m " ," I -4 2 d ", \ -" P 4/m m m " ," P 4/m c c " ," P 4/n b m " ," P 4/n n c ", \ -" P 4/m b m " ," P 4/m n c " ," P 4/n m m " ," P 4/n c c ", \ -" P 42/m m c" ," P 42/m c m" ," P 42/n b c" ," P 42/n n m", \ -" P 42/m b c" ," P 42/m n m" ," P 42/n m c" ," P 42/n c m", \ -" I 4/m m m " ," I 4/m c m " ," I 41/a m d" ," I 41/a c d", \ -# RHOMBOHEDRAL SPACE GROUPS -" P 3 " ," P 31 " ," P 32 " ," R 3 ", \ -" P -3 " ," R -3 " ," P 3 1 2 " ," P 3 2 1 ", \ -" P 31 1 2 " ," P 31 2 1 " ," P 32 1 2 " ," P 32 2 1 ", \ -" R 3 2 " ," P 3 m 1 " ," P 3 1 m " ," P 3 c 1 ", \ -" P 3 1 c " ," R 3 m " ," R 3 c " ," P -3 1 m ", \ -" P -3 1 c " ," P -3 m 1 " ," P -3 c 1 " ," R -3 m ", \ -" R -3 c ", \ -# HEXAGONAL SPACE GROUPS -" P 6 " ," P 61 " ," P 65 " ," P 62 ", \ -" P 64 " ," P 63 " ," P -6 " ," P 6/m ", \ -" P 63/m " ," P 6 2 2 " ," P 61 2 2 " ," P 65 2 2 ", \ -" P 62 2 2 " ," P 64 2 2 " ," P 63 2 2 " ," P 6 m m ", \ -" P 6 c c " ," P 63 c m " ," P 63 m c " ," P -6 m 2 ", \ -" P -6 c 2 " ," P -6 2 m " ," P -6 2 c " ," P 6/m m m ", \ -" P 6/m c c " ," P 63/m c m" ," P 63/m m c", \ -#CUBIC SPACE GROUPS -" P 2 3 " ," F 2 3 " ," I 2 3 " ," P 21 3 ", \ -" I 21 3 " ," P m 3 " ," P n 3 " ," F m 3 ", \ -" F d 3 " ," I m 3 " ," P a 3 " ," I a 3 ", \ -" P 4 3 2 " ," P 42 3 2 " ," F 4 3 2 " ," F 41 3 2 ", \ -" I 4 3 2 " ," P 43 3 2 " ," P 41 3 2 " ," I 41 3 2 ", \ -" P -4 3 m " ," F -4 3 m " ," I -4 3 m " ," P -4 3 n ", \ -" F -4 3 c " ," I -4 3 d " ," P m 3 m " ," P n 3 n ", \ -" P m 3 n " ," P n 3 m " ," F m 3 m " ," F m 3 c ", \ -" F d 3 m " ," F d 3 c " ," I m 3 m " ," I a 3 d ", \ -# TRIGONAL GROUPS RHOMBOHEDRAL SETTING -" R 3 |146" ," R -3 |148" ," R 3 2 |155" ," R 3 m |160", \ -" R 3 c |161" ," R -3 m|166" ," R -3 c|167"] +pstr_spacegroup = [ + " P 1 ", " P -1 ", \ + # MONOCLINIC SPACE GROUPS + " P 2 ", " P 21 ", " C 2 ", " P m ", \ + " P c ", " C m ", " C c ", " P 2/m ", \ + " P 21/m ", " C 2/m ", " P 2/c ", " P 21/c ", \ + " C 2/c ", \ + # ORTHORHOMBIC SPACE GROUPS + " P 2 2 2 ", " P 2 2 21 ", " P 21 21 2 ", " P 21 21 21", \ + " C 2 2 21 ", " C 2 2 2 ", " F 2 2 2 ", " I 2 2 2 ", \ + " I 21 21 21", " P m m 2 ", " P m c 21 ", " P c c 2 ", \ + " P m a 2 ", " P c a 21 ", " P n c 2 ", " P m n 21 ", \ + " P b a 2 ", " P n a 21 ", " P n n 2 ", " C m m 2 ", \ + " C m c 21 ", " C c c 2 ", " A m m 2 ", " A b m 2 ", \ + " A m a 2 ", " A b a 2 ", " F m m 2 ", " F d d 2 ", \ + " I m m 2 ", " I b a 2 ", " I m a 2 ", " P m m m ", \ + " P n n n ", " P c c m ", " P b a n ", " P m m a ", \ + " P n n a ", " P m n a ", " P c c a ", " P b a m ", \ + " P c c n ", " P b c m ", " P n n m ", " P m m n ", \ + " P b c n ", " P b c a ", " P n m a ", " C m c m ", \ + " C m c a ", " C m m m ", " C c c m ", " C m m a ", \ + " C c c a ", " F m m m ", " F d d d ", " I m m m ", \ + " I b a m ", " I b c a ", " I m m a ", \ + # TETRAGONAL SPACE GROUPS + " P 4 ", " P 41 ", " P 42 ", " P 43 ", \ + " I 4 ", " I 41 ", " P -4 ", " I -4 ", \ + " P 4/m ", " P 42/m ", " P 4/n ", " P 42/n ", \ + " I 4/m ", " I 41/a ", " P 4 2 2 ", " P 4 21 2 ", \ + " P 41 2 2 ", " P 41 21 2 ", " P 42 2 2 ", " P 42 21 2 ", \ + " P 43 2 2 ", " P 43 21 2 ", " I 4 2 2 ", " I 41 2 2 ", \ + " P 4 m m ", " P 4 b m ", " P 42 c m ", " P 42 n m ", \ + " P 4 c c ", " P 4 n c ", " P 42 m c ", " P 42 b c ", \ + " I 4 m m ", " I 4 c m ", " I 41 m d ", " I 41 c d ", \ + " P -4 2 m ", " P -4 2 c ", " P -4 21 m ", " P -4 21 c ", \ + " P -4 m 2 ", " P -4 c 2 ", " P -4 b 2 ", " P -4 n 2 ", \ + " I -4 m 2 ", " I -4 c 2 ", " I -4 2 m ", " I -4 2 d ", \ + " P 4/m m m ", " P 4/m c c ", " P 4/n b m ", " P 4/n n c ", \ + " P 4/m b m ", " P 4/m n c ", " P 4/n m m ", " P 4/n c c ", \ + " P 42/m m c", " P 42/m c m", " P 42/n b c", " P 42/n n m", \ + " P 42/m b c", " P 42/m n m", " P 42/n m c", " P 42/n c m", \ + " I 4/m m m ", " I 4/m c m ", " I 41/a m d", " I 41/a c d", \ + # RHOMBOHEDRAL SPACE GROUPS + " P 3 ", " P 31 ", " P 32 ", " R 3 ", \ + " P -3 ", " R -3 ", " P 3 1 2 ", " P 3 2 1 ", \ + " P 31 1 2 ", " P 31 2 1 ", " P 32 1 2 ", " P 32 2 1 ", \ + " R 3 2 ", " P 3 m 1 ", " P 3 1 m ", " P 3 c 1 ", \ + " P 3 1 c ", " R 3 m ", " R 3 c ", " P -3 1 m ", \ + " P -3 1 c ", " P -3 m 1 ", " P -3 c 1 ", " R -3 m ", \ + " R -3 c ", \ + # HEXAGONAL SPACE GROUPS + " P 6 ", " P 61 ", " P 65 ", " P 62 ", \ + " P 64 ", " P 63 ", " P -6 ", " P 6/m ", \ + " P 63/m ", " P 6 2 2 ", " P 61 2 2 ", " P 65 2 2 ", \ + " P 62 2 2 ", " P 64 2 2 ", " P 63 2 2 ", " P 6 m m ", \ + " P 6 c c ", " P 63 c m ", " P 63 m c ", " P -6 m 2 ", \ + " P -6 c 2 ", " P -6 2 m ", " P -6 2 c ", " P 6/m m m ", \ + " P 6/m c c ", " P 63/m c m", " P 63/m m c", \ + # CUBIC SPACE GROUPS + " P 2 3 ", " F 2 3 ", " I 2 3 ", " P 21 3 ", \ + " I 21 3 ", " P m 3 ", " P n 3 ", " F m 3 ", \ + " F d 3 ", " I m 3 ", " P a 3 ", " I a 3 ", \ + " P 4 3 2 ", " P 42 3 2 ", " F 4 3 2 ", " F 41 3 2 ", \ + " I 4 3 2 ", " P 43 3 2 ", " P 41 3 2 ", " I 41 3 2 ", \ + " P -4 3 m ", " F -4 3 m ", " I -4 3 m ", " P -4 3 n ", \ + " F -4 3 c ", " I -4 3 d ", " P m 3 m ", " P n 3 n ", \ + " P m 3 n ", " P n 3 m ", " F m 3 m ", " F m 3 c ", \ + " F d 3 m ", " F d 3 c ", " I m 3 m ", " I a 3 d ", \ + # TRIGONAL GROUPS RHOMBOHEDRAL SETTING + " R 3 |146", " R -3 |148", " R 3 2 |155", " R 3 m |160", \ + " R 3 c |161", " R -3 m|166", " R -3 c|167"] -xtal_dict = {1:'cubic', 2:'tetragonal', 3:'orthorhombic', 4:'hexagonal', 5:'trigonal', 6:'monoclinic', 7:'triclinic'} -xtal_sys_dict = {'cubic':1, 'tetragonal':2, 'orthorhombic':3, 'hexagonal':4, 'trigonal':5, 'monoclinic':6, 'triclinic':7} -pstr_pointgroup = [' 1',' -1',' 2',' m',' 2/m',' 222', \ -' mm2',' mmm',' 4',' -4',' 4/m',' 422', \ -' 4mm',' -42m','4/mmm',' 3',' -3',' 32', \ -' 3m',' -3m',' 6',' -6',' 6/m',' 622', \ -' 6mm',' -6m2','6/mmm',' 23',' m3',' 432', \ -' -43m',' m-3m',' 532',' 822',' 1022',' 1222' ] +xtal_dict = {1: 'cubic', 2: 'tetragonal', 3: 'orthorhombic', + 4: 'hexagonal', 5: 'trigonal', 6: 'monoclinic', + 7: 'triclinic'} +xtal_sys_dict = {'cubic': 1, 'tetragonal': 2, 'orthorhombic': 3, + 'hexagonal': 4, 'trigonal': 5, 'monoclinic': 6, + 'triclinic': 7} + + + +pstr_pointgroup = [ + ' 1', ' -1', ' 2', ' m', ' 2/m', ' 222', + ' mm2', ' mmm', ' 4', ' -4', ' 4/m', ' 422', + ' 4mm', ' -42m', '4/mmm', ' 3', ' -3', ' 32', + ' 3m', ' -3m', ' 6', ' -6', ' 6/m', ' 622', + ' 6mm', ' -6m2', '6/mmm', ' 23', ' m3', ' 432', + ' -43m', ' m-3m', ' 532', ' 822', ' 1022', ' 1222'] + + + +TRIG = [146, 148, 155, 160, 161, 166, 167] -TRIG = [146,148,155,160,161,166,167] # symbols and Z for all elements pstr_Elements = ' ------------------------------------ Periodic Table of the Elements --------------------------------------' + "\n" \ @@ -120,52 +131,73 @@ '89:Ac 90:Th 91:Pa 92:U' + "\n" \ ' ----------------------------------------------------------------------------------------------------------' -sitesym = ['222 ',' -1 ','222/n ',' -1 ','mm2/n ',' -1 ', \ -'222 ',' -1 ','222 ',' -1 ','-4 ',' -1 ', \ -'-4 ',' -1 ','-4 ',' -1 ','422 ','2/m ', \ -'422/n ',' -1 ','-4m2 ','2/m ','-4/ncn ',' -1 ', \ -'-4121/c',' -1 ','-42m ','2/m ','-4m2/n ',' -1 ', \ -'-4cg ','2/m ','-4m2 ','2/m ','-4c21 ',' -1 ', \ -'23 ',' -3 ','23 ',' -3 ','432 ',' -3 ', \ -'-43m ','-3m ','-43m ','-3m ','23 ',' -3 '] +'''this dictionary contains the space group numbers +which have two origin choices the two values are the +site symmetries of the origin. There are 24 such space +groups''' -tworig = [48,50,59,68,70,85,86,88,125,126,129,130,133,134,137,138,\ -141,142,201,203,222,224,227,228] + +two_origin_choice = { + 48: ['222', '-1'], + 50: ['222/n', '-1'], + 59: ['mm2/n', '-1'], + 68: ['222', '-1'], + 70: ['222', '-1'], + 85: ['-4', '-1'], + 86: ['-4', '-1'], + 88: ['-4', '-1'], + 125: ['422', '2/m'], + 126: ['422/n', '-1'], + 129: ['-4m2', '2/m'], + 130: ['-4/ncn', '-1'], + 133: ['-4121/c', '-1'], + 134: ['-42m', '2/m'], + 137: ['-4m2/n', '-1'], + 138: ['-4cg', '2/m'], + 141: ['-4m2', '2/m'], + 142: ['-4c21', '-1'], + 201: ['23', '-3'], + 203: ['23', '-3'], + 222: ['432', '-3'], + 224: ['-43m', '-3m'], + 227: ['-43m', '-3m'], + 228: ['23', '-3'], +} def PrintPossibleSG(xtal_sys): - if(xtal_sys == 1): - sgmax = 230 - sgmin = 195 - elif(xtal_sys == 2): - sgmax = 142 - sgmin = 75 - elif(xtal_sys == 3): - sgmax = 74 - sgmin = 16 - elif(xtal_sys == 4): - sgmax = 194 - sgmin = 168 - elif(xtal_sys == 5): - sgmax = 167 - sgmin = 143 - elif(xtal_sys == 6): - sgmax = 15 - sgmin = 3 - elif(xtal_sys == 7): - sgmax = 2 - sgmin = 1 + if(xtal_sys == 1): + sgmax = 230 + sgmin = 195 + elif(xtal_sys == 2): + sgmax = 142 + sgmin = 75 + elif(xtal_sys == 3): + sgmax = 74 + sgmin = 16 + elif(xtal_sys == 4): + sgmax = 194 + sgmin = 168 + elif(xtal_sys == 5): + sgmax = 167 + sgmin = 143 + elif(xtal_sys == 6): + sgmax = 15 + sgmin = 3 + elif(xtal_sys == 7): + sgmax = 2 + sgmin = 1 - for i in range(sgmin,sgmax+1): - j = i - sgmin + 1 - pstr = str(i) + ":" + pstr_spacegroup[i-1] + "\t" - if(j % 4 == 0 or j == sgmax): - print(pstr) - else: - print(pstr, end='') - print("\n") + for i in range(sgmin, sgmax+1): + j = i - sgmin + 1 + pstr = str(i) + ":" + pstr_spacegroup[i-1] + "\t" + if(j % 4 == 0 or j == sgmax): + print(pstr) + else: + print(pstr, end='') + print("\n") - return sgmin, sgmax + return sgmin, sgmax # @@ -1237,6 +1269,7 @@ def PrintPossibleSG(xtal_sys): 230 I a -3 d """ + def _buildDict(hstr): """build the dictionaries from the notation string @@ -1247,7 +1280,7 @@ def _buildDict(hstr): for the rhombohedral lattices so that they use the hexagonal convention. """ - d = dict() + d = dict() di = dict() hs = hstr.split('\n') for l in hs: diff --git a/hexrd/material/symmetry.py b/hexrd/material/symmetry.py index ed8830fb..0d614f7f 100644 --- a/hexrd/material/symmetry.py +++ b/hexrd/material/symmetry.py @@ -31,9 +31,9 @@ import numpy as np from numba import njit -from numpy import array, sqrt, pi, \ - vstack, c_, dot, \ - argmax +from numpy import (array, sqrt, pi, + vstack, c_, dot, + argmax) # from hexrd.rotations import quatOfAngleAxis, quatProductMatrix, fixQuat from hexrd import rotations as rot @@ -41,8 +41,9 @@ from hexrd.utils.decorators import memoize # Imports in case others are importing from here -from hexrd.rotations import toFundamentalRegion, ltypeOfLaueGroup, \ - quatOfLaueGroup +from hexrd.rotations import (toFundamentalRegion, + ltypeOfLaueGroup, + quatOfLaueGroup) # ============================================================================= @@ -78,6 +79,7 @@ def GeneratorString(sgnum): return constants.SYM_GL[sg] + def MakeGenerators(genstr, setting): t = 'aOOO' @@ -98,7 +100,7 @@ def MakeGenerators(genstr, setting): if(n > 0): for i in range(n): istart = 2 + i * 4 - istop = 2 + (i+1) * 4 + istop = 2 + (i+1) * 4 t = genstr[istart:istop] @@ -113,32 +115,39 @@ def MakeGenerators(genstr, setting): if(genstr[istop] != '0'): if(setting != 0): t = genstr[istop+1:istop+4] - trans = np.array([constants.SYM_GENERATORS[t[0]],\ - constants.SYM_GENERATORS[t[1]],\ - constants.SYM_GENERATORS[t[2]] - ]) - for i in range(genmat.shape[0]): - genmat[i,0:3,3] -= trans + t = 'a' + t # get the translation without any rotation + sym = np.squeeze(SYM_fillgen(t, sgn=-1)) + sym2 = np.squeeze(SYM_fillgen(t)) + for i in range(1, genmat.shape[0]): + generator = np.dot(sym2, np.dot( + np.squeeze(genmat[i, :, :]), + sym)) + frac = np.modf(generator[0:3, 3])[0] + frac[frac < 0.] += 1. + frac[np.abs(frac) < 1E-5] = 0.0 + frac[np.abs(frac-1.0) < 1E-5] = 0.0 + generator[0:3, 3] = frac + genmat[i, :, :] = generator return genmat, centrosymmetric -def SYM_fillgen(t): - mat = np.zeros([4,4]) - mat[3,3] = 1. - mat[0:3,0:3] = constants.SYM_GENERATORS[t[0]] - mat[0:3,3] = np.array([constants.SYM_GENERATORS[t[1]],\ - constants.SYM_GENERATORS[t[2]],\ - constants.SYM_GENERATORS[t[3]] - ]) +def SYM_fillgen(t, sgn=1): + mat = np.zeros([4, 4]) + mat[3, 3] = 1. + + mat[0:3, 0:3] = constants.SYM_GENERATORS[t[0]] + mat[0:3, 3] = sgn*np.array([constants.SYM_GENERATORS[t[1]], + constants.SYM_GENERATORS[t[2]], + constants.SYM_GENERATORS[t[3]] + ]) - mat = np.broadcast_to(mat, [1,4,4]) + mat = np.broadcast_to(mat, [1, 4, 4]) return mat @memoize(maxsize=20) def GenerateSGSym(sgnum, setting=0): - ''' get the generators for a space group using the generator string @@ -165,22 +174,22 @@ def GenerateSGSym(sgnum, setting=0): k1 = 0 while k1 < nsym: - g1 = np.squeeze(SYM_SG[k1,:,:]) + g1 = np.squeeze(SYM_SG[k1, :, :]) k2 = k1 while k2 < nsym: - g2 = np.squeeze(SYM_SG[k2,:,:]) + g2 = np.squeeze(SYM_SG[k2, :, :]) gnew = np.dot(g1, g2) # only fractional parts - frac = np.modf(gnew[0:3,3])[0] + frac = np.modf(gnew[0:3, 3])[0] frac[frac < 0.] += 1. frac[np.abs(frac) < 1E-5] = 0.0 frac[np.abs(frac-1.0) < 1E-5] = 0.0 - gnew[0:3,3] = frac + gnew[0:3, 3] = frac if(isnew(gnew, SYM_SG)): - gnew = np.broadcast_to(gnew, [1,4,4]) + gnew = np.broadcast_to(gnew, [1, 4, 4]) SYM_SG = np.concatenate((SYM_SG, gnew)) nsym += 1 @@ -195,12 +204,12 @@ def GenerateSGSym(sgnum, setting=0): SYM_PG_d_laue = GeneratePGSym_Laue(SYM_PG_d) for s in SYM_PG_d: - if(np.allclose(-np.eye(3),s)): + if(np.allclose(-np.eye(3), s)): centrosymmetric = True - return SYM_SG, SYM_PG_d, SYM_PG_d_laue, centrosymmetric, symmorphic + def GeneratePGSym(SYM_SG): ''' calculate the direct space point group symmetries @@ -214,19 +223,20 @@ def GeneratePGSym(SYM_SG): nsgsym = SYM_SG.shape[0] # first fill the identity rotation - SYM_PG_d = SYM_SG[0,0:3,0:3] - SYM_PG_d = np.broadcast_to(SYM_PG_d,[1,3,3]) - - for i in range(1,nsgsym): - g = SYM_SG[i,:,:] - t = g[0:3,3] - g = g[0:3,0:3] - if(isnew(g,SYM_PG_d)): - g = np.broadcast_to(g,[1,3,3]) + SYM_PG_d = SYM_SG[0, 0:3, 0:3] + SYM_PG_d = np.broadcast_to(SYM_PG_d, [1, 3, 3]) + + for i in range(1, nsgsym): + g = SYM_SG[i, :, :] + t = g[0:3, 3] + g = g[0:3, 0:3] + if(isnew(g, SYM_PG_d)): + g = np.broadcast_to(g, [1, 3, 3]) SYM_PG_d = np.concatenate((SYM_PG_d, g)) return SYM_PG_d.astype(np.int32) + def GeneratePGSym_Laue(SYM_PG_d): ''' generate the laue group symmetry for the given set of @@ -240,7 +250,7 @@ def GeneratePGSym_Laue(SYM_PG_d): first check if the group already has the inversion symmetry ''' for s in SYM_PG_d: - if(np.allclose(s,-np.eye(3))): + if(np.allclose(s, -np.eye(3))): return SYM_PG_d ''' @@ -248,7 +258,7 @@ def GeneratePGSym_Laue(SYM_PG_d): add the inversion symmetry ''' SYM_PG_d_laue = SYM_PG_d - g = np.broadcast_to(-np.eye(3).astype(np.int32),[1,3,3]) + g = np.broadcast_to(-np.eye(3).astype(np.int32), [1, 3, 3]) SYM_PG_d_laue = np.concatenate((SYM_PG_d_laue, g)) ''' @@ -258,14 +268,14 @@ def GeneratePGSym_Laue(SYM_PG_d): nsym = SYM_PG_d_laue.shape[0] k1 = 0 while k1 < nsym: - g1 = np.squeeze(SYM_PG_d_laue[k1,:,:]) + g1 = np.squeeze(SYM_PG_d_laue[k1, :, :]) k2 = k1 while k2 < nsym: - g2 = np.squeeze(SYM_PG_d_laue[k2,:,:]) + g2 = np.squeeze(SYM_PG_d_laue[k2, :, :]) gnew = np.dot(g1, g2) if(isnew(gnew, SYM_PG_d_laue)): - gnew = np.broadcast_to(gnew, [1,3,3]) + gnew = np.broadcast_to(gnew, [1, 3, 3]) SYM_PG_d_laue = np.concatenate((SYM_PG_d_laue, gnew)) nsym += 1 @@ -287,6 +297,7 @@ def isnew(mat, sym_mats): return False return True + def latticeType(sgnum): if(sgnum <= 2): @@ -301,7 +312,7 @@ def latticeType(sgnum): return 'trigonal' elif(sgnum > 167 and sgnum <= 194): return 'hexagonal' - elif(sgnum > 194 and sgnum <=230): + elif(sgnum > 194 and sgnum <= 230): return 'cubic' else: raise RuntimeError('symmetry.latticeType: unknown space group number') @@ -315,14 +326,15 @@ def MakeGenerators_PGSYM(pggenstr): for any point group. this is needed for the coloring routines ''' ngen = int(pggenstr[0]) - SYM_GEN_PG = np.zeros([ngen,3,3]) + SYM_GEN_PG = np.zeros([ngen, 3, 3]) for i in range(ngen): s = pggenstr[i+1] - SYM_GEN_PG[i,:,:] = constants.SYM_GENERATORS[s] + SYM_GEN_PG[i, :, :] = constants.SYM_GENERATORS[s] return SYM_GEN_PG + def GeneratePGSYM(pgsym): ''' @AUTHOR Saransh Singh, Lawrence Livermore National Lab, saransh1@llnl.gov @@ -337,7 +349,6 @@ def GeneratePGSYM(pgsym): generate the powers of the group ''' - ''' now go through the group actions and see if its a new matrix if it is then add it to the group @@ -345,14 +356,14 @@ def GeneratePGSYM(pgsym): nsym = SYM_GEN_PG.shape[0] k1 = 0 while k1 < nsym: - g1 = np.squeeze(SYM_GEN_PG[k1,:,:]) + g1 = np.squeeze(SYM_GEN_PG[k1, :, :]) k2 = k1 while k2 < nsym: - g2 = np.squeeze(SYM_GEN_PG[k2,:,:]) + g2 = np.squeeze(SYM_GEN_PG[k2, :, :]) gnew = np.dot(g1, g2) if(isnew(gnew, SYM_GEN_PG)): - gnew = np.broadcast_to(gnew, [1,3,3]) + gnew = np.broadcast_to(gnew, [1, 3, 3]) SYM_GEN_PG = np.concatenate((SYM_GEN_PG, gnew)) nsym += 1 diff --git a/hexrd/material/unitcell.py b/hexrd/material/unitcell.py index 9f3b3f24..26f3a2fe 100644 --- a/hexrd/material/unitcell.py +++ b/hexrd/material/unitcell.py @@ -42,7 +42,7 @@ def _calcstar(v, sym, mat): if dist < 1E-3: isnew = False break - if(isnew): + if isnew: vp = np.atleast_2d(vp) vsym = np.vstack((vsym, vp)) @@ -158,7 +158,7 @@ def calcmatrices(self): [a*c*cb, b*c*ca, c**2]]) self._vol = np.sqrt(np.linalg.det(self.dmt)) - if(self.vol < 1e-5): + if self.vol < 1e-5: warnings.warn('unitcell volume is suspiciously small') ''' @@ -198,32 +198,34 @@ def calcmatrices(self): choices are 'd' (direct), 'r' (reciprocal) and 'c' (cartesian)''' def TransSpace(self, v_in, inspace, outspace): - if(inspace == 'd'): - if(outspace == 'r'): + if inspace == outspace: + return v_in + if inspace == 'd': + if outspace == 'r': v_out = np.dot(v_in, self.dmt) - elif(outspace == 'c'): + elif outspace == 'c': v_out = np.dot(self.dsm, v_in) else: raise ValueError( - 'inspace in ''d'' but outspace can''t be identified') + 'inspace in "d" but outspace can\'t be identified') - elif(inspace == 'r'): - if(outspace == 'd'): + elif inspace == 'r': + if outspace == 'd': v_out = np.dot(v_in, self.rmt) - elif(outspace == 'c'): + elif outspace == 'c': v_out = np.dot(self.rsm, v_in) else: raise ValueError( - 'inspace in ''r'' but outspace can''t be identified') + 'inspace in "r" but outspace can\'t be identified') - elif(inspace == 'c'): - if(outspace == 'r'): - v_out = np.dot(v_in, self.rsm) - elif(outspace == 'd'): + elif inspace == 'c': + if outspace == 'r': v_out = np.dot(v_in, self.dsm) + elif outspace == 'd': + v_out = np.dot(v_in, self.rsm) else: raise ValueError( - 'inspace in ''c'' but outspace can''t be identified') + 'inspace in "c" but outspace can\'t be identified') else: raise ValueError('incorrect inspace argument') @@ -234,11 +236,11 @@ def TransSpace(self, v_in, inspace, outspace): def CalcDot(self, u, v, space): - if(space == 'd'): + if space == 'd': dot = np.dot(u, np.dot(self.dmt, v)) - elif(space == 'r'): + elif space == 'r': dot = np.dot(u, np.dot(self.rmt, v)) - elif(space == 'c'): + elif space == 'c': dot = np.dot(u, v) else: raise ValueError('space is unidentified') @@ -247,13 +249,13 @@ def CalcDot(self, u, v, space): def CalcLength(self, u, space): - if(space == 'd'): + if space == 'd': mat = self.dmt # vlen = np.sqrt(np.dot(u, np.dot(self.dmt, u))) - elif(space == 'r'): + elif space == 'r': mat = self.rmt # vlen = np.sqrt(np.dot(u, np.dot(self.rmt, u))) - elif(space == 'c'): + elif space == 'c': mat = np.eye(3) # vlen = np.linalg.norm(u) else: @@ -297,7 +299,7 @@ def CalcAngle(self, u, v, space): def CalcCross(self, p, q, inspace, outspace, vol_divide=False): iv = 0 - if(vol_divide): + if vol_divide: vol = self.vol else: vol = 1.0 @@ -306,50 +308,50 @@ def CalcCross(self, p, q, inspace, outspace, vol_divide=False): p[2]*q[0]-p[0]*q[2], p[0]*q[1]-p[1]*q[0]]) - if(inspace == 'd'): + if inspace == 'd': ''' cross product vector is in reciprocal space and can be converted to direct or cartesian space ''' pxq *= vol - if(outspace == 'r'): + if outspace == 'r': pass - elif(outspace == 'd'): + elif outspace == 'd': pxq = self.TransSpace(pxq, 'r', 'd') - elif(outspace == 'c'): + elif outspace == 'c': pxq = self.TransSpace(pxq, 'r', 'c') else: raise ValueError( 'inspace is ''d'' but outspace is unidentified') - elif(inspace == 'r'): + elif inspace == 'r': ''' cross product vector is in direct space and can be converted to any other space ''' pxq /= vol - if(outspace == 'r'): + if outspace == 'r': pxq = self.TransSpace(pxq, 'd', 'r') - elif(outspace == 'd'): + elif outspace == 'd': pass - elif(outspace == 'c'): + elif outspace == 'c': pxq = self.TransSpace(pxq, 'd', 'c') else: raise ValueError( 'inspace is ''r'' but outspace is unidentified') - elif(inspace == 'c'): + elif inspace == 'c': ''' cross product is already in cartesian space so no volume factor is involved. can be converted to any other space too ''' - if(outspace == 'r'): + if outspace == 'r': pxq = self.TransSpace(pxq, 'c', 'r') - elif(outspace == 'd'): + elif outspace == 'd': pxq = self.TransSpace(pxq, 'c', 'd') - elif(outspace == 'c'): + elif outspace == 'c': pass else: raise ValueError( @@ -398,7 +400,7 @@ def GenerateCartesianPGSym(self): self.SYM_PG_c = np.array(self.SYM_PG_c) self.SYM_PG_c[np.abs(self.SYM_PG_c) < eps] = 0. - if(self._pointGroup == self._laueGroup): + if self._pointGroup == self._laueGroup: self.SYM_PG_c_laue = self.SYM_PG_c else: for sop in self.SYM_PG_d_laue: @@ -422,8 +424,7 @@ def GenerateCartesianPGSym(self): supergroup_laue = self._supergroup_laue sym_supergroup_laue = symmetry.GeneratePGSYM(supergroup_laue) - if((self.latticeType == 'monoclinic' or - self.latticeType == 'triclinic')): + if self.latticeType in ('monoclinic', 'triclinic'): ''' for monoclinic groups c2 and c2h, the supergroups are orthorhombic, so no need to convert from direct to @@ -462,7 +463,7 @@ def GenerateCartesianPGSym(self): not be rotated as they have the c* axis already aligned with the z-axis SS 12/10/2020 ''' - if(self.latticeType == 'monoclinic'): + if self.latticeType == 'monoclinic': om = np.array([[1., 0., 0.], [0., 0., 1.], [0., -1., 0.]]) @@ -480,7 +481,7 @@ def GenerateCartesianPGSym(self): triclinic group c1!! SS 12/10/2020 ''' - if(self._pointGroup == 'c1'): + if self._pointGroup == 'c1': om = np.array([[1., 0., 0.], [0., 0., 1.], [0., -1., 0.]]) for i, s in enumerate(self.SYM_PG_supergroup): @@ -536,7 +537,7 @@ def CalcOrbit(self, v, reduceToUC=True): break # if its new add this to the list - if(isnew): + if isnew: asym_pos = np.vstack((asym_pos, rr)) n += 1 @@ -549,21 +550,21 @@ def CalcStar(self, v, space, applyLaue=False): this function calculates the symmetrically equivalent hkls (or uvws) for the reciprocal (or direct) point group symmetry. ''' - if(space == 'd'): + if space == 'd': mat = self.dmt.astype(np.float64) - if(applyLaue): + if applyLaue: sym = self.SYM_PG_d_laue.astype(np.float64) else: sym = self.SYM_PG_d.astype(np.float64) - elif(space == 'r'): + elif space == 'r': mat = self.rmt.astype(np.float64) - if(applyLaue): + if applyLaue: sym = self.SYM_PG_r_laue.astype(np.float64) else: sym = self.SYM_PG_r.astype(np.float64) - elif(space == 'c'): + elif space == 'c': mat = np.eye(3) - if(applyLaue): + if applyLaue: sym = self.SYM_PG_c_laue.astype(np.float64) else: sym = self.SYM_PG_c.astype(np.float64) @@ -750,7 +751,11 @@ def InitializeInterpTable(self): f_anomalous_data = [] self.pe_cs = {} - data = importlib.resources.open_binary(hexrd.resources, 'Anomalous.h5') + data = ( + importlib.resources.files(hexrd.resources) + .joinpath('Anomalous.h5') + .open('rb') + ) with h5py.File(data, 'r') as fid: for i in range(0, self.atom_ntype): @@ -917,7 +922,7 @@ def ChooseSymmetric(self, hkllist, InversionSymmetry=True): laue = InversionSymmetry for i, g in enumerate(hkllist): - if(mask[i]): + if mask[i]: geqv = self.CalcStar(g, 'r', applyLaue=laue) @@ -995,11 +1000,11 @@ def getHKLs(self, dmin): for g in hkl_allowed: # ignore [0 0 0] as it is the direct beam - if(np.sum(np.abs(g)) != 0): + if np.sum(np.abs(g)) != 0: dspace = 1./self.CalcLength(g, 'r') - if(dspace >= dmin): + if dspace >= dmin: hkl_dsp.append(g) ''' @@ -1037,7 +1042,7 @@ def Required_C(self, C): return np.array([C[x] for x in _StiffnessDict[self._laueGroup][0]]) def MakeStiffnessMatrix(self, inp_Cvals): - if(len(inp_Cvals) != len(_StiffnessDict[self._laueGroup][0])): + if len(inp_Cvals) != len(_StiffnessDict[self._laueGroup][0]): x = len(_StiffnessDict[self._laueGroup][0]) msg = (f"number of constants entered is not correct." f" need a total of {x} independent constants.") @@ -1079,16 +1084,16 @@ def inside_spheretriangle(self, conn, dir3, hemisphere, switch): first get vertices of the triangles in the ''' vertex = self.sphere_sector.vertices[switch] - # if(switch == 'pg'): + # if switch == 'pg': # vertex = self.sphere_sector.vertices - # elif(switch == 'laue'): + # elif switch == 'laue': # vertex = self.sphere_sector.vertices_laue - # elif(switch == 'super'): + # elif switch == 'super': # vertex = self.sphere_sector.vertices_supergroup - # elif(switch == 'superlaue'): + # elif switch == 'superlaue': # vertex = self.sphere_sector.vertices_supergroup_laue A = np.atleast_2d(vertex[:, conn[0]]).T @@ -1107,29 +1112,29 @@ def inside_spheretriangle(self, conn, dir3, hemisphere, switch): determinant can be very small positive or negative number ''' - if(np.abs(d1) < eps): + if np.abs(d1) < eps: d1 = 0. - if(np.abs(d2) < eps): + if np.abs(d2) < eps: d2 = 0. - if(np.abs(d3) < eps): + if np.abs(d3) < eps: d3 = 0. ss = np.unique(np.sign([d1, d2, d3])) - if(hemisphere == 'upper'): - if(np.all(ss >= 0.)): + if hemisphere == 'upper': + if np.all(ss >= 0.): mask.append(True) else: mask.append(False) - elif(hemisphere == 'both'): - if(len(ss) == 1): + elif hemisphere == 'both': + if len(ss) == 1: mask.append(True) - elif(len(ss) == 2): - if(0 in ss): + elif len(ss) == 2: + if 0 in ss: mask.append(True) else: mask.append(False) - elif(len(ss) == 3): + elif len(ss) == 3: mask.append(False) mask = np.array(mask) @@ -1161,7 +1166,7 @@ def reduce_dirvector(self, dir3, switch='pg'): ''' idx = np.arange(dir3.shape[0], dtype=np.int32) dir3 = np.ascontiguousarray(np.atleast_2d(dir3)) - if(dir3.ndim != 2): + if dir3.ndim != 2: raise RuntimeError("reduce_dirvector: invalid shape of dir3 array") ''' @@ -1173,10 +1178,10 @@ def reduce_dirvector(self, dir3, switch='pg'): ''' eps = constants.sqrt_epsf - if(np.all(np.abs(np.linalg.norm(dir3, axis=1) - 1.0) < eps)): + if np.all(np.abs(np.linalg.norm(dir3, axis=1) - 1.0) < eps): dir3n = dir3 else: - if(np.all(np.linalg.norm(dir3) > eps)): + if np.all(np.linalg.norm(dir3) > eps): dir3n = dir3/np.tile(np.linalg.norm(dir3, axis=1), [3, 1]).T else: raise RuntimeError( @@ -1200,30 +1205,30 @@ def reduce_dirvector(self, dir3, switch='pg'): ntriangle = self.sphere_sector.ntriangle[switch] connectivity = self.sphere_sector.connectivity[switch] - if(switch == 'pg'): + if switch == 'pg': sym = self.SYM_PG_c - elif(switch == 'super'): + elif switch == 'super': sym = self.SYM_PG_supergroup - elif(switch == 'laue'): + elif switch == 'laue': sym = self.SYM_PG_c_laue - elif(switch == 'superlaue'): + elif switch == 'superlaue': sym = self.SYM_PG_supergroup_laue for sop in sym: - if(dir3_copy.size != 0): + if dir3_copy.size != 0: dir3_sym = np.dot(sop, dir3_copy.T).T mask = np.zeros(dir3_sym.shape[0]).astype(bool) - if(ntriangle == 0): - if(hemisphere == 'both'): + if ntriangle == 0: + if hemisphere == 'both': mask = np.ones(dir3_sym.shape[0], dtype=bool) - elif(hemisphere == 'upper'): + elif hemisphere == 'upper': mask = dir3_sym[:, 2] >= 0. else: for ii in range(ntriangle): @@ -1232,8 +1237,8 @@ def reduce_dirvector(self, dir3, switch='pg'): hemisphere, switch) mask = np.logical_or(mask, tmpmask) - if(np.sum(mask) > 0): - if(dir3_reduced.size != 0): + if np.sum(mask) > 0: + if dir3_reduced.size != 0: dir3_reduced = np.vstack( (dir3_reduced, dir3_sym[mask, :])) idx_red = np.hstack((idx_red, idx[mask])) @@ -1269,7 +1274,7 @@ class which correctly color the orientations for this crystal class. the replaceing barycenter to pi-theta) ''' - if(laueswitch == True): + if laueswitch: ''' this is the case where we color orientations based on the laue group of the crystal. this is always going to be the case with x-ray which @@ -1280,7 +1285,7 @@ class which correctly color the orientations for this crystal class. the dir3, switch='superlaue') switch = 'superlaue' - elif(laueswitch == False): + else: ''' follow the logic in the function description ''' @@ -1317,7 +1322,7 @@ def color_orientations(self, ''' first make sure that the rotation matric is size nx3x3 ''' - if(rmats.ndim == 2): + if rmats.ndim == 2: rmats = np.atleast_3d(rmats).T else: assert rmats.ndim == 3, "rotations matrices need to \ @@ -1358,7 +1363,7 @@ def convert_lp_to_valunits(self, lp): """ lp_valunit = [] for i in range(6): - if(i < 3): + if i < 3: lp_valunit.append( valWUnit('lp', 'length', lp[i], 'nm')) @@ -1415,7 +1420,7 @@ def lparms(self, lp): self.calcmatrices() self.init_max_g_index() self.CalcMaxGIndex() - if(hasattr(self, 'numat')): + if hasattr(self, 'numat'): self.CalcDensity() @property @@ -1541,7 +1546,7 @@ def U(self): def U(self, Uarr): self._U = Uarr self.aniU = False - if(Uarr.ndim > 1): + if Uarr.ndim > 1: self.aniU = True self.calcBetaij() @@ -1569,9 +1574,9 @@ def sgnum(self): @sgnum.setter def sgnum(self, val): - if(not(isinstance(val, int))): + if not(isinstance(val, int)): raise ValueError('space group should be integer') - if(not((val >= 1) and (val <= 230))): + if not((val >= 1) and (val <= 230)): raise ValueError('space group number should be between 1 and 230.') self._sym_sgnum = val diff --git a/hexrd/matrixutil.py b/hexrd/matrixutil.py index 64ce134e..d48c9f87 100644 --- a/hexrd/matrixutil.py +++ b/hexrd/matrixutil.py @@ -182,6 +182,7 @@ def vecMVToSymm(A, scale=True): convert from Mandel-Voigt vector to symmetric matrix representation (JVB) """ + A = np.atleast_1d(A).flatten() if scale: fac = sqr2 else: diff --git a/hexrd/xrdutil/utils.py b/hexrd/xrdutil/utils.py index 1e9961db..2245dbdb 100644 --- a/hexrd/xrdutil/utils.py +++ b/hexrd/xrdutil/utils.py @@ -262,109 +262,7 @@ def make_polar_net( return np.vstack(pts) -def validateAngleRanges( - angList: Union[np.ndarray, list[float]], - startAngs: Union[np.ndarray, list[float]], - stopAngs: Union[np.ndarray, list[float]], - ccw: bool = True, -) -> np.ndarray: - """ - Indetify angles that fall within specified ranges. - - A better way to go. find out if an angle is in the range - CCW or CW from start to stop - - There is, of course an ambigutiy if the start and stop angle are - the same; we treat them as implying 2*pi - """ - angList = np.atleast_1d(angList).flatten() - startAngs = np.atleast_1d(startAngs).flatten() - stopAngs = np.atleast_1d(stopAngs).flatten() - - if len(startAngs) != len(stopAngs): - raise ValueError("start and stop angles must have same length") - - # to avoid warnings in >=, <= later down, mark nans; - # need these to trick output to False in the case of nan input - nan_mask = np.isnan(angList) - - reflInRange = np.zeros(angList.shape, dtype=bool) - - # bin length for chunking - binLen = np.pi / 2.0 - - # in plane vectors defining wedges - x0 = np.vstack([np.cos(startAngs), np.sin(startAngs)]) - x1 = np.vstack([np.cos(stopAngs), np.sin(stopAngs)]) - - # dot products - dp = np.sum(x0 * x1, axis=0) - if np.any(dp >= 1.0 - sqrt_epsf) and len(startAngs) > 1: - # ambiguous case - raise RuntimeError( - "Improper usage; " - + "at least one of your ranges is already 360 degrees!" - ) - elif dp[0] >= 1.0 - sqrt_epsf and len(startAngs) == 1: - # trivial case! - reflInRange = np.ones(angList.shape, dtype=bool) - reflInRange[nan_mask] = False - else: - # solve for arc lengths - # ...note: no zeros should have made it here - a = x0[0, :] * x1[1, :] - x0[1, :] * x1[0, :] - b = x0[0, :] * x1[0, :] + x0[1, :] * x1[1, :] - phi = np.arctan2(b, a) - - arclen = 0.5 * np.pi - phi # these are clockwise - cw_phis = arclen < 0 - arclen[cw_phis] += 2 * np.pi # all positive (CW) now - if not ccw: - arclen = 2 * np.pi - arclen - - if sum(arclen) > 2 * np.pi: - raise RuntimeWarning( - "Specified angle ranges sum to > 360 degrees, " - + "which is suspect..." - ) - - # check that there are no more than dp = np.zeros(len(startAngs)) - for i in range(len(startAngs)): - # number or subranges using 'binLen' - numSubranges = int(np.ceil(arclen[i] / binLen)) - - # check remainder - binrem = np.remainder(arclen[i], binLen) - finalBinLen = binLen if binrem == 0 else binrem - - # if clockwise, negate bin length - if not ccw: - binLen = -binLen - finalBinLen = -finalBinLen - - # Create sub ranges on the fly to avoid ambiguity in dot product - # for wedges >= 180 degrees - subRanges = np.array( - [startAngs[i] + binLen * j for j in range(numSubranges)] - + [startAngs[i] + binLen * (numSubranges - 1) + finalBinLen] - ) - - for k in range(numSubranges): - zStart = _zproject(angList, subRanges[k]) - zStop = _zproject(angList, subRanges[k + 1]) - if ccw: - zStart[nan_mask] = 999.0 - zStop[nan_mask] = -999.0 - reflInRange = reflInRange | np.logical_and( - zStart <= 0, zStop >= 0 - ) - else: - zStart[nan_mask] = -999.0 - zStop[nan_mask] = 999.0 - reflInRange = reflInRange | np.logical_and( - zStart >= 0, zStop <= 0 - ) - return reflInRange +validateAngleRanges = xfcapi.validate_angle_ranges @deprecated(removal_date='2025-01-01') diff --git a/setup.py b/setup.py index 88a5b63e..ae91a962 100644 --- a/setup.py +++ b/setup.py @@ -34,37 +34,36 @@ # Add it to the requirements automatically for intel computers. install_reqs.append('tbb') -# This will determine which compiler is being used to build the C modules +# Determine which compiler is being used to build the C/C++ modules compiler_type = distutils.ccompiler.get_default_compiler() if compiler_type in ("unix", "mingw32"): - compiler_optimize_flags = ['-O3', '-ftree-vectorize'] + compiler_flags = ['-O3', '-ftree-vectorize', '-Wall', '-funroll-loops'] + if not sys.platform.startswith('win'): + compiler_flags.append('-fPIC') elif compiler_type == "msvc": - compiler_optimize_flags = ['/Ox', '/GL'] + compiler_flags = ['/Ox', '/GL', '/std:c++14'] else: - compiler_optimize_flags = [] - + compiler_flags = [] -# extension for convolution from astropy +# Extension for convolution from astropy def get_convolution_extensions(): c_convolve_pkgdir = Path('hexrd') / 'convolution' src_files = [str(c_convolve_pkgdir / 'src/convolve.c')] - extra_compile_args = ['-UNDEBUG'] - if not sys.platform.startswith('win'): - extra_compile_args.append('-fPIC') - extra_compile_args += compiler_optimize_flags + extra_compile_args = ['-UNDEBUG'] + compiler_flags # Add '-Rpass-missed=.*' to ``extra_compile_args`` when compiling with # clang to report missed optimizations - _convolve_ext = Extension(name='hexrd.convolution._convolve', - sources=src_files, - extra_compile_args=extra_compile_args, - include_dirs=[numpy.get_include()], - language='c') + _convolve_ext = Extension( + name='hexrd.convolution._convolve', + sources=src_files, + extra_compile_args=extra_compile_args, + include_dirs=[numpy.get_include()], + language='c' + ) return [_convolve_ext] - def get_include_path(library_name): env_var_hint = os.getenv(f"{library_name.upper()}_INCLUDE_DIR") if env_var_hint is not None and os.path.exists(env_var_hint): @@ -97,7 +96,6 @@ def get_include_path(library_name): # It should exist now return full_path - def get_pybind11_include_path(): # If we can import pybind11, use that include path try: @@ -110,16 +108,10 @@ def get_pybind11_include_path(): # Otherwise, we will download the source and include that return get_include_path('pybind11') - def get_cpp_extensions(): cpp_transform_pkgdir = Path('hexrd') / 'transforms/cpp_sublibrary' src_files = [str(cpp_transform_pkgdir / 'src/inverse_distortion.cpp')] - extra_compile_args = ['-O3', '-Wall', '-shared', '-std=c++14', - '-funroll-loops'] - if not sys.platform.startswith('win'): - extra_compile_args.append('-fPIC') - # Define include directories include_dirs = [ get_include_path('xsimd'), @@ -131,14 +123,13 @@ def get_cpp_extensions(): inverse_distortion_ext = Extension( name='hexrd.extensions.inverse_distortion', sources=src_files, - extra_compile_args=extra_compile_args, + extra_compile_args=compiler_flags+['-std=c++14'], include_dirs=include_dirs, language='c++', ) return [inverse_distortion_ext] - def get_old_xfcapi_extension_modules(): # for transforms srclist = ['transforms_CAPI.c', 'transforms_CFUNC.c'] @@ -147,23 +138,21 @@ def get_old_xfcapi_extension_modules(): 'hexrd.extensions._transforms_CAPI', sources=srclist, include_dirs=[np_include_dir], - extra_compile_args=compiler_optimize_flags, + extra_compile_args=compiler_flags, ) return [transforms_mod] - def get_new_xfcapi_extension_modules(): transforms_mod = Extension( 'hexrd.extensions._new_transforms_capi', sources=['hexrd/transforms/new_capi/module.c'], include_dirs=[np_include_dir], - extra_compile_args=compiler_optimize_flags, + extra_compile_args=compiler_flags, ) return [transforms_mod] - def get_extension_modules(): # Flatten the lists return [item for sublist in ( @@ -173,7 +162,6 @@ def get_extension_modules(): get_cpp_extensions(), ) for item in sublist] - ext_modules = get_extension_modules() # use entry_points, not scripts: diff --git a/tests/data/calibration_expected.npy b/tests/data/calibration_expected.npy new file mode 100644 index 00000000..587793d8 Binary files /dev/null and b/tests/data/calibration_expected.npy differ diff --git a/tests/planedata/test_exclusion.py b/tests/planedata/test_exclusion.py index bc8954b5..f85c0fc4 100644 --- a/tests/planedata/test_exclusion.py +++ b/tests/planedata/test_exclusion.py @@ -24,10 +24,7 @@ def test_exclusion(): assert pd.nHKLs == 2 assert pd.getNhklRef() == 3 - try: - # Exclude a range - pd.exclusions = [[0, 2]] - assert pd.nHKLs == 1 - assert pd.getNhklRef() == 3 - except NotImplementedError: - pass + # Exclude a range + pd.exclusions = [[0, 2]] + assert pd.nHKLs == 1 + assert pd.getNhklRef() == 3 diff --git a/tests/planedata/test_misc.py b/tests/planedata/test_misc.py index 9a168643..78a97703 100644 --- a/tests/planedata/test_misc.py +++ b/tests/planedata/test_misc.py @@ -1,3 +1,4 @@ +import os import numpy as np from hexrd.material.crystallography import PlaneData @@ -16,6 +17,8 @@ def test_misc(): 0.5, 0.01, ) + # Set ACK_DEPRECATED environment variable to true + os.environ['ACK_DEPRECATED'] = 'true' assert pd.laueGroup == 'C2h' assert pd.getLatticeType() == 'monoclinic' diff --git a/tests/planedata/test_with_data.py b/tests/planedata/test_with_data.py index 7b1e3ae6..0b4e9d53 100644 --- a/tests/planedata/test_with_data.py +++ b/tests/planedata/test_with_data.py @@ -1,3 +1,4 @@ +import os import numpy as np import pytest @@ -23,7 +24,8 @@ def materials(test_data_dir): for mat_name in material_names: # Load {test_data_dir}/materials/{mat_name}.cif mat = Material( - mat_name, str(test_data_dir) + "/materials/" + mat_name + ".cif" + mat_name, str(test_data_dir) + "/materials/" + mat_name + ".cif", + sgsetting=0 ) mats[mat_name] = mat.planeData return mats @@ -56,6 +58,7 @@ def flatten(xs): def test_plane_data_with_data(test_data_dir, materials): data = np.load(test_data_dir / 'plane_data_test.npy', allow_pickle=True) + os.environ['ACK_DEPRECATED'] = 'true' for obj in data: pd = materials[obj['mat_name']] assertEqualNumpyArrays(pd.hkls, obj['hkls']) @@ -84,6 +87,7 @@ def test_plane_data_with_data(test_data_dir, materials): # def test_with_data_maker(test_data_dir, materials): # data = [] +# os.environ['ACK_DEPRECATED'] = 'true' # for mat_name, pd in materials.items(): # obj = {"mat_name": mat_name} # obj['hkls'] = pd.hkls diff --git a/tests/test_calibration.py b/tests/test_calibration.py new file mode 100644 index 00000000..1939c96b --- /dev/null +++ b/tests/test_calibration.py @@ -0,0 +1,169 @@ +import copy + +import h5py +import numpy as np +import yaml + +import pytest + +from hexrd.material.material import load_materials_hdf5 +from hexrd.instrument.hedm_instrument import HEDMInstrument + +from hexrd.fitting.calibration import ( + InstrumentCalibrator, + LaueCalibrator, + PowderCalibrator, +) + + +@pytest.fixture +def calibration_dir(example_repo_path): + return example_repo_path / 'tardis' / 'simulated' + + +def test_calibration(calibration_dir, test_data_dir): + # Load the material + with h5py.File(calibration_dir / 'materials.h5', 'r') as rf: + materials = load_materials_hdf5(rf) + + # Load the picks + with open( + calibration_dir / 'uncalibrated_tardis.yml', 'r', encoding='utf-8' + ) as rf: + conf = yaml.safe_load(rf) + + instrument = HEDMInstrument(conf) + + # Load the images + img_npz = np.load(calibration_dir / 'tardis_images.npz') + img_dict = {k: img_npz[k] for k in img_npz} + + # Load picks + picks = np.load(calibration_dir / 'picks.npy', allow_pickle=True) + + euler_convention = {'axes_order': 'zxz', 'extrinsic': False} + + # Create the calibrators + calibrators = [] + for pick_data in picks: + if pick_data['type'] == 'powder': + kwargs = { + 'instr': instrument, + 'material': materials[pick_data['material']], + 'img_dict': img_dict, + 'default_refinements': pick_data['default_refinements'], + 'tth_distortion': pick_data['tth_distortion'], + 'calibration_picks': pick_data['picks'], + } + calibrators.append(PowderCalibrator(**kwargs)) + + elif pick_data['type'] == 'laue': + # gpflags = [i[1] for i in pick_data['refinements']] + kwargs = { + 'instr': instrument, + 'material': materials[pick_data['material']], + 'grain_params': pick_data['options']['crystal_params'], + 'default_refinements': pick_data['default_refinements'], + 'min_energy': pick_data['options']['min_energy'], + 'max_energy': pick_data['options']['max_energy'], + 'calibration_picks': pick_data['picks'], + 'euler_convention': euler_convention, + } + calibrators.append(LaueCalibrator(**kwargs)) + + calibrator = InstrumentCalibrator( + *calibrators, + engineering_constraints='TARDIS', + euler_convention=euler_convention, + ) + + tilt_angle_names = [ + [ + f'IMAGE_PLATE_{n}_euler_z', + f'IMAGE_PLATE_{n}_euler_xp', + f'IMAGE_PLATE_{n}_euler_zpp', + ] + for n in [2, 4] + ] + all_tilt_angle_names = tilt_angle_names[0] + tilt_angle_names[1] + + # The tilts are already ideal. Do not refine. + orig_params = copy.deepcopy(calibrator.params) + for name in all_tilt_angle_names: + calibrator.params[name].vary = False + + # The lattice parameter is actually perfect. Adjust it a little + # So that it can be corrected. + calibrator.params['diamond_a'].value = 3.58 + + x0 = calibrator.params.valuesdict() + result = calibrator.run_calibration({}) + x1 = result.params.valuesdict() + + # Parse the data + tvec_names = [ + [ + f'IMAGE_PLATE_{n}_tvec_x', + f'IMAGE_PLATE_{n}_tvec_y', + f'IMAGE_PLATE_{n}_tvec_z', + ] + for n in [2, 4] + ] + + tvecs = { + 'old': [np.array([x0[k] for k in vec_names]) for vec_names in tvec_names], + 'new': [np.array([x1[k] for k in vec_names]) for vec_names in tvec_names], + } + + grain_param_names = [f'LiF_grain_param_{n}' for n in range(12)] + grain_params = { + 'old': np.array([x0[name] for name in grain_param_names]), + 'new': np.array([x1[name] for name in grain_param_names]), + } + + diamond_a_vals = { + 'old': x0['diamond_a'], + 'new': x1['diamond_a'], + } + + expected = np.load( + test_data_dir / 'calibration_expected.npy', allow_pickle=True + ) + + assert_errors_are_better( + tvecs, grain_params, diamond_a_vals, expected.item() + ) + + +def assert_errors_are_better( + tvecs, grain_params, diamond_a_vals, expected +): + """ + Make sure error has decreased during fitting + """ + # What fraction of the old error we need to have (at worst) for the + # test to pass. For now, just make sure the error decreased. + ERROR_TOLERANCE = 1 + + tvec_error_2_old = np.linalg.norm(tvecs['old'][0] - expected['tvec_2']) + tvec_error_2_new = np.linalg.norm(tvecs['new'][0] - expected['tvec_2']) + tvec_error_4_old = np.linalg.norm(tvecs['old'][1] - expected['tvec_4']) + tvec_error_4_new = np.linalg.norm(tvecs['new'][1] - expected['tvec_4']) + + assert tvec_error_2_new < tvec_error_2_old * ERROR_TOLERANCE + assert tvec_error_4_new < tvec_error_4_old * ERROR_TOLERANCE + + grain_param_error_old = np.linalg.norm( + grain_params['old'] - expected['grain_params'] + ) + grain_param_error_new = np.linalg.norm( + grain_params['new'] - expected['grain_params'] + ) + assert grain_param_error_new < grain_param_error_old * ERROR_TOLERANCE + + diamond_a_error_old = np.abs(diamond_a_vals['old'] - expected['diamond_a']) + diamond_a_error_new = np.abs(diamond_a_vals['new'] - expected['diamond_a']) + + # The old diamond setting was actually perfect, but we let it refine + # The new diamond error should be less than 2% + assert diamond_a_error_new < diamond_a_error_old * ERROR_TOLERANCE diff --git a/tests/unitcell/test_vec_math.py b/tests/unitcell/test_vec_math.py new file mode 100644 index 00000000..794ad4da --- /dev/null +++ b/tests/unitcell/test_vec_math.py @@ -0,0 +1,103 @@ +from pytest import fixture +from hexrd.material import Material, unitcell +import numpy as np + + +@fixture +def cell() -> unitcell.unitcell: + return Material().unitcell + + +def get_space(): + return np.random.choice(['d', 'r', 'c']) + + +def test_calc_dot(cell: unitcell.unitcell): + np.random.seed(0) + for _ in range(100): + v1 = np.random.rand(3) * 10 - 5 + v2 = np.random.rand(3) * 10 - 5 + space = get_space() + v1c = cell.TransSpace(v1, space, 'c') + v2c = cell.TransSpace(v2, space, 'c') + assert np.allclose(np.dot(v1c, v2c), cell.CalcDot(v1, v2, space)) + + +def test_calc_length(cell: unitcell.unitcell): + np.random.seed(0) + for _ in range(100): + v1 = np.random.rand(3) * 10 - 5 + space = get_space() + v1c = cell.TransSpace(v1, space, 'c') + assert np.allclose(np.linalg.norm(v1c), cell.CalcLength(v1, space)) + + +def test_calc_angle(cell: unitcell.unitcell): + np.random.seed(0) + for _ in range(100): + v1 = np.random.rand(3) * 10 - 5 + v2 = np.random.rand(3) * 10 - 5 + space = get_space() + v1c = cell.TransSpace(v1, space, 'c') + v2c = cell.TransSpace(v2, space, 'c') + norms = np.linalg.norm(v1c) * np.linalg.norm(v2c) + assert np.allclose( + np.arccos(np.dot(v1c, v2c) / norms), cell.CalcAngle(v1, v2, space) + ) + + +def test_calc_cross(cell: unitcell.unitcell): + np.random.seed(0) + for _ in range(100): + v1 = np.random.rand(3) * 10 - 5 + v2 = np.random.rand(3) * 10 - 5 + inspace = get_space() + outspace = get_space() + v1c = cell.TransSpace(v1, inspace, 'c') + v2c = cell.TransSpace(v2, inspace, 'c') + expected = cell.TransSpace(np.cross(v1c, v2c), 'c', outspace) + assert np.allclose( + expected, cell.CalcCross(v1, v2, inspace, outspace, True) + ) + + +def test_norm_vec(cell: unitcell.unitcell): + np.random.seed(0) + for _ in range(100): + v1 = np.random.rand(3) * 10 - 5 + space = get_space() + norm_v1 = cell.NormVec(v1, space) + assert np.allclose(1, cell.CalcLength(norm_v1, space)) + # Make sure we don't change the direction + assert np.allclose( + v1 / np.linalg.norm(v1), norm_v1 / np.linalg.norm(norm_v1) + ) + + +def test_trans_space(cell: unitcell.unitcell): + np.random.seed(0) + for _ in range(100): + v1 = np.random.rand(3) * 10 - 5 + inspace = get_space() + outspace = get_space() + v1_out = cell.TransSpace(v1, inspace, outspace) + assert np.allclose(v1, cell.TransSpace(v1_out, outspace, inspace)) + + +def test_calc_star(cell: unitcell.unitcell): + """ + Just ensuring that the outspace doesn't matter + """ + np.random.seed(0) + for _ in range(100): + v1 = np.random.rand(3) * 10 - 5 + space = np.random.choice(['d', 'r']) + v1c = cell.TransSpace(v1, space, 'c') + assert np.allclose( + cell.CalcStar(v1, space, False), + cell.TransSpace(cell.CalcStar(v1c, 'c', False), 'c', space), + ) + assert np.allclose( + cell.CalcStar(v1, space, True), + cell.TransSpace(cell.CalcStar(v1c, 'c', True), 'c', space), + )