source: trunk/GSASIIlattice.py @ 3568

Last change on this file since 3568 was 3568, checked in by vondreele, 3 years ago

finish interface with k-SUBGROUPSMAG output - selection of space group & transformation to make magnetic phase

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 98.4 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIlattice: Unit cells*
4---------------------------
5
6Perform lattice-related computations
7
8Note that *g* is the reciprocal lattice tensor, and *G* is its inverse,
9:math:`G = g^{-1}`, where
10
11  .. math::
12
13   G = \\left( \\begin{matrix}
14   a^2 & a b\\cos\gamma & a c\\cos\\beta \\\\
15   a b\\cos\\gamma & b^2 & b c \cos\\alpha \\\\
16   a c\\cos\\beta &  b c \\cos\\alpha & c^2
17   \\end{matrix}\\right)
18
19The "*A* tensor" terms are defined as
20:math:`A = (\\begin{matrix} G_{11} & G_{22} & G_{33} & 2G_{12} & 2G_{13} & 2G_{23}\\end{matrix})` and *A* can be used in this fashion:
21:math:`d^* = \sqrt {A_1 h^2 + A_2 k^2 + A_3 l^2 + A_4 hk + A_5 hl + A_6 kl}`, where
22*d* is the d-spacing, and :math:`d^*` is the reciprocal lattice spacing,
23:math:`Q = 2 \\pi d^* = 2 \\pi / d`
24'''
25########### SVN repository information ###################
26# $Date: 2018-08-23 20:40:51 +0000 (Thu, 23 Aug 2018) $
27# $Author: vondreele $
28# $Revision: 3568 $
29# $URL: trunk/GSASIIlattice.py $
30# $Id: GSASIIlattice.py 3568 2018-08-23 20:40:51Z vondreele $
31########### SVN repository information ###################
32from __future__ import division, print_function
33import math
34import copy
35import sys
36import random as ran
37import numpy as np
38import numpy.linalg as nl
39import GSASIIpath
40import GSASIImath as G2mth
41import GSASIIspc as G2spc
42import GSASIIElem as G2elem
43GSASIIpath.SetVersionNumber("$Revision: 3568 $")
44# trig functions in degrees
45sind = lambda x: np.sin(x*np.pi/180.)
46asind = lambda x: 180.*np.arcsin(x)/np.pi
47tand = lambda x: np.tan(x*np.pi/180.)
48atand = lambda x: 180.*np.arctan(x)/np.pi
49atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
50cosd = lambda x: np.cos(x*np.pi/180.)
51acosd = lambda x: 180.*np.arccos(x)/np.pi
52rdsq2d = lambda x,p: round(1.0/np.sqrt(x),p)
53rpd = np.pi/180.
54RSQ2PI = 1./np.sqrt(2.*np.pi)
55SQ2 = np.sqrt(2.)
56RSQPI = 1./np.sqrt(np.pi)
57R2pisq = 1./(2.*np.pi**2)
58nxs = np.newaxis
59
60def sec2HMS(sec):
61    """Convert time in sec to H:M:S string
62   
63    :param sec: time in seconds
64    :return: H:M:S string (to nearest 100th second)
65   
66    """
67    H = int(sec//3600)
68    M = int(sec//60-H*60)
69    S = sec-3600*H-60*M
70    return '%d:%2d:%.2f'%(H,M,S)
71   
72def rotdMat(angle,axis=0):
73    """Prepare rotation matrix for angle in degrees about axis(=0,1,2)
74
75    :param angle: angle in degrees
76    :param axis:  axis (0,1,2 = x,y,z) about which for the rotation
77    :return: rotation matrix - 3x3 numpy array
78
79    """
80    if axis == 2:
81        return np.array([[cosd(angle),-sind(angle),0],[sind(angle),cosd(angle),0],[0,0,1]])
82    elif axis == 1:
83        return np.array([[cosd(angle),0,-sind(angle)],[0,1,0],[sind(angle),0,cosd(angle)]])
84    else:
85        return np.array([[1,0,0],[0,cosd(angle),-sind(angle)],[0,sind(angle),cosd(angle)]])
86       
87def rotdMat4(angle,axis=0):
88    """Prepare rotation matrix for angle in degrees about axis(=0,1,2) with scaling for OpenGL
89
90    :param angle: angle in degrees
91    :param axis:  axis (0,1,2 = x,y,z) about which for the rotation
92    :return: rotation matrix - 4x4 numpy array (last row/column for openGL scaling)
93
94    """
95    Mat = rotdMat(angle,axis)
96    return np.concatenate((np.concatenate((Mat,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0)
97   
98def fillgmat(cell):
99    """Compute lattice metric tensor from unit cell constants
100
101    :param cell: tuple with a,b,c,alpha, beta, gamma (degrees)
102    :return: 3x3 numpy array
103
104    """
105    a,b,c,alp,bet,gam = cell
106    g = np.array([
107        [a*a,  a*b*cosd(gam),  a*c*cosd(bet)],
108        [a*b*cosd(gam),  b*b,  b*c*cosd(alp)],
109        [a*c*cosd(bet) ,b*c*cosd(alp),   c*c]])
110    return g
111           
112def cell2Gmat(cell):
113    """Compute real and reciprocal lattice metric tensor from unit cell constants
114
115    :param cell: tuple with a,b,c,alpha, beta, gamma (degrees)
116    :return: reciprocal (G) & real (g) metric tensors (list of two numpy 3x3 arrays)
117
118    """
119    g = fillgmat(cell)
120    G = nl.inv(g)       
121    return G,g
122
123def A2Gmat(A,inverse=True):
124    """Fill real & reciprocal metric tensor (G) from A.
125
126    :param A: reciprocal metric tensor elements as [G11,G22,G33,2*G12,2*G13,2*G23]
127    :param bool inverse: if True return both G and g; else just G
128    :return: reciprocal (G) & real (g) metric tensors (list of two numpy 3x3 arrays)
129
130    """
131    G = np.array([
132        [A[0],  A[3]/2.,  A[4]/2.], 
133        [A[3]/2.,A[1],    A[5]/2.], 
134        [A[4]/2.,A[5]/2.,    A[2]]])
135    if inverse:
136        g = nl.inv(G)
137        return G,g
138    else:
139        return G
140
141def Gmat2A(G):
142    """Extract A from reciprocal metric tensor (G)
143
144    :param G: reciprocal maetric tensor (3x3 numpy array
145    :return: A = [G11,G22,G33,2*G12,2*G13,2*G23]
146
147    """
148    return [G[0][0],G[1][1],G[2][2],2.*G[0][1],2.*G[0][2],2.*G[1][2]]
149   
150def cell2A(cell):
151    """Obtain A = [G11,G22,G33,2*G12,2*G13,2*G23] from lattice parameters
152
153    :param cell: [a,b,c,alpha,beta,gamma] (degrees)
154    :return: G reciprocal metric tensor as 3x3 numpy array
155
156    """
157    G,g = cell2Gmat(cell)
158    return Gmat2A(G)
159
160def A2cell(A):
161    """Compute unit cell constants from A
162
163    :param A: [G11,G22,G33,2*G12,2*G13,2*G23] G - reciprocal metric tensor
164    :return: a,b,c,alpha, beta, gamma (degrees) - lattice parameters
165
166    """
167    G,g = A2Gmat(A)
168    return Gmat2cell(g)
169
170def Gmat2cell(g):
171    """Compute real/reciprocal lattice parameters from real/reciprocal metric tensor (g/G)
172    The math works the same either way.
173
174    :param g (or G): real (or reciprocal) metric tensor 3x3 array
175    :return: a,b,c,alpha, beta, gamma (degrees) (or a*,b*,c*,alpha*,beta*,gamma* degrees)
176
177    """
178    oldset = np.seterr('raise')
179    a = np.sqrt(max(0,g[0][0]))
180    b = np.sqrt(max(0,g[1][1]))
181    c = np.sqrt(max(0,g[2][2]))
182    alp = acosd(g[2][1]/(b*c))
183    bet = acosd(g[2][0]/(a*c))
184    gam = acosd(g[0][1]/(a*b))
185    np.seterr(**oldset)
186    return a,b,c,alp,bet,gam
187
188def invcell2Gmat(invcell):
189    """Compute real and reciprocal lattice metric tensor from reciprocal
190       unit cell constants
191       
192    :param invcell: [a*,b*,c*,alpha*, beta*, gamma*] (degrees)
193    :return: reciprocal (G) & real (g) metric tensors (list of two 3x3 arrays)
194
195    """
196    G = fillgmat(invcell)
197    g = nl.inv(G)
198    return G,g
199
200def cellDijFill(pfx,phfx,SGData,parmDict): 
201    '''Returns the filled-out reciprocal cell (A) terms
202    from the parameter dictionaries corrected for Dij.
203
204    :param str pfx: parameter prefix ("n::", where n is a phase number)
205    :param dict SGdata: a symmetry object
206    :param dict parmDict: a dictionary of parameters
207
208    :returns: A,sigA where each is a list of six terms with the A terms
209    '''
210    if SGData['SGLaue'] in ['-1',]:
211        A = [parmDict[pfx+'A0']+parmDict[phfx+'D11'],parmDict[pfx+'A1']+parmDict[phfx+'D22'],
212             parmDict[pfx+'A2']+parmDict[phfx+'D33'],
213             parmDict[pfx+'A3']+parmDict[phfx+'D12'],parmDict[pfx+'A4']+parmDict[phfx+'D13'],
214             parmDict[pfx+'A5']+parmDict[phfx+'D23']]
215    elif SGData['SGLaue'] in ['2/m',]:
216        if SGData['SGUniq'] == 'a':
217            A = [parmDict[pfx+'A0']+parmDict[phfx+'D11'],parmDict[pfx+'A1']+parmDict[phfx+'D22'],
218                 parmDict[pfx+'A2']+parmDict[phfx+'D33'],0,0,parmDict[pfx+'A5']+parmDict[phfx+'D23']]
219        elif SGData['SGUniq'] == 'b':
220            A = [parmDict[pfx+'A0']+parmDict[phfx+'D11'],parmDict[pfx+'A1']+parmDict[phfx+'D22'],
221                 parmDict[pfx+'A2']+parmDict[phfx+'D33'],0,parmDict[pfx+'A4']+parmDict[phfx+'D13'],0]
222        else:
223            A = [parmDict[pfx+'A0']+parmDict[phfx+'D11'],parmDict[pfx+'A1']+parmDict[phfx+'D22'],
224                 parmDict[pfx+'A2']+parmDict[phfx+'D33'],parmDict[pfx+'A3']+parmDict[phfx+'D12'],0,0]
225    elif SGData['SGLaue'] in ['mmm',]:
226        A = [parmDict[pfx+'A0']+parmDict[phfx+'D11'],parmDict[pfx+'A1']+parmDict[phfx+'D22'],
227             parmDict[pfx+'A2']+parmDict[phfx+'D33'],0,0,0]
228    elif SGData['SGLaue'] in ['4/m','4/mmm']:
229        A = [parmDict[pfx+'A0']+parmDict[phfx+'D11'],parmDict[pfx+'A0']+parmDict[phfx+'D11'],
230             parmDict[pfx+'A2']+parmDict[phfx+'D33'],0,0,0]
231    elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
232        A = [parmDict[pfx+'A0']+parmDict[phfx+'D11'],parmDict[pfx+'A0']+parmDict[phfx+'D11'],
233             parmDict[pfx+'A2']+parmDict[phfx+'D33'],parmDict[pfx+'A0']+parmDict[phfx+'D11'],0,0]
234    elif SGData['SGLaue'] in ['3R', '3mR']:
235        A = [parmDict[pfx+'A0']+parmDict[phfx+'D11'],parmDict[pfx+'A0']+parmDict[phfx+'D11'],
236            parmDict[pfx+'A0']+parmDict[phfx+'D11'],
237            parmDict[pfx+'A3']+parmDict[phfx+'D23'],parmDict[pfx+'A3']+parmDict[phfx+'D23'],
238            parmDict[pfx+'A3']+parmDict[phfx+'D23']]
239    elif SGData['SGLaue'] in ['m3m','m3']:
240        A = [parmDict[pfx+'A0']+parmDict[phfx+'D11'],parmDict[pfx+'A0']+parmDict[phfx+'D11'],
241             parmDict[pfx+'A0']+parmDict[phfx+'D11'],0,0,0]
242    return A
243   
244def prodMGMT(G,Mat):
245    '''Transform metric tensor by matrix
246   
247    :param G: array metric tensor
248    :param Mat: array transformation matrix
249    :return: array new metric tensor
250   
251    '''
252    return np.inner(Mat,np.inner(G,Mat).T)
253   
254def TransformCell(cell,Trans):
255    '''Transform lattice parameters by matrix
256   
257    :param cell: list a,b,c,alpha,beta,gamma,(volume)
258    :param Trans: array transformation matrix
259    :return: array transformed a,b,c,alpha,beta,gamma,volume
260   
261    '''
262    newCell = np.zeros(7)
263    g = cell2Gmat(cell)[1]
264    newg = prodMGMT(g,Trans)
265    newCell[:6] = Gmat2cell(newg)
266    newCell[6] = calc_V(cell2A(newCell[:6]))
267    return newCell
268   
269def TransformXYZ(XYZ,Trans,Vec):
270    return np.inner(XYZ,Trans)+Vec
271   
272def TransformU6(U6,Trans):
273    Uij = np.inner(Trans,np.inner(U6toUij(U6),Trans).T)/nl.det(Trans)
274    return UijtoU6(Uij)
275   
276def TransformPhase(oldPhase,newPhase,Trans,Uvec,Vvec,ifMag):
277    '''Transform atoms from oldPhase to newPhase
278    M' is inv(M)
279    does X' = (X-U)M'+V transformation for coordinates and U' = MUM/det(M)
280    for anisotropic thermal parameters
281   
282    :param oldPhase: dict G2 phase info for old phase
283    :param newPhase: dict G2 phase info for new phase; with new cell & space group
284            atoms are from oldPhase & will be transformed
285    :param Trans: lattice transformation matrix M
286    :param Uvec: array parent coordinates transformation vector U
287    :param Vvec: array child coordinate transformation vector V
288    :param ifMag: bool True if convert to magnetic phase;
289        if True all nonmagnetic atoms will be removed
290       
291    '''
292   
293    cx,ct,cs,cia = oldPhase['General']['AtomPtrs']
294    cm = 0
295    if oldPhase['General']['Type'] == 'magnetic':
296        cm = cx+4
297    oAmat,oBmat = cell2AB(oldPhase['General']['Cell'][1:7])
298    nAmat,nBmat = cell2AB(newPhase['General']['Cell'][1:7])
299    SGData = newPhase['General']['SGData']
300    invTrans = nl.inv(Trans)
301    newAtoms,atCodes = FillUnitCell(oldPhase)
302    Unit =[abs(int(max(unit))-1) for unit in Trans]
303    for i,unit in enumerate(Unit):
304        if unit > 0:
305            for j in range(unit):
306                moreAtoms = copy.deepcopy(newAtoms)
307                moreCodes = []
308                for atom,code in zip(moreAtoms,atCodes):
309                    atom[cx+i] += 1.
310                    if '+' in code:
311                        cell = list(eval(code.split('+')[1]))
312                        ops = code.split('+')[0]
313                    else:
314                        cell = [0,0,0]
315                        ops = code
316                    cell[i] += 1
317                    moreCodes.append('%s+%d,%d,%d'%(ops,cell[0],cell[1],cell[2])) 
318                newAtoms += moreAtoms
319                atCodes += moreCodes
320    if ifMag:
321        cia += 3
322        cs += 3
323        newPhase['General']['Type'] = 'magnetic'
324        newPhase['General']['AtomPtrs'] = [cx,ct,cs,cia]
325        magAtoms = []
326        magatCodes = []
327        Landeg = 2.0
328        for iat,atom in enumerate(newAtoms):
329            if len(G2elem.GetMFtable([atom[ct],],[Landeg,])):
330                magAtoms.append(atom[:cx+4]+[0.,0.,0.]+atom[cx+4:])
331                magatCodes.append(atCodes[iat])
332        newAtoms = magAtoms
333        atCodes = magatCodes
334        newPhase['Draw Atoms'] = []
335    for atom in newAtoms:
336        atom[cx:cx+3] = TransformXYZ(atom[cx:cx+3]-Uvec,invTrans.T,Vvec)%1.
337        if atom[cia] == 'A':
338            atom[cia+2:cia+8] = TransformU6(atom[cia+2:cia+8],Trans)
339        atom[cs:cs+2] = G2spc.SytSym(atom[cx:cx+3],SGData)[:2]
340        atom[cia+8] = ran.randint(0,sys.maxsize)
341        if cm:
342            mag = np.sqrt(np.sum(np.array(atom[cm:cm+3])**2))
343            if mag:
344                mom = np.inner(np.array(atom[cm:cm+3]),oBmat)
345                mom = np.inner(mom,invTrans)
346                mom = np.inner(mom,nAmat)
347                mom /= np.sqrt(np.sum(mom**2))
348                atom[cm:cm+3] = mom*mag
349    newPhase['Atoms'] = newAtoms
350    newPhase['Atoms'],atCodes = GetUnique(newPhase,atCodes)
351    newPhase['Drawing'] = []
352    newPhase['ranId'] = ran.randint(0,sys.maxsize)
353    return newPhase,atCodes
354   
355def FillUnitCell(Phase):
356    Atoms = copy.deepcopy(Phase['Atoms'])
357    atomData = []
358    atCodes = []
359    SGData = Phase['General']['SGData']
360    SpnFlp = SGData.get('SpnFlp',[])
361    Amat,Bmat = cell2AB(Phase['General']['Cell'][1:7])
362    cx,ct,cs,cia = Phase['General']['AtomPtrs']
363    cm = 0
364    if Phase['General']['Type'] == 'magnetic':
365        cm = cx+4
366    for iat,atom in enumerate(Atoms):
367        XYZ = np.array(atom[cx:cx+3])
368        xyz = XYZ%1.
369        if atom[cia] == 'A':
370            Uij = atom[cia+2:cia+8]
371            result = G2spc.GenAtom(xyz,SGData,False,Uij,True)
372            for item in result:
373                if item[0][2] >= .95: item[0][2] -= 1.
374                atom[cx:cx+3] = item[0]
375                atom[cia+2:cia+8] = item[1]
376                if cm:
377                    Opr = abs(item[2])%100
378                    M = SGData['SGOps'][Opr-1][0]
379                    opNum = G2spc.GetOpNum(item[2],SGData)
380                    mom = np.inner(np.array(atom[cm:cm+3]),Bmat)
381                    atom[cm:cm+3] = np.inner(np.inner(mom,M),Amat)*nl.det(M)*SpnFlp[opNum-1]
382                atCodes.append('%d:%s'%(iat,str(item[2])))
383                atomData.append(atom[:cia+9])  #not SS stuff
384        else:
385            result = G2spc.GenAtom(xyz,SGData,False,Move=True)
386            for item in result:
387                if item[0][2] >= .95: item[0][2] -= 1.
388                atom[cx:cx+3] = item[0]
389                if cm:
390                    Opr = abs(item[1])%100
391                    M = SGData['SGOps'][Opr-1][0]
392                    opNum = G2spc.GetOpNum(item[1],SGData)
393                    mom = np.inner(np.array(atom[cm:cm+3]),Bmat)
394                    atom[cm:cm+3] = np.inner(np.inner(mom,M),Amat)*nl.det(M)*SpnFlp[opNum-1]
395                atCodes.append('%d:%s'%(iat,str(item[1])))
396                atomData.append(atom[:cia+9])  #not SS stuff
397           
398    return atomData,atCodes
399       
400def GetUnique(Phase,atCodes):
401   
402    def noDuplicate(xyzA,XYZ):
403        if True in [np.allclose(xyzA,xyzB,atol=0.0002) for xyzB in XYZ]:
404            return False
405        return True
406
407    cx,ct,cs,cia = Phase['General']['AtomPtrs']
408    SGData = Phase['General']['SGData']
409    Atoms = Phase['Atoms']
410    Ind = len(Atoms)
411    newAtoms = []
412    newAtCodes = []
413    Indx = {}
414    XYZ = {}
415    for ind in range(Ind):
416        XYZ[ind] = np.array(Atoms[ind][cx:cx+3])%1.
417        Indx[ind] = True
418    for ind in range(Ind):
419        if Indx[ind]:
420            xyz = XYZ[ind]
421            for jnd in range(Ind):
422                if Atoms[ind][ct-1] == Atoms[jnd][ct-1]:
423                    if ind != jnd and Indx[jnd]:                       
424                        Equiv = G2spc.GenAtom(XYZ[jnd],SGData,Move=True)
425                        xyzs = np.array([equiv[0] for equiv in Equiv])
426                        Indx[jnd] = noDuplicate(xyz,xyzs)
427    Ind = []
428    for ind in Indx:
429        if Indx[ind]:
430            newAtoms.append(Atoms[ind])
431            newAtCodes.append(atCodes[ind])
432    return newAtoms,newAtCodes
433           
434def calc_rVsq(A):
435    """Compute the square of the reciprocal lattice volume (1/V**2) from A'
436
437    """
438    G,g = A2Gmat(A)
439    rVsq = nl.det(G)
440    if rVsq < 0:
441        return 1
442    return rVsq
443   
444def calc_rV(A):
445    """Compute the reciprocal lattice volume (V*) from A
446    """
447    return np.sqrt(calc_rVsq(A))
448   
449def calc_V(A):
450    """Compute the real lattice volume (V) from A
451    """
452    return 1./calc_rV(A)
453
454def A2invcell(A):
455    """Compute reciprocal unit cell constants from A
456    returns tuple with a*,b*,c*,alpha*, beta*, gamma* (degrees)
457    """
458    G,g = A2Gmat(A)
459    return Gmat2cell(G)
460   
461def Gmat2AB(G):
462    """Computes orthogonalization matrix from reciprocal metric tensor G
463
464    :returns: tuple of two 3x3 numpy arrays (A,B)
465
466       * A for crystal to Cartesian transformations (A*x = np.inner(A,x) = X)
467       * B (= inverse of A) for Cartesian to crystal transformation (B*X = np.inner(B,X) = x)
468
469    """
470    cellstar = Gmat2cell(G)
471    g = nl.inv(G)
472    cell = Gmat2cell(g)
473    A = np.zeros(shape=(3,3))
474    # from Giacovazzo (Fundamentals 2nd Ed.) p.75
475    A[0][0] = cell[0]                # a
476    A[0][1] = cell[1]*cosd(cell[5])  # b cos(gamma)
477    A[0][2] = cell[2]*cosd(cell[4])  # c cos(beta)
478    A[1][1] = cell[1]*sind(cell[5])  # b sin(gamma)
479    A[1][2] = -cell[2]*cosd(cellstar[3])*sind(cell[4]) # - c cos(alpha*) sin(beta)
480    A[2][2] = 1./cellstar[2]         # 1/c*
481    B = nl.inv(A)
482    return A,B
483   
484def cell2AB(cell):
485    """Computes orthogonalization matrix from unit cell constants
486
487    :param tuple cell: a,b,c, alpha, beta, gamma (degrees)
488    :returns: tuple of two 3x3 numpy arrays (A,B)
489       A for crystal to Cartesian transformations A*x = np.inner(A,x) = X
490       B (= inverse of A) for Cartesian to crystal transformation B*X = np.inner(B,X) = x
491    """
492    G,g = cell2Gmat(cell) 
493    cellstar = Gmat2cell(G)
494    A = np.zeros(shape=(3,3))
495    # from Giacovazzo (Fundamentals 2nd Ed.) p.75
496    A[0][0] = cell[0]                # a
497    A[0][1] = cell[1]*cosd(cell[5])  # b cos(gamma)
498    A[0][2] = cell[2]*cosd(cell[4])  # c cos(beta)
499    A[1][1] = cell[1]*sind(cell[5])  # b sin(gamma)
500    A[1][2] = -cell[2]*cosd(cellstar[3])*sind(cell[4]) # - c cos(alpha*) sin(beta)
501    A[2][2] = 1./cellstar[2]         # 1/c*
502    B = nl.inv(A)
503    return A,B
504   
505def HKL2SpAng(H,cell,SGData):
506    """Computes spherical coords for hkls; view along 001
507
508    :param array H: arrays of hkl
509    :param tuple cell: a,b,c, alpha, beta, gamma (degrees)
510    :param dict SGData: space group dictionary
511    :returns: arrays of r,phi,psi (radius,inclination,azimuth) about 001
512    """
513    A,B = cell2AB(cell)
514    xH = np.inner(B.T,H)
515    r = np.sqrt(np.sum(xH**2,axis=0))
516    phi = acosd(xH[2]/r)
517    psi = atan2d(xH[1],xH[0])
518    phi = np.where(phi>90.,180.-phi,phi)
519#    GSASIIpath.IPyBreak()
520    return r,phi,psi
521   
522def U6toUij(U6):
523    """Fill matrix (Uij) from U6 = [U11,U22,U33,U12,U13,U23]
524    NB: there is a non numpy version in GSASIIspc: U2Uij
525
526    :param list U6: 6 terms of u11,u22,...
527    :returns:
528        Uij - numpy [3][3] array of uij
529    """
530    U = np.array([
531        [U6[0],  U6[3],  U6[4]], 
532        [U6[3],  U6[1],  U6[5]], 
533        [U6[4],  U6[5],  U6[2]]])
534    return U
535
536def UijtoU6(U):
537    """Fill vector [U11,U22,U33,U12,U13,U23] from Uij
538    NB: there is a non numpy version in GSASIIspc: Uij2U
539    """
540    U6 = np.array([U[0][0],U[1][1],U[2][2],U[0][1],U[0][2],U[1][2]])
541    return U6
542
543def betaij2Uij(betaij,G):
544    """
545    Convert beta-ij to Uij tensors
546   
547    :param beta-ij - numpy array [beta-ij]
548    :param G: reciprocal metric tensor
549    :returns: Uij: numpy array [Uij]
550    """
551    ast = np.sqrt(np.diag(G))   #a*, b*, c*
552    Mast = np.multiply.outer(ast,ast)   
553    return R2pisq*UijtoU6(U6toUij(betaij)/Mast)
554   
555def Uij2betaij(Uij,G):
556    """
557    Convert Uij to beta-ij tensors -- stub for eventual completion
558   
559    :param Uij: numpy array [Uij]
560    :param G: reciprocal metric tensor
561    :returns: beta-ij - numpy array [beta-ij]
562    """
563    pass
564   
565def cell2GS(cell):
566    ''' returns Uij to betaij conversion matrix'''
567    G,g = cell2Gmat(cell)
568    GS = G
569    GS[0][1] = GS[1][0] = math.sqrt(GS[0][0]*GS[1][1])
570    GS[0][2] = GS[2][0] = math.sqrt(GS[0][0]*GS[2][2])
571    GS[1][2] = GS[2][1] = math.sqrt(GS[1][1]*GS[2][2])
572    return GS   
573   
574def Uij2Ueqv(Uij,GS,Amat):
575    ''' returns 1/3 trace of diagonalized U matrix'''
576    U = np.multiply(U6toUij(Uij),GS)
577    U = np.inner(Amat,np.inner(U,Amat).T)
578    E,R = nl.eigh(U)
579    return np.sum(E)/3.
580       
581def CosAngle(U,V,G):
582    """ calculate cos of angle between U & V in generalized coordinates
583    defined by metric tensor G
584
585    :param U: 3-vectors assume numpy arrays, can be multiple reflections as (N,3) array
586    :param V: 3-vectors assume numpy arrays, only as (3) vector
587    :param G: metric tensor for U & V defined space assume numpy array
588    :returns:
589        cos(phi)
590    """
591    u = (U.T/np.sqrt(np.sum(np.inner(U,G)*U,axis=1))).T
592    v = V/np.sqrt(np.inner(V,np.inner(G,V)))
593    cosP = np.inner(u,np.inner(G,v))
594    return cosP
595   
596def CosSinAngle(U,V,G):
597    """ calculate sin & cos of angle between U & V in generalized coordinates
598    defined by metric tensor G
599
600    :param U: 3-vectors assume numpy arrays
601    :param V: 3-vectors assume numpy arrays
602    :param G: metric tensor for U & V defined space assume numpy array
603    :returns:
604        cos(phi) & sin(phi)
605    """
606    u = U/np.sqrt(np.inner(U,np.inner(G,U)))
607    v = V/np.sqrt(np.inner(V,np.inner(G,V)))
608    cosP = np.inner(u,np.inner(G,v))
609    sinP = np.sqrt(max(0.0,1.0-cosP**2))
610    return cosP,sinP
611   
612def criticalEllipse(prob):
613    """
614    Calculate critical values for probability ellipsoids from probability
615    """
616    if not ( 0.01 <= prob < 1.0):
617        return 1.54 
618    coeff = np.array([6.44988E-09,4.16479E-07,1.11172E-05,1.58767E-04,0.00130554,
619        0.00604091,0.0114921,-0.040301,-0.6337203,1.311582])
620    llpr = math.log(-math.log(prob))
621    return np.polyval(coeff,llpr)
622   
623def CellBlock(nCells):
624    """
625    Generate block of unit cells n*n*n on a side; [0,0,0] centered, n = 2*nCells+1
626    currently only works for nCells = 0 or 1 (not >1)
627    """
628    if nCells:
629        N = 2*nCells+1
630        N2 = N*N
631        N3 = N*N*N
632        cellArray = []
633        A = np.array(range(N3))
634        cellGen = np.array([A//N2-1,A//N%N-1,A%N-1]).T
635        for cell in cellGen:
636            cellArray.append(cell)
637        return cellArray
638    else:
639        return [0,0,0]
640       
641def CellAbsorption(ElList,Volume):
642    '''Compute unit cell absorption
643
644    :param dict ElList: dictionary of element contents including mu and
645      number of atoms be cell
646    :param float Volume: unit cell volume
647    :returns: mu-total/Volume
648    '''
649    muT = 0
650    for El in ElList:
651        muT += ElList[El]['mu']*ElList[El]['FormulaNo']
652    return muT/Volume
653   
654#Permutations and Combinations
655# Four routines: combinations,uniqueCombinations, selections & permutations
656#These taken from Python Cookbook, 2nd Edition. 19.15 p724-726
657#   
658def _combinators(_handle, items, n):
659    """ factored-out common structure of all following combinators """
660    if n==0:
661        yield [ ]
662        return
663    for i, item in enumerate(items):
664        this_one = [ item ]
665        for cc in _combinators(_handle, _handle(items, i), n-1):
666            yield this_one + cc
667def combinations(items, n):
668    """ take n distinct items, order matters """
669    def skipIthItem(items, i):
670        return items[:i] + items[i+1:]
671    return _combinators(skipIthItem, items, n)
672def uniqueCombinations(items, n):
673    """ take n distinct items, order is irrelevant """
674    def afterIthItem(items, i):
675        return items[i+1:]
676    return _combinators(afterIthItem, items, n)
677def selections(items, n):
678    """ take n (not necessarily distinct) items, order matters """
679    def keepAllItems(items, i):
680        return items
681    return _combinators(keepAllItems, items, n)
682def permutations(items):
683    """ take all items, order matters """
684    return combinations(items, len(items))
685
686#reflection generation routines
687#for these: H = [h,k,l]; A is as used in calc_rDsq; G - inv metric tensor, g - metric tensor;
688#           cell - a,b,c,alp,bet,gam in A & deg
689   
690def Pos2dsp(Inst,pos):
691    ''' convert powder pattern position (2-theta or TOF, musec) to d-spacing
692    '''
693    if 'C' in Inst['Type'][0] or 'PKS' in Inst['Type'][0]:
694        wave = G2mth.getWave(Inst)
695        return wave/(2.0*sind((pos-Inst.get('Zero',[0,0])[1])/2.0))
696    else:   #'T'OF - ignore difB
697        return TOF2dsp(Inst,pos)
698       
699def TOF2dsp(Inst,Pos):
700    ''' convert powder pattern TOF, musec to d-spacing by successive approximation
701    Pos can be numpy array
702    '''
703    def func(d,pos,Inst):       
704        return (pos-Inst['difA'][1]*d**2-Inst['Zero'][1]-Inst['difB'][1]/d)/Inst['difC'][1]
705    dsp0 = np.ones_like(Pos)
706    N = 0
707    while True:      #successive approximations
708        dsp = func(dsp0,Pos,Inst)
709        if np.allclose(dsp,dsp0,atol=0.000001):
710            return dsp
711        dsp0 = dsp
712        N += 1
713        if N > 10:
714            return dsp
715   
716def Dsp2pos(Inst,dsp):
717    ''' convert d-spacing to powder pattern position (2-theta or TOF, musec)
718    '''
719    if 'C' in Inst['Type'][0] or 'PKS' in Inst['Type'][0]:
720        wave = G2mth.getWave(Inst)
721        val = min(0.995,wave/(2.*dsp))  #set max at 168deg
722        pos = 2.0*asind(val)+Inst.get('Zero',[0,0])[1]             
723    else:   #'T'OF
724        pos = Inst['difC'][1]*dsp+Inst['Zero'][1]+Inst['difA'][1]*dsp**2+Inst.get('difB',[0,0,False])[1]/dsp
725    return pos
726   
727def getPeakPos(dataType,parmdict,dsp):
728    ''' convert d-spacing to powder pattern position (2-theta or TOF, musec)
729    '''
730    if 'C' in dataType:
731        pos = 2.0*asind(parmdict['Lam']/(2.*dsp))+parmdict['Zero']
732    else:   #'T'OF
733        pos = parmdict['difC']*dsp+parmdict['difA']*dsp**2+parmdict['difB']/dsp+parmdict['Zero']
734    return pos
735                   
736def calc_rDsq(H,A):
737    'needs doc string'
738    rdsq = H[0]*H[0]*A[0]+H[1]*H[1]*A[1]+H[2]*H[2]*A[2]+H[0]*H[1]*A[3]+H[0]*H[2]*A[4]+H[1]*H[2]*A[5]
739    return rdsq
740   
741def calc_rDsq2(H,G):
742    'needs doc string'
743    return np.inner(H,np.inner(G,H))
744   
745def calc_rDsqSS(H,A,vec):
746    'needs doc string'
747    rdsq = calc_rDsq(H[:3]+(H[3]*vec).T,A)
748    return rdsq
749       
750def calc_rDsqZ(H,A,Z,tth,lam):
751    'needs doc string'
752    rdsq = calc_rDsq(H,A)+Z*sind(tth)*2.0*rpd/lam**2
753    return rdsq
754       
755def calc_rDsqZSS(H,A,vec,Z,tth,lam):
756    'needs doc string'
757    rdsq = calc_rDsq(H[:3]+(H[3][:,np.newaxis]*vec).T,A)+Z*sind(tth)*2.0*rpd/lam**2
758    return rdsq
759       
760def calc_rDsqT(H,A,Z,tof,difC):
761    'needs doc string'
762    rdsq = calc_rDsq(H,A)+Z/difC
763    return rdsq
764       
765def calc_rDsqTSS(H,A,vec,Z,tof,difC):
766    'needs doc string'
767    rdsq = calc_rDsq(H[:3]+(H[3][:,np.newaxis]*vec).T,A)+Z/difC
768    return rdsq
769   
770def PlaneIntercepts(Amat,Bmat,H,phase,stack):
771    ''' find unit cell intercepts for a stack of hkl planes
772    '''
773    Steps = range(-1,2,2)
774    if stack:
775        Steps = range(-10,10,1)
776    Stack = []
777    Ux = np.array([[0,0],[1,0],[1,1],[0,1]])
778    for step in Steps:
779        HX = []
780        for i in [0,1,2]:
781            if H[i]:
782               h,k,l = [(i+1)%3,(i+2)%3,(i+3)%3]
783               for j in [0,1,2,3]:
784                    hx = [0,0,0]
785                    intcpt = (phase/360.+step-H[h]*Ux[j,0]-H[k]*Ux[j,1])/H[l]
786                    if 0. <= intcpt <= 1.:                       
787                        hx[h] = Ux[j,0]
788                        hx[k] = Ux[j,1]
789                        hx[l] = intcpt
790                        HX.append(hx)
791        if len(HX)> 2:
792            HX = np.array(HX)
793            DX = np.inner(HX-HX[0],Amat)
794            D = np.sqrt(np.sum(DX**2,axis=1))
795            Dsort = np.argsort(D)
796            HX = HX[Dsort]
797            DX = DX[Dsort]
798            D = D[Dsort]
799            DX[1:,:] = DX[1:,:]/D[1:,nxs]
800            A = 2.*np.ones(HX.shape[0])
801            A[1:] = [np.dot(DX[1],dx) for dx in DX[1:]]
802            HX = HX[np.argsort(A)]
803#            GSASIIpath.IPyBreak()
804            Stack.append(HX)
805    return Stack
806       
807def MaxIndex(dmin,A):
808    'needs doc string'
809    Hmax = [0,0,0]
810    try:
811        cell = A2cell(A)
812    except:
813        cell = [1.,1.,1.,90.,90.,90.]
814    for i in range(3):
815        Hmax[i] = int(round(cell[i]/dmin))
816    return Hmax
817   
818def transposeHKLF(transMat,Super,refList):
819    ''' Apply transformation matrix to hkl(m)
820    param: transmat: 3x3 or 4x4 array
821    param: Super: 0 or 1 for extra index
822    param: refList list of h,k,l,....
823    return: newRefs transformed list of h',k',l',,,
824    return: badRefs list of noninteger h',k',l'...
825    '''
826    newRefs = np.copy(refList)
827    badRefs = []
828    for H in newRefs:
829        newH = np.inner(transMat,H[:3+Super])
830        H[:3+Super] = np.rint(newH)
831        if not np.allclose(newH,H[:3+Super],atol=0.01):
832            badRefs.append(newH)
833    return newRefs,badRefs
834   
835def sortHKLd(HKLd,ifreverse,ifdup,ifSS=False):
836    '''sort reflection list on d-spacing; can sort in either order
837
838    :param HKLd: a list of [h,k,l,d,...];
839    :param ifreverse: True for largest d first
840    :param ifdup: True if duplicate d-spacings allowed
841    :return: sorted reflection list
842    '''
843    T = []
844    N = 3
845    if ifSS:
846        N = 4
847    for i,H in enumerate(HKLd):
848        if ifdup:
849            T.append((H[N],i))
850        else:
851            T.append(H[N])           
852    D = dict(zip(T,HKLd))
853    T.sort()
854    if ifreverse:
855        T.reverse()
856    X = []
857    okey = ''
858    for key in T: 
859        if key != okey: X.append(D[key])    #remove duplicate d-spacings
860        okey = key
861    return X
862   
863def SwapIndx(Axis,H):
864    'needs doc string'
865    if Axis in [1,-1]:
866        return H
867    elif Axis in [2,-3]:
868        return [H[1],H[2],H[0]]
869    else:
870        return [H[2],H[0],H[1]]
871       
872def Rh2Hx(Rh):
873    'needs doc string'
874    Hx = [0,0,0]
875    Hx[0] = Rh[0]-Rh[1]
876    Hx[1] = Rh[1]-Rh[2]
877    Hx[2] = np.sum(Rh)
878    return Hx
879   
880def Hx2Rh(Hx):
881    'needs doc string'
882    Rh = [0,0,0]
883    itk = -Hx[0]+Hx[1]+Hx[2]
884    if itk%3 != 0:
885        return 0        #error - not rhombohedral reflection
886    else:
887        Rh[1] = itk//3
888        Rh[0] = Rh[1]+Hx[0]
889        Rh[2] = Rh[1]-Hx[1]
890        if Rh[0] < 0:
891            for i in range(3):
892                Rh[i] = -Rh[i]
893        return Rh
894       
895def CentCheck(Cent,H):
896    'needs doc string'
897    h,k,l = H
898    if Cent == 'A' and (k+l)%2:
899        return False
900    elif Cent == 'B' and (h+l)%2:
901        return False
902    elif Cent == 'C' and (h+k)%2:
903        return False
904    elif Cent == 'I' and (h+k+l)%2:
905        return False
906    elif Cent == 'F' and ((h+k)%2 or (h+l)%2 or (k+l)%2):
907        return False
908    elif Cent == 'R' and (-h+k+l)%3:
909        return False
910    else:
911        return True
912                                   
913def GetBraviasNum(center,system):
914    """Determine the Bravais lattice number, as used in GenHBravais
915   
916    :param center: one of: 'P', 'C', 'I', 'F', 'R' (see SGLatt from GSASIIspc.SpcGroup)
917    :param system: one of 'cubic', 'hexagonal', 'tetragonal', 'orthorhombic', 'trigonal' (for R)
918      'monoclinic', 'triclinic' (see SGSys from GSASIIspc.SpcGroup)
919    :return: a number between 0 and 13
920      or throws a ValueError exception if the combination of center, system is not found (i.e. non-standard)
921
922    """
923    if center.upper() == 'F' and system.lower() == 'cubic':
924        return 0
925    elif center.upper() == 'I' and system.lower() == 'cubic':
926        return 1
927    elif center.upper() == 'P' and system.lower() == 'cubic':
928        return 2
929    elif center.upper() == 'R' and system.lower() == 'trigonal':
930        return 3
931    elif center.upper() == 'P' and system.lower() == 'hexagonal':
932        return 4
933    elif center.upper() == 'I' and system.lower() == 'tetragonal':
934        return 5
935    elif center.upper() == 'P' and system.lower() == 'tetragonal':
936        return 6
937    elif center.upper() == 'F' and system.lower() == 'orthorhombic':
938        return 7
939    elif center.upper() == 'I' and system.lower() == 'orthorhombic':
940        return 8
941    elif center.upper() == 'C' and system.lower() == 'orthorhombic':
942        return 9
943    elif center.upper() == 'P' and system.lower() == 'orthorhombic':
944        return 10
945    elif center.upper() == 'C' and system.lower() == 'monoclinic':
946        return 11
947    elif center.upper() == 'P' and system.lower() == 'monoclinic':
948        return 12
949    elif center.upper() == 'P' and system.lower() == 'triclinic':
950        return 13
951    raise ValueError('non-standard Bravais lattice center=%s, cell=%s' % (center,system))
952
953def GenHBravais(dmin,Bravais,A):
954    """Generate the positionally unique powder diffraction reflections
955     
956    :param dmin: minimum d-spacing in A
957    :param Bravais: lattice type (see GetBraviasNum). Bravais is one of:
958   
959            * 0 F cubic
960            * 1 I cubic
961            * 2 P cubic
962            * 3 R hexagonal (trigonal not rhombohedral)
963            * 4 P hexagonal
964            * 5 I tetragonal
965            * 6 P tetragonal
966            * 7 F orthorhombic
967            * 8 I orthorhombic
968            * 9 C orthorhombic
969            * 10 P orthorhombic
970            * 11 C monoclinic
971            * 12 P monoclinic
972            * 13 P triclinic
973           
974    :param A: reciprocal metric tensor elements as [G11,G22,G33,2*G12,2*G13,2*G23]
975    :return: HKL unique d list of [h,k,l,d,-1] sorted with largest d first
976           
977    """
978    if Bravais in [9,11]:
979        Cent = 'C'
980    elif Bravais in [1,5,8]:
981        Cent = 'I'
982    elif Bravais in [0,7]:
983        Cent = 'F'
984    elif Bravais in [3]:
985        Cent = 'R'
986    else:
987        Cent = 'P'
988    Hmax = MaxIndex(dmin,A)
989    dminsq = 1./(dmin**2)
990    HKL = []
991    if Bravais == 13:                       #triclinic
992        for l in range(-Hmax[2],Hmax[2]+1):
993            for k in range(-Hmax[1],Hmax[1]+1):
994                hmin = 0
995                if (k < 0): hmin = 1
996                if (k ==0 and l < 0): hmin = 1
997                for h in range(hmin,Hmax[0]+1):
998                    H=[h,k,l]
999                    rdsq = calc_rDsq(H,A)
1000                    if 0 < rdsq <= dminsq:
1001                        HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
1002    elif Bravais in [11,12]:                #monoclinic - b unique
1003        Hmax = SwapIndx(2,Hmax)
1004        for h in range(Hmax[0]+1):
1005            for k in range(-Hmax[1],Hmax[1]+1):
1006                lmin = 0
1007                if k < 0:lmin = 1
1008                for l in range(lmin,Hmax[2]+1):
1009                    [h,k,l] = SwapIndx(-2,[h,k,l])
1010                    H = []
1011                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
1012                    if H:
1013                        rdsq = calc_rDsq(H,A)
1014                        if 0 < rdsq <= dminsq:
1015                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
1016                    [h,k,l] = SwapIndx(2,[h,k,l])
1017    elif Bravais in [7,8,9,10]:            #orthorhombic
1018        for h in range(Hmax[0]+1):
1019            for k in range(Hmax[1]+1):
1020                for l in range(Hmax[2]+1):
1021                    H = []
1022                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
1023                    if H:
1024                        rdsq = calc_rDsq(H,A)
1025                        if 0 < rdsq <= dminsq:
1026                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
1027    elif Bravais in [5,6]:                  #tetragonal
1028        for l in range(Hmax[2]+1):
1029            for k in range(Hmax[1]+1):
1030                for h in range(k,Hmax[0]+1):
1031                    H = []
1032                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
1033                    if H:
1034                        rdsq = calc_rDsq(H,A)
1035                        if 0 < rdsq <= dminsq:
1036                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
1037    elif Bravais in [3,4]:
1038        lmin = 0
1039        if Bravais == 3: lmin = -Hmax[2]                  #hexagonal/trigonal
1040        for l in range(lmin,Hmax[2]+1):
1041            for k in range(Hmax[1]+1):
1042                hmin = k
1043                if l < 0: hmin += 1
1044                for h in range(hmin,Hmax[0]+1):
1045                    H = []
1046                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
1047                    if H:
1048                        rdsq = calc_rDsq(H,A)
1049                        if 0 < rdsq <= dminsq:
1050                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
1051
1052    else:                                   #cubic
1053        for l in range(Hmax[2]+1):
1054            for k in range(l,Hmax[1]+1):
1055                for h in range(k,Hmax[0]+1):
1056                    H = []
1057                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
1058                    if H:
1059                        rdsq = calc_rDsq(H,A)
1060                        if 0 < rdsq <= dminsq:
1061                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
1062    return sortHKLd(HKL,True,False)
1063   
1064def getHKLmax(dmin,SGData,A):
1065    'finds maximum allowed hkl for given A within dmin'
1066    SGLaue = SGData['SGLaue']
1067    if SGLaue in ['3R','3mR']:        #Rhombohedral axes
1068        Hmax = [0,0,0]
1069        cell = A2cell(A)
1070        aHx = cell[0]*math.sqrt(2.0*(1.0-cosd(cell[3])))
1071        cHx = cell[0]*math.sqrt(3.0*(1.0+2.0*cosd(cell[3])))
1072        Hmax[0] = Hmax[1] = int(round(aHx/dmin))
1073        Hmax[2] = int(round(cHx/dmin))
1074        #print Hmax,aHx,cHx
1075    else:                           # all others
1076        Hmax = MaxIndex(dmin,A)
1077    return Hmax
1078   
1079def GenHLaue(dmin,SGData,A):
1080    """Generate the crystallographically unique powder diffraction reflections
1081    for a lattice and Bravais type
1082   
1083    :param dmin: minimum d-spacing
1084    :param SGData: space group dictionary with at least
1085   
1086        * 'SGLaue': Laue group symbol: one of '-1','2/m','mmm','4/m','6/m','4/mmm','6/mmm', '3m1', '31m', '3', '3R', '3mR', 'm3', 'm3m'
1087        * 'SGLatt': lattice centering: one of 'P','A','B','C','I','F'
1088        * 'SGUniq': code for unique monoclinic axis one of 'a','b','c' (only if 'SGLaue' is '2/m') otherwise an empty string
1089       
1090    :param A: reciprocal metric tensor elements as [G11,G22,G33,2*G12,2*G13,2*G23]
1091    :return: HKL = list of [h,k,l,d] sorted with largest d first and is unique
1092            part of reciprocal space ignoring anomalous dispersion
1093           
1094    """
1095    import math
1096    SGLaue = SGData['SGLaue']
1097    SGLatt = SGData['SGLatt']
1098    SGUniq = SGData['SGUniq']
1099    #finds maximum allowed hkl for given A within dmin
1100    Hmax = getHKLmax(dmin,SGData,A)
1101       
1102    dminsq = 1./(dmin**2)
1103    HKL = []
1104    if SGLaue == '-1':                       #triclinic
1105        for l in range(-Hmax[2],Hmax[2]+1):
1106            for k in range(-Hmax[1],Hmax[1]+1):
1107                hmin = 0
1108                if (k < 0) or (k ==0 and l < 0): hmin = 1
1109                for h in range(hmin,Hmax[0]+1):
1110                    H = []
1111                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
1112                    if H:
1113                        rdsq = calc_rDsq(H,A)
1114                        if 0 < rdsq <= dminsq:
1115                            HKL.append([h,k,l,1./math.sqrt(rdsq)])
1116    elif SGLaue == '2/m':                #monoclinic
1117        axisnum = 1 + ['a','b','c'].index(SGUniq)
1118        Hmax = SwapIndx(axisnum,Hmax)
1119        for h in range(Hmax[0]+1):
1120            for k in range(-Hmax[1],Hmax[1]+1):
1121                lmin = 0
1122                if k < 0:lmin = 1
1123                for l in range(lmin,Hmax[2]+1):
1124                    [h,k,l] = SwapIndx(-axisnum,[h,k,l])
1125                    H = []
1126                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
1127                    if H:
1128                        rdsq = calc_rDsq(H,A)
1129                        if 0 < rdsq <= dminsq:
1130                            HKL.append([h,k,l,1./math.sqrt(rdsq)])
1131                    [h,k,l] = SwapIndx(axisnum,[h,k,l])
1132    elif SGLaue in ['mmm','4/m','6/m']:            #orthorhombic
1133        for l in range(Hmax[2]+1):
1134            for h in range(Hmax[0]+1):
1135                kmin = 1
1136                if SGLaue == 'mmm' or h ==0: kmin = 0
1137                for k in range(kmin,Hmax[1]+1):
1138                    H = []
1139                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
1140                    if H:
1141                        rdsq = calc_rDsq(H,A)
1142                        if 0 < rdsq <= dminsq:
1143                            HKL.append([h,k,l,1./math.sqrt(rdsq)])
1144    elif SGLaue in ['4/mmm','6/mmm']:                  #tetragonal & hexagonal
1145        for l in range(Hmax[2]+1):
1146            for h in range(Hmax[0]+1):
1147                for k in range(h+1):
1148                    H = []
1149                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
1150                    if H:
1151                        rdsq = calc_rDsq(H,A)
1152                        if 0 < rdsq <= dminsq:
1153                            HKL.append([h,k,l,1./math.sqrt(rdsq)])
1154    elif SGLaue in ['3m1','31m','3','3R','3mR']:                  #trigonals
1155        for l in range(-Hmax[2],Hmax[2]+1):
1156            hmin = 0
1157            if l < 0: hmin = 1
1158            for h in range(hmin,Hmax[0]+1):
1159                if SGLaue in ['3R','3']:
1160                    kmax = h
1161                    kmin = -int((h-1.)/2.)
1162                else:
1163                    kmin = 0
1164                    kmax = h
1165                    if SGLaue in ['3m1','3mR'] and l < 0: kmax = h-1
1166                    if SGLaue == '31m' and l < 0: kmin = 1
1167                for k in range(kmin,kmax+1):
1168                    H = []
1169                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
1170                    if SGLaue in ['3R','3mR']:
1171                        H = Hx2Rh(H)
1172                    if H:
1173                        rdsq = calc_rDsq(H,A)
1174                        if 0 < rdsq <= dminsq:
1175                            HKL.append([H[0],H[1],H[2],1./math.sqrt(rdsq)])
1176    else:                                   #cubic
1177        for h in range(Hmax[0]+1):
1178            for k in range(h+1):
1179                lmin = 0
1180                lmax = k
1181                if SGLaue =='m3':
1182                    lmax = h-1
1183                    if h == k: lmax += 1
1184                for l in range(lmin,lmax+1):
1185                    H = []
1186                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
1187                    if H:
1188                        rdsq = calc_rDsq(H,A)
1189                        if 0 < rdsq <= dminsq:
1190                            HKL.append([h,k,l,1./math.sqrt(rdsq)])
1191    return sortHKLd(HKL,True,True)
1192   
1193def GenPfHKLs(nMax,SGData,A):   
1194    """Generate the unique pole figure reflections for a lattice and Bravais type.
1195    Min d-spacing=1.0A & no more than nMax returned
1196   
1197    :param nMax: maximum number of hkls returned
1198    :param SGData: space group dictionary with at least
1199   
1200        * 'SGLaue': Laue group symbol: one of '-1','2/m','mmm','4/m','6/m','4/mmm','6/mmm', '3m1', '31m', '3', '3R', '3mR', 'm3', 'm3m'
1201        * 'SGLatt': lattice centering: one of 'P','A','B','C','I','F'
1202        * 'SGUniq': code for unique monoclinic axis one of 'a','b','c' (only if 'SGLaue' is '2/m') otherwise an empty string
1203       
1204    :param A: reciprocal metric tensor elements as [G11,G22,G33,2*G12,2*G13,2*G23]
1205    :return: HKL = list of 'h k l' strings sorted with largest d first; no duplicate zones
1206           
1207    """
1208    HKL = np.array(GenHLaue(1.0,SGData,A)).T[:3].T     #strip d-spacings
1209    N = min(nMax,len(HKL))
1210    return ['%d %d %d'%(h[0],h[1],h[2]) for h in HKL[:N]]       
1211
1212def GenSSHLaue(dmin,SGData,SSGData,Vec,maxH,A):
1213    'needs a doc string'
1214    HKLs = []
1215    vec = np.array(Vec)
1216    vstar = np.sqrt(calc_rDsq(vec,A))     #find extra needed for -n SS reflections
1217    dvec = 1./(maxH*vstar+1./dmin)
1218    HKL = GenHLaue(dvec,SGData,A)       
1219    SSdH = [vec*h for h in range(-maxH,maxH+1)]
1220    SSdH = dict(zip(range(-maxH,maxH+1),SSdH))
1221    for h,k,l,d in HKL:
1222        ext = G2spc.GenHKLf([h,k,l],SGData)[0]  #h,k,l must be integral values here
1223        if not ext and d >= dmin:
1224            HKLs.append([h,k,l,0,d])
1225        for dH in SSdH:
1226            if dH:
1227                DH = SSdH[dH]
1228                H = [h+DH[0],k+DH[1],l+DH[2]]
1229                d = 1./np.sqrt(calc_rDsq(H,A))
1230                if d >= dmin:
1231                    HKLM = np.array([h,k,l,dH])
1232                    if G2spc.checkSSLaue([h,k,l,dH],SGData,SSGData) and G2spc.checkSSextc(HKLM,SSGData):
1233                        HKLs.append([h,k,l,dH,d])   
1234    return HKLs
1235   
1236def LaueUnique2(SGData,refList):
1237    ''' Impose Laue symmetry on hkl
1238   
1239    :param SGData: space group data from 'P '+Laue
1240    :param HKLF: np.array([[h,k,l,...]]) reflection set to be converted
1241   
1242    :return: HKLF new reflection array with imposed Laue symmetry
1243    '''
1244    for ref in refList:
1245        H = ref[:3]
1246        Uniq = G2spc.GenHKLf(H,SGData)[2]
1247        Uniq = G2mth.sortArray(G2mth.sortArray(G2mth.sortArray(Uniq,2),1),0)
1248        ref[:3] = Uniq[-1]
1249    return refList
1250   
1251def LaueUnique(Laue,HKLF):
1252    ''' Impose Laue symmetry on hkl
1253   
1254    :param str Laue: Laue symbol, as below
1255   
1256      centrosymmetric Laue groups::
1257       
1258            ['-1','2/m','112/m','2/m11','mmm','-42m','-4m2','4/mmm','-3',
1259            '-31m','-3m1','6/m','6/mmm','m3','m3m']
1260     
1261      noncentrosymmetric Laue groups::
1262     
1263           ['1','2','211','112','m','m11','11m','222','mm2','m2m','2mm',
1264           '4','-4','422','4mm','3','312','321','31m','3m1','6','-6',
1265           '622','6mm','-62m','-6m2','23','432','-43m']
1266     
1267    :param HKLF: np.array([[h,k,l,...]]) reflection set to be converted
1268   
1269    :returns: HKLF new reflection array with imposed Laue symmetry
1270    '''
1271   
1272    HKLFT = HKLF.T
1273    mat41 = np.array([[0,1,0],[-1,0,0],[0,0,1]])    #hkl -> k,-h,l
1274    mat43 = np.array([[0,-1,0],[1,0,0],[0,0,1]])    #hkl -> -k,h,l
1275    mat4bar = np.array([[0,-1,0],[1,0,0],[0,0,-1]]) #hkl -> k,-h,-l
1276    mat31 = np.array([[-1,-1,0],[1,0,0],[0,0,1]])   #hkl -> ihl = -h-k,h,l
1277    mat32 = np.array([[0,1,0],[-1,-1,0],[0,0,1]])   #hkl -> kil = k,-h-k,l
1278    matd3 = np.array([[0,1,0],[0,0,1],[1,0,0]])     #hkl -> k,l,h
1279    matd3q = np.array([[0,0,-1],[-1,0,0],[0,1,0]])  #hkl -> -l,-h,k
1280    matd3t = np.array([[0,0,-1],[1,0,0],[0,-1,0]])  #hkl -> -l,h,-k
1281    mat6 = np.array([[1,1,0],[-1,0,0],[0,0,1]])     #hkl -> h+k,-h,l really 65
1282    matdm = np.array([[0,1,0],[1,0,0],[0,0,1]])     #hkl -> k,h,l
1283    matdmp = np.array([[-1,-1,0],[0,1,0],[0,0,1]])  #hkl -> -h-k,k,l
1284    matkm = np.array([[-1,0,0],[1,1,0],[0,0,1]])    #hkl -> -h,h+k,l
1285    matd2 = np.array([[0,1,0],[1,0,0],[0,0,-1]])    #hkl -> k,h,-l
1286    matdm3 = np.array([[1,0,0],[0,0,1],[0,1,0]])    #hkl -> h,l,k
1287    mat2d43 = np.array([[0,1,0],[1,0,0],[0,0,1]])   #hkl -> k,-h,l
1288    matk2 = np.array([[-1,0,0],[1,1,0],[0,0,-1]])   #hkl -> -h,-i,-l
1289    #triclinic
1290    if Laue == '1': #ok
1291        pass
1292    elif Laue == '-1':  #ok
1293        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
1294        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]<0),HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
1295        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
1296    #monoclinic
1297    #noncentrosymmetric - all ok
1298    elif Laue == '2': 
1299        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3])
1300        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3])
1301    elif Laue == '1 1 2':
1302        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1303        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]<0),HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1304    elif Laue == '2 1 1':   
1305        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1306        HKLFT[:3] = np.where((HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1307    elif Laue == 'm':
1308        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1309    elif Laue == 'm 1 1':
1310        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1311    elif Laue == '1 1 m':
1312        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1313    #centrosymmetric - all ok
1314    elif Laue == '2/m 1 1':       
1315        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1316        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1317        HKLFT[:3] = np.where((HKLFT[2]*HKLFT[0]==0)&(HKLFT[1]<0),HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1318    elif Laue == '2/m':
1319        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3])
1320        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1321        HKLFT[:3] = np.where((HKLFT[0]*HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1322    elif Laue == '1 1 2/m':
1323        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1324        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1325        HKLFT[:3] = np.where((HKLFT[1]*HKLFT[2]==0)&(HKLFT[0]<0),HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1326    #orthorhombic
1327    #noncentrosymmetric - all OK
1328    elif Laue == '2 2 2':
1329        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1330        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1331        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1332        HKLFT[:3] = np.where((HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1333    elif Laue == 'm m 2':
1334        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1335        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1336    elif Laue == '2 m m': 
1337        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1338        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1339    elif Laue == 'm 2 m':
1340        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1341        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1342    #centrosymmetric - all ok
1343    elif Laue == 'm m m':
1344        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1345        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1346        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1347    #tetragonal
1348    #noncentrosymmetric - all ok
1349    elif Laue == '4':
1350        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1351        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1352        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]>0),np.squeeze(np.inner(HKLF[:,:3],mat41[nxs,:,:])).T,HKLFT[:3])
1353    elif Laue == '-4': 
1354        HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])     
1355        HKLFT[:3] = np.where(HKLFT[0]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1356        HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])     
1357        HKLFT[:3] = np.where(HKLFT[1]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1358        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1359    elif Laue == '4 2 2':
1360        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1361        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1362        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1363        HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]<HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1364        HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])   #in lieu od 2-fold
1365    elif Laue == '4 m m':
1366        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1367        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1368        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1369        HKLFT[:3] = np.where(HKLFT[0]<HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1370    elif Laue == '-4 2 m':
1371        HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])     
1372        HKLFT[:3] = np.where(HKLFT[0]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1373        HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])     
1374        HKLFT[:3] = np.where(HKLFT[1]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1375        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1376        HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1377        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1378    elif Laue == '-4 m 2':
1379        HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1380        HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])     
1381        HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]<=0),np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1382        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]<0),HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])     
1383        HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]==0),np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1384        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3]) 
1385        HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[0]>HKLFT[1]),np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1386    #centrosymmetric - all ok
1387    elif Laue == '4/m':
1388        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1389        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1390        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1391        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]>0),np.squeeze(np.inner(HKLF[:,:3],mat41[nxs,:,:])).T,HKLFT[:3])
1392    elif Laue == '4/m m m':
1393        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1394        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1395        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])       
1396        HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],mat41[nxs,:,:])).T,HKLFT[:3])
1397        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1398    #trigonal - all hex cell
1399    #noncentrosymmetric - all ok
1400    elif Laue == '3':
1401        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1402        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1403        HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1404    elif Laue == '3 1 2':
1405        HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matk2[nxs,:,:])).T,HKLFT[:3])
1406        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1407        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1408        HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1409        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],matk2[nxs,:,:])).T,HKLFT[:3])
1410    elif Laue == '3 2 1':
1411        HKLFT[:3] = np.where(HKLFT[0]<=-2*HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1412        HKLFT[:3] = np.where(HKLFT[1]<-2*HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1413        HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1414        HKLFT[:3] = np.where((HKLFT[2]>0)&(HKLFT[1]==HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1415        HKLFT[:3] = np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T
1416        HKLFT[:3] = np.where((HKLFT[0]!=0)&(HKLFT[2]>0)&(HKLFT[0]==-2*HKLFT[1]),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1417    elif Laue == '3 1 m':
1418        HKLFT[:3] = np.where(HKLFT[0]>=HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1419        HKLFT[:3] = np.where(2*HKLFT[1]<-HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1420        HKLFT[:3] = np.where(HKLFT[1]>-2*HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matdmp[nxs,:,:])).T,HKLFT[:3])
1421        HKLFT[:3] = np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T
1422    elif Laue == '3 m 1':
1423        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1424        HKLFT[:3] = np.where((HKLFT[1]+HKLFT[0])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1425        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],matkm[nxs,:,:])).T,HKLFT[:3])
1426    #centrosymmetric
1427    elif Laue == '-3':  #ok
1428        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
1429        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1430        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1431        HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1432        HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[0]<0),-np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1433        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],-mat31[nxs,:,:])).T,HKLFT[:3])   
1434    elif Laue == '-3 m 1':  #ok
1435        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1436        HKLFT[:3] = np.where((HKLFT[1]+HKLFT[0])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1437        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],matkm[nxs,:,:])).T,HKLFT[:3])
1438        HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1439        HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]<HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1440    elif Laue == '-3 1 m':  #ok
1441        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
1442        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1443        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1444        HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1445        HKLFT[:3] = np.where(HKLFT[0]<=0,np.squeeze(np.inner(HKLF[:,:3],-mat31[nxs,:,:])).T,HKLFT[:3])   
1446        HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1447    #hexagonal
1448    #noncentrosymmetric
1449    elif Laue == '6':   #ok
1450        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1451        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1452        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1453        HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1454    elif Laue == '-6':  #ok
1455        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1456        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1457        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1458        HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1459    elif Laue == '6 2 2':   #ok
1460        HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1461        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1462        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1463        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1464        HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1465        HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[0]>HKLFT[1]),np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1466    elif Laue == '6 m m':   #ok
1467        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1468        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1469        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1470        HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1471        HKLFT[:3] = np.where(HKLFT[0]>HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1472    elif Laue == '-6 m 2':  #ok
1473        HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matk2[nxs,:,:])).T,HKLFT[:3])
1474        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1475        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1476        HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1477        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],matk2[nxs,:,:])).T,HKLFT[:3])
1478        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1479    elif Laue == '-6 2 m':  #ok
1480        HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1481        HKLFT[:3] = np.where(HKLFT[0]<=-2*HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1482        HKLFT[:3] = np.where(HKLFT[1]<-2*HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1483        HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1484        HKLFT[:3] = np.where((HKLFT[2]>0)&(HKLFT[1]==HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1485        HKLFT[:3] = np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T
1486        HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1487        HKLFT[:3] = np.where(HKLFT[0]>HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1488    #centrosymmetric
1489    elif Laue == '6/m': #ok
1490        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1491        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1492        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1493        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1494        HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1495    elif Laue == '6/m m m': #ok
1496        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1497        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1498        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1499        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1500        HKLFT[:3] = np.where(HKLFT[0]>HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm.T[nxs,:,:])).T,HKLFT[:3])
1501    #cubic - all ok
1502    #noncentrosymmetric -
1503    elif Laue == '2 3': 
1504        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1505        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1506        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1507        HKLFT[:3] = np.where((HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1508        HKLFT[:3] = np.where((HKLFT[2]>=0)&((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
1509        HKLFT[:3] = np.where((HKLFT[2]>=0)&((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
1510        HKLFT[:3] = np.where((HKLFT[2]<0)&((HKLFT[0]>-HKLFT[2])|(HKLFT[1]>-HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3t[nxs,:,:])).T,HKLFT[:3])
1511        HKLFT[:3] = np.where((HKLFT[2]<0)&((HKLFT[0]>-HKLFT[2])|(HKLFT[1]>=-HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3t[nxs,:,:])).T,HKLFT[:3])
1512        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3])       
1513    elif Laue == '4 3 2':   
1514        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1515        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1516        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1517        HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]<HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1518        HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])   #in lieu od 2-fold
1519        HKLFT[:3] = np.where((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2]),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
1520        HKLFT[:3] = np.where((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2]),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
1521        HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat2d43[nxs,:,:])).T,HKLFT[:3])
1522    elif Laue == '-4 3 m': 
1523        HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])     
1524        HKLFT[:3] = np.where(HKLFT[0]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1525        HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])     
1526        HKLFT[:3] = np.where(HKLFT[1]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1527        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1528        HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1529        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1530        HKLFT[:3] = np.where((HKLFT[2]>=0)&((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
1531        HKLFT[:3] = np.where((HKLFT[2]>=0)&((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
1532        HKLFT[:3] = np.where((HKLFT[2]>=0)&(HKLFT[1]<HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1533        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3]) 
1534        HKLFT[:3] = np.where((HKLFT[0]<0)&(HKLFT[2]<-HKLFT[0])&(HKLFT[1]>HKLFT[2]),np.squeeze(np.inner(HKLF[:,:3],matd3q[nxs,:,:])).T,HKLFT[:3])
1535        HKLFT[:3] = np.where((HKLFT[0]<0)&(HKLFT[2]>=-HKLFT[0])&(HKLFT[1]>HKLFT[2]),np.squeeze(np.inner(HKLF[:,:3],matdm3[nxs,:,:])).T,HKLFT[:3])
1536    #centrosymmetric
1537    elif Laue == 'm 3':
1538        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1539        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1540        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])           
1541        HKLFT[:3] = np.where((HKLFT[2]>=0)&((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
1542        HKLFT[:3] = np.where((HKLFT[2]>=0)&((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
1543    elif Laue == 'm 3 m':
1544        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1545        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1546        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])           
1547        HKLFT[:3] = np.where((HKLFT[2]>=0)&((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
1548        HKLFT[:3] = np.where((HKLFT[2]>=0)&((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
1549        HKLFT[:3] = np.where(HKLFT[0]>HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1550    return HKLFT.T
1551       
1552
1553#Spherical harmonics routines
1554def OdfChk(SGLaue,L,M):
1555    'needs doc string'
1556    if not L%2 and abs(M) <= L:
1557        if SGLaue == '0':                      #cylindrical symmetry
1558            if M == 0: return True
1559        elif SGLaue == '-1':
1560            return True
1561        elif SGLaue == '2/m':
1562            if not abs(M)%2: return True
1563        elif SGLaue == 'mmm':
1564            if not abs(M)%2 and M >= 0: return True
1565        elif SGLaue == '4/m':
1566            if not abs(M)%4: return True
1567        elif SGLaue == '4/mmm':
1568            if not abs(M)%4 and M >= 0: return True
1569        elif SGLaue in ['3R','3']:
1570            if not abs(M)%3: return True
1571        elif SGLaue in ['3mR','3m1','31m']:
1572            if not abs(M)%3 and M >= 0: return True
1573        elif SGLaue == '6/m':
1574            if not abs(M)%6: return True
1575        elif SGLaue == '6/mmm':
1576            if not abs(M)%6 and M >= 0: return True
1577        elif SGLaue == 'm3':
1578            if M > 0:
1579                if L%12 == 2:
1580                    if M <= L//12: return True
1581                else:
1582                    if M <= L//12+1: return True
1583        elif SGLaue == 'm3m':
1584            if M > 0:
1585                if L%12 == 2:
1586                    if M <= L//12: return True
1587                else:
1588                    if M <= L//12+1: return True
1589    return False
1590       
1591def GenSHCoeff(SGLaue,SamSym,L,IfLMN=True):
1592    'needs doc string'
1593    coeffNames = []
1594    for iord in [2*i+2 for i in range(L//2)]:
1595        for m in [i-iord for i in range(2*iord+1)]:
1596            if OdfChk(SamSym,iord,m):
1597                for n in [i-iord for i in range(2*iord+1)]:
1598                    if OdfChk(SGLaue,iord,n):
1599                        if IfLMN:
1600                            coeffNames.append('C(%d,%d,%d)'%(iord,m,n))
1601                        else:
1602                            coeffNames.append('C(%d,%d)'%(iord,n))
1603    return coeffNames
1604   
1605def CrsAng(H,cell,SGData):
1606    'needs doc string'
1607    a,b,c,al,be,ga = cell
1608    SQ3 = 1.732050807569
1609    H1 = np.array([1,0,0])
1610    H2 = np.array([0,1,0])
1611    H3 = np.array([0,0,1])
1612    H4 = np.array([1,1,1])
1613    G,g = cell2Gmat(cell)
1614    Laue = SGData['SGLaue']
1615    Naxis = SGData['SGUniq']
1616    if len(H.shape) == 1:
1617        DH = np.inner(H,np.inner(G,H))
1618    else:
1619        DH = np.array([np.inner(h,np.inner(G,h)) for h in H])
1620    if Laue == '2/m':
1621        if Naxis == 'a':
1622            DR = np.inner(H1,np.inner(G,H1))
1623            DHR = np.inner(H,np.inner(G,H1))
1624        elif Naxis == 'b':
1625            DR = np.inner(H2,np.inner(G,H2))
1626            DHR = np.inner(H,np.inner(G,H2))
1627        else:
1628            DR = np.inner(H3,np.inner(G,H3))
1629            DHR = np.inner(H,np.inner(G,H3))
1630    elif Laue in ['R3','R3m']:
1631        DR = np.inner(H4,np.inner(G,H4))
1632        DHR = np.inner(H,np.inner(G,H4))
1633    else:
1634        DR = np.inner(H3,np.inner(G,H3))
1635        DHR = np.inner(H,np.inner(G,H3))
1636    DHR /= np.sqrt(DR*DH)
1637    phi = np.where(DHR <= 1.0,acosd(DHR),0.0)
1638    if Laue == '-1':
1639        BA = H.T[1]*a/(b-H.T[0]*cosd(ga))
1640        BB = H.T[0]*sind(ga)**2
1641    elif Laue == '2/m':
1642        if Naxis == 'a':
1643            BA = H.T[2]*b/(c-H.T[1]*cosd(al))
1644            BB = H.T[1]*sind(al)**2
1645        elif Naxis == 'b':
1646            BA = H.T[0]*c/(a-H.T[2]*cosd(be))
1647            BB = H.T[2]*sind(be)**2
1648        else:
1649            BA = H.T[1]*a/(b-H.T[0]*cosd(ga))
1650            BB = H.T[0]*sind(ga)**2
1651    elif Laue in ['mmm','4/m','4/mmm']:
1652        BA = H.T[1]*a
1653        BB = H.T[0]*b
1654    elif Laue in ['3R','3mR']:
1655        BA = H.T[0]+H.T[1]-2.0*H.T[2]
1656        BB = SQ3*(H.T[0]-H.T[1])
1657    elif Laue in ['m3','m3m']:
1658        BA = H.T[1]
1659        BB = H.T[0]
1660    else:
1661        BA = H.T[0]+2.0*H.T[1]
1662        BB = SQ3*H.T[0]
1663    beta = atan2d(BA,BB)
1664    return phi,beta
1665   
1666def SamAng(Tth,Gangls,Sangl,IFCoup):
1667    """Compute sample orientation angles vs laboratory coord. system
1668
1669    :param Tth:        Signed theta                                   
1670    :param Gangls:     Sample goniometer angles phi,chi,omega,azmuth 
1671    :param Sangl:      Sample angle zeros om-0, chi-0, phi-0         
1672    :param IFCoup:     True if omega & 2-theta coupled in CW scan
1673    :returns: 
1674        psi,gam:    Sample odf angles                             
1675        dPSdA,dGMdA:    Angle zero derivatives
1676    """                         
1677   
1678    if IFCoup:
1679        GSomeg = sind(Gangls[2]+Tth)
1680        GComeg = cosd(Gangls[2]+Tth)
1681    else:
1682        GSomeg = sind(Gangls[2])
1683        GComeg = cosd(Gangls[2])
1684    GSTth = sind(Tth)
1685    GCTth = cosd(Tth)     
1686    GSazm = sind(Gangls[3])
1687    GCazm = cosd(Gangls[3])
1688    GSchi = sind(Gangls[1])
1689    GCchi = cosd(Gangls[1])
1690    GSphi = sind(Gangls[0]+Sangl[2])
1691    GCphi = cosd(Gangls[0]+Sangl[2])
1692    SSomeg = sind(Sangl[0])
1693    SComeg = cosd(Sangl[0])
1694    SSchi = sind(Sangl[1])
1695    SCchi = cosd(Sangl[1])
1696    AT = -GSTth*GComeg+GCTth*GCazm*GSomeg
1697    BT = GSTth*GSomeg+GCTth*GCazm*GComeg
1698    CT = -GCTth*GSazm*GSchi
1699    DT = -GCTth*GSazm*GCchi
1700   
1701    BC1 = -AT*GSphi+(CT+BT*GCchi)*GCphi
1702    BC2 = DT-BT*GSchi
1703    BC3 = AT*GCphi+(CT+BT*GCchi)*GSphi
1704     
1705    BC = BC1*SComeg*SCchi+BC2*SComeg*SSchi-BC3*SSomeg     
1706    psi = acosd(BC)
1707   
1708    BD = 1.0-BC**2
1709    C = np.where(BD>1.e-6,rpd/np.sqrt(BD),0.)
1710    dPSdA = [-C*(-BC1*SSomeg*SCchi-BC2*SSomeg*SSchi-BC3*SComeg),
1711        -C*(-BC1*SComeg*SSchi+BC2*SComeg*SCchi),
1712        -C*(-BC1*SSomeg-BC3*SComeg*SCchi)]
1713     
1714    BA = -BC1*SSchi+BC2*SCchi
1715    BB = BC1*SSomeg*SCchi+BC2*SSomeg*SSchi+BC3*SComeg
1716    gam = atan2d(BB,BA)
1717
1718    BD = (BA**2+BB**2)/rpd
1719
1720    dBAdO = 0
1721    dBAdC = -BC1*SCchi-BC2*SSchi
1722    dBAdF = BC3*SSchi
1723   
1724    dBBdO = BC1*SComeg*SCchi+BC2*SComeg*SSchi-BC3*SSomeg
1725    dBBdC = -BC1*SSomeg*SSchi+BC2*SSomeg*SCchi
1726    dBBdF = BC1*SComeg-BC3*SSomeg*SCchi
1727   
1728    dGMdA = np.where(BD > 1.e-6,[(BA*dBBdO-BB*dBAdO)/BD,(BA*dBBdC-BB*dBAdC)/BD, \
1729        (BA*dBBdF-BB*dBAdF)/BD],[np.zeros_like(BD),np.zeros_like(BD),np.zeros_like(BD)])
1730       
1731    return psi,gam,dPSdA,dGMdA
1732
1733BOH = {
1734'L=2':[[],[],[]],
1735'L=4':[[0.30469720,0.36418281],[],[]],
1736'L=6':[[-0.14104740,0.52775103],[],[]],
1737'L=8':[[0.28646862,0.21545346,0.32826995],[],[]],
1738'L=10':[[-0.16413497,0.33078546,0.39371345],[],[]],
1739'L=12':[[0.26141975,0.27266871,0.03277460,0.32589402],
1740    [0.09298802,-0.23773812,0.49446631,0.0],[]],
1741'L=14':[[-0.17557309,0.25821932,0.27709173,0.33645360],[],[]],
1742'L=16':[[0.24370673,0.29873515,0.06447688,0.00377,0.32574495],
1743    [0.12039646,-0.25330128,0.23950998,0.40962508,0.0],[]],
1744'L=18':[[-0.16914245,0.17017340,0.34598142,0.07433932,0.32696037],
1745    [-0.06901768,0.16006562,-0.24743528,0.47110273,0.0],[]],
1746'L=20':[[0.23067026,0.31151832,0.09287682,0.01089683,0.00037564,0.32573563],
1747    [0.13615420,-0.25048007,0.12882081,0.28642879,0.34620433,0.0],[]],
1748'L=22':[[-0.16109560,0.10244188,0.36285175,0.13377513,0.01314399,0.32585583],
1749    [-0.09620055,0.20244115,-0.22389483,0.17928946,0.42017231,0.0],[]],
1750'L=24':[[0.22050742,0.31770654,0.11661736,0.02049853,0.00150861,0.00003426,0.32573505],
1751    [0.13651722,-0.21386648,0.00522051,0.33939435,0.10837396,0.32914497,0.0],
1752    [0.05378596,-0.11945819,0.16272298,-0.26449730,0.44923956,0.0,0.0]],
1753'L=26':[[-0.15435003,0.05261630,0.35524646,0.18578869,0.03259103,0.00186197,0.32574594],
1754    [-0.11306511,0.22072681,-0.18706142,0.05439948,0.28122966,0.35634355,0.0],[]],
1755'L=28':[[0.21225019,0.32031716,0.13604702,0.03132468,0.00362703,0.00018294,0.00000294,0.32573501],
1756    [0.13219496,-0.17206256,-0.08742608,0.32671661,0.17973107,0.02567515,0.32619598,0.0],
1757    [0.07989184,-0.16735346,0.18839770,-0.20705337,0.12926808,0.42715602,0.0,0.0]],
1758'L=30':[[-0.14878368,0.01524973,0.33628434,0.22632587,0.05790047,0.00609812,0.00022898,0.32573594],
1759    [-0.11721726,0.20915005,-0.11723436,-0.07815329,0.31318947,0.13655742,0.33241385,0.0],
1760    [-0.04297703,0.09317876,-0.11831248,0.17355132,-0.28164031,0.42719361,0.0,0.0]],
1761'L=32':[[0.20533892,0.32087437,0.15187897,0.04249238,0.00670516,0.00054977,0.00002018,0.00000024,0.32573501],
1762    [0.12775091,-0.13523423,-0.14935701,0.28227378,0.23670434,0.05661270,0.00469819,0.32578978,0.0],
1763    [0.09703829,-0.19373733,0.18610682,-0.14407046,0.00220535,0.26897090,0.36633402,0.0,0.0]],
1764'L=34':[[-0.14409234,-0.01343681,0.31248977,0.25557722,0.08571889,0.01351208,0.00095792,0.00002550,0.32573508],
1765    [-0.11527834,0.18472133,-0.04403280,-0.16908618,0.27227021,0.21086614,0.04041752,0.32688152,0.0],
1766    [-0.06773139,0.14120811,-0.15835721,0.18357456,-0.19364673,0.08377174,0.43116318,0.0,0.0]]
1767}
1768
1769Lnorm = lambda L: 4.*np.pi/(2.0*L+1.)
1770
1771def GetKcl(L,N,SGLaue,phi,beta):
1772    'needs doc string'
1773    import pytexture as ptx
1774    if SGLaue in ['m3','m3m']:
1775        if 'array' in str(type(phi)) and np.any(phi.shape):
1776            Kcl = np.zeros_like(phi)
1777        else:
1778            Kcl = 0.
1779        for j in range(0,L+1,4):
1780            im = j//4
1781            if 'array' in str(type(phi)) and np.any(phi.shape):
1782                pcrs = ptx.pyplmpsi(L,j,len(phi),phi)[0]
1783            else:
1784                pcrs = ptx.pyplmpsi(L,j,1,phi)[0]
1785            Kcl += BOH['L=%d'%(L)][N-1][im]*pcrs*cosd(j*beta)       
1786    else:
1787        if 'array' in str(type(phi)) and np.any(phi.shape):
1788            pcrs = ptx.pyplmpsi(L,N,len(phi),phi)[0]
1789        else:
1790            pcrs = ptx.pyplmpsi(L,N,1,phi)[0]
1791        pcrs *= RSQ2PI
1792        if N:
1793            pcrs *= SQ2
1794        if SGLaue in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']:
1795            if SGLaue in ['3mR','3m1','31m']: 
1796                if N%6 == 3:
1797                    Kcl = pcrs*sind(N*beta)
1798                else:
1799                    Kcl = pcrs*cosd(N*beta)
1800            else:
1801                Kcl = pcrs*cosd(N*beta)
1802        else:
1803            Kcl = pcrs*(cosd(N*beta)+sind(N*beta))
1804    return Kcl
1805   
1806def GetKsl(L,M,SamSym,psi,gam):
1807    'needs doc string'
1808    import pytexture as ptx
1809    if 'array' in str(type(psi)) and np.any(psi.shape):
1810        psrs,dpdps = ptx.pyplmpsi(L,M,len(psi),psi)
1811    else:
1812        psrs,dpdps = ptx.pyplmpsi(L,M,1,psi)
1813    psrs *= RSQ2PI
1814    dpdps *= RSQ2PI
1815    if M:
1816        psrs *= SQ2
1817        dpdps *= SQ2
1818    if SamSym in ['mmm',]:
1819        dum = cosd(M*gam)
1820        Ksl = psrs*dum
1821        dKsdp = dpdps*dum
1822        dKsdg = -psrs*M*sind(M*gam)
1823    else:
1824        dum = cosd(M*gam)+sind(M*gam)
1825        Ksl = psrs*dum
1826        dKsdp = dpdps*dum
1827        dKsdg = psrs*M*(-sind(M*gam)+cosd(M*gam))
1828    return Ksl,dKsdp,dKsdg
1829   
1830def GetKclKsl(L,N,SGLaue,psi,phi,beta):
1831    """
1832    This is used for spherical harmonics description of preferred orientation;
1833        cylindrical symmetry only (M=0) and no sample angle derivatives returned
1834    """
1835    import pytexture as ptx
1836    Ksl,x = ptx.pyplmpsi(L,0,1,psi)
1837    Ksl *= RSQ2PI
1838    if SGLaue in ['m3','m3m']:
1839        Kcl = 0.0
1840        for j in range(0,L+1,4):
1841            im = j//4
1842            pcrs,dum = ptx.pyplmpsi(L,j,1,phi)
1843            Kcl += BOH['L=%d'%(L)][N-1][im]*pcrs*cosd(j*beta)       
1844    else:
1845        pcrs,dum = ptx.pyplmpsi(L,N,1,phi)
1846        pcrs *= RSQ2PI
1847        if N:
1848            pcrs *= SQ2
1849        if SGLaue in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']:
1850            if SGLaue in ['3mR','3m1','31m']: 
1851                if N%6 == 3:
1852                    Kcl = pcrs*sind(N*beta)
1853                else:
1854                    Kcl = pcrs*cosd(N*beta)
1855            else:
1856                Kcl = pcrs*cosd(N*beta)
1857        else:
1858            Kcl = pcrs*(cosd(N*beta)+sind(N*beta))
1859    return Kcl*Ksl,Lnorm(L)
1860   
1861def Glnh(Start,SHCoef,psi,gam,SamSym):
1862    'needs doc string'
1863    import pytexture as ptx
1864
1865    if Start:
1866        ptx.pyqlmninit()
1867        Start = False
1868    Fln = np.zeros(len(SHCoef))
1869    for i,term in enumerate(SHCoef):
1870        l,m,n = eval(term.strip('C'))
1871        pcrs,dum = ptx.pyplmpsi(l,m,1,psi)
1872        pcrs *= RSQPI
1873        if m == 0:
1874            pcrs /= SQ2
1875        if SamSym in ['mmm',]:
1876            Ksl = pcrs*cosd(m*gam)
1877        else:
1878            Ksl = pcrs*(cosd(m*gam)+sind(m*gam))
1879        Fln[i] = SHCoef[term]*Ksl*Lnorm(l)
1880    ODFln = dict(zip(SHCoef.keys(),list(zip(SHCoef.values(),Fln))))
1881    return ODFln
1882
1883def Flnh(Start,SHCoef,phi,beta,SGData):
1884    'needs doc string'
1885    import pytexture as ptx
1886   
1887    if Start:
1888        ptx.pyqlmninit()
1889        Start = False
1890    Fln = np.zeros(len(SHCoef))
1891    for i,term in enumerate(SHCoef):
1892        l,m,n = eval(term.strip('C'))
1893        if SGData['SGLaue'] in ['m3','m3m']:
1894            Kcl = 0.0
1895            for j in range(0,l+1,4):
1896                im = j//4
1897                pcrs,dum = ptx.pyplmpsi(l,j,1,phi)
1898                Kcl += BOH['L='+str(l)][n-1][im]*pcrs*cosd(j*beta)       
1899        else:                #all but cubic
1900            pcrs,dum = ptx.pyplmpsi(l,n,1,phi)
1901            pcrs *= RSQPI
1902            if n == 0:
1903                pcrs /= SQ2
1904            if SGData['SGLaue'] in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']:
1905               if SGData['SGLaue'] in ['3mR','3m1','31m']: 
1906                   if n%6 == 3:
1907                       Kcl = pcrs*sind(n*beta)
1908                   else:
1909                       Kcl = pcrs*cosd(n*beta)
1910               else:
1911                   Kcl = pcrs*cosd(n*beta)
1912            else:
1913                Kcl = pcrs*(cosd(n*beta)+sind(n*beta))
1914        Fln[i] = SHCoef[term]*Kcl*Lnorm(l)
1915    ODFln = dict(zip(SHCoef.keys(),list(zip(SHCoef.values(),Fln))))
1916    return ODFln
1917   
1918def polfcal(ODFln,SamSym,psi,gam):
1919    '''Perform a pole figure computation.
1920    Note that the the number of gam values must either be 1 or must
1921    match psi. Updated for numpy 1.8.0
1922    '''
1923    import pytexture as ptx
1924    PolVal = np.ones_like(psi)
1925    for term in ODFln:
1926        if abs(ODFln[term][1]) > 1.e-3:
1927            l,m,n = eval(term.strip('C'))
1928            psrs,dum = ptx.pyplmpsi(l,m,len(psi),psi)
1929            if SamSym in ['-1','2/m']:
1930                if m:
1931                    Ksl = RSQPI*psrs*(cosd(m*gam)+sind(m*gam))
1932                else:
1933                    Ksl = RSQPI*psrs/SQ2
1934            else:
1935                if m:
1936                    Ksl = RSQPI*psrs*cosd(m*gam)
1937                else:
1938                    Ksl = RSQPI*psrs/SQ2
1939            PolVal += ODFln[term][1]*Ksl
1940    return PolVal
1941   
1942def invpolfcal(ODFln,SGData,phi,beta):
1943    'needs doc string'
1944    import pytexture as ptx
1945   
1946    invPolVal = np.ones_like(beta)
1947    for term in ODFln:
1948        if abs(ODFln[term][1]) > 1.e-3:
1949            l,m,n = eval(term.strip('C'))
1950            if SGData['SGLaue'] in ['m3','m3m']:
1951                Kcl = 0.0
1952                for j in range(0,l+1,4):
1953                    im = j//4
1954                    pcrs,dum = ptx.pyplmpsi(l,j,len(beta),phi)
1955                    Kcl += BOH['L=%d'%(l)][n-1][im]*pcrs*cosd(j*beta)       
1956            else:                #all but cubic
1957                pcrs,dum = ptx.pyplmpsi(l,n,len(beta),phi)
1958                pcrs *= RSQPI
1959                if n == 0:
1960                    pcrs /= SQ2
1961                if SGData['SGLaue'] in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']:
1962                   if SGData['SGLaue'] in ['3mR','3m1','31m']: 
1963                       if n%6 == 3:
1964                           Kcl = pcrs*sind(n*beta)
1965                       else:
1966                           Kcl = pcrs*cosd(n*beta)
1967                   else:
1968                       Kcl = pcrs*cosd(n*beta)
1969                else:
1970                    Kcl = pcrs*(cosd(n*beta)+sind(n*beta))
1971            invPolVal += ODFln[term][1]*Kcl
1972    return invPolVal
1973   
1974   
1975def textureIndex(SHCoef):
1976    'needs doc string'
1977    Tindx = 1.0
1978    for term in SHCoef:
1979        l = eval(term.strip('C'))[0]
1980        Tindx += SHCoef[term]**2/(2.0*l+1.)
1981    return Tindx
1982   
1983# self-test materials follow.
1984selftestlist = []
1985'''Defines a list of self-tests'''
1986selftestquiet = True
1987def _ReportTest():
1988    'Report name and doc string of current routine when ``selftestquiet`` is False'
1989    if not selftestquiet:
1990        import inspect
1991        caller = inspect.stack()[1][3]
1992        doc = eval(caller).__doc__
1993        if doc is not None:
1994            print('testing '+__file__+' with '+caller+' ('+doc+')')
1995        else:
1996            print('testing '+__file__()+" with "+caller)
1997NeedTestData = True
1998def TestData():
1999    array = np.array
2000    global NeedTestData
2001    NeedTestData = False
2002    global CellTestData
2003    # output from uctbx computed on platform darwin on 2010-05-28
2004    CellTestData = [
2005# cell, g, G, cell*, V, V*
2006  [(4, 4, 4, 90, 90, 90), 
2007   array([[  1.60000000e+01,   9.79717439e-16,   9.79717439e-16],
2008       [  9.79717439e-16,   1.60000000e+01,   9.79717439e-16],
2009       [  9.79717439e-16,   9.79717439e-16,   1.60000000e+01]]), array([[  6.25000000e-02,   3.82702125e-18,   3.82702125e-18],
2010       [  3.82702125e-18,   6.25000000e-02,   3.82702125e-18],
2011       [  3.82702125e-18,   3.82702125e-18,   6.25000000e-02]]), (0.25, 0.25, 0.25, 90.0, 90.0, 90.0), 64.0, 0.015625],
2012# cell, g, G, cell*, V, V*
2013  [(4.0999999999999996, 5.2000000000000002, 6.2999999999999998, 100, 80, 130), 
2014   array([[ 16.81      , -13.70423184,   4.48533243],
2015       [-13.70423184,  27.04      ,  -5.6887143 ],
2016       [  4.48533243,  -5.6887143 ,  39.69      ]]), array([[ 0.10206349,  0.05083339, -0.00424823],
2017       [ 0.05083339,  0.06344997,  0.00334956],
2018       [-0.00424823,  0.00334956,  0.02615544]]), (0.31947376387537696, 0.25189277536327803, 0.16172643497798223, 85.283666420376008, 94.716333579624006, 50.825714168082683), 100.98576357983838, 0.0099023858863968445],
2019# cell, g, G, cell*, V, V*
2020  [(3.5, 3.5, 6, 90, 90, 120), 
2021   array([[  1.22500000e+01,  -6.12500000e+00,   1.28587914e-15],
2022       [ -6.12500000e+00,   1.22500000e+01,   1.28587914e-15],
2023       [  1.28587914e-15,   1.28587914e-15,   3.60000000e+01]]), array([[  1.08843537e-01,   5.44217687e-02,   3.36690552e-18],
2024       [  5.44217687e-02,   1.08843537e-01,   3.36690552e-18],
2025       [  3.36690552e-18,   3.36690552e-18,   2.77777778e-02]]), (0.32991443953692895, 0.32991443953692895, 0.16666666666666669, 90.0, 90.0, 60.000000000000021), 63.652867178156257, 0.015710211406520427],
2026  ]
2027    global CoordTestData
2028    CoordTestData = [
2029# cell, ((frac, ortho),...)
2030  ((4,4,4,90,90,90,), [
2031 ((0.10000000000000001, 0.0, 0.0),(0.40000000000000002, 0.0, 0.0)),
2032 ((0.0, 0.10000000000000001, 0.0),(2.4492935982947065e-17, 0.40000000000000002, 0.0)),
2033 ((0.0, 0.0, 0.10000000000000001),(2.4492935982947065e-17, -2.4492935982947065e-17, 0.40000000000000002)),
2034 ((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(0.40000000000000013, 0.79999999999999993, 1.2)),
2035 ((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(0.80000000000000016, 1.2, 0.40000000000000002)),
2036 ((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(1.2, 0.80000000000000004, 0.40000000000000002)),
2037 ((0.5, 0.5, 0.5),(2.0, 1.9999999999999998, 2.0)),
2038]),
2039# cell, ((frac, ortho),...)
2040  ((4.1,5.2,6.3,100,80,130,), [
2041 ((0.10000000000000001, 0.0, 0.0),(0.40999999999999998, 0.0, 0.0)),
2042 ((0.0, 0.10000000000000001, 0.0),(-0.33424955703700043, 0.39834311042186865, 0.0)),
2043 ((0.0, 0.0, 0.10000000000000001),(0.10939835193016617, -0.051013289294572106, 0.6183281045774256)),
2044 ((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(0.069695941716497567, 0.64364635296002093, 1.8549843137322766)),
2045 ((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(-0.073350319180835066, 1.1440160419710339, 0.6183281045774256)),
2046 ((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(0.67089923785616512, 0.74567293154916525, 0.6183281045774256)),
2047 ((0.5, 0.5, 0.5),(0.92574397446582857, 1.7366491056364828, 3.0916405228871278)),
2048]),
2049# cell, ((frac, ortho),...)
2050  ((3.5,3.5,6,90,90,120,), [
2051 ((0.10000000000000001, 0.0, 0.0),(0.35000000000000003, 0.0, 0.0)),
2052 ((0.0, 0.10000000000000001, 0.0),(-0.17499999999999993, 0.3031088913245536, 0.0)),
2053 ((0.0, 0.0, 0.10000000000000001),(3.6739403974420595e-17, -3.6739403974420595e-17, 0.60000000000000009)),
2054 ((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(2.7675166561703527e-16, 0.60621778264910708, 1.7999999999999998)),
2055 ((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(0.17500000000000041, 0.90932667397366063, 0.60000000000000009)),
2056 ((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(0.70000000000000018, 0.6062177826491072, 0.60000000000000009)),
2057 ((0.5, 0.5, 0.5),(0.87500000000000067, 1.5155444566227676, 3.0)),
2058]),
2059]
2060    global LaueTestData             #generated by GSAS
2061    LaueTestData = {
2062    'R 3 m':[(4.,4.,6.,90.,90.,120.),((1,0,1,6),(1,0,-2,6),(0,0,3,2),(1,1,0,6),(2,0,-1,6),(2,0,2,6),
2063        (1,1,3,12),(1,0,4,6),(2,1,1,12),(2,1,-2,12),(3,0,0,6),(1,0,-5,6),(2,0,-4,6),(3,0,-3,6),(3,0,3,6),
2064        (0,0,6,2),(2,2,0,6),(2,1,4,12),(2,0,5,6),(3,1,-1,12),(3,1,2,12),(1,1,6,12),(2,2,3,12),(2,1,-5,12))],
2065    'R 3':[(4.,4.,6.,90.,90.,120.),((1,0,1,6),(1,0,-2,6),(0,0,3,2),(1,1,0,6),(2,0,-1,6),(2,0,2,6),(1,1,3,6),
2066        (1,1,-3,6),(1,0,4,6),(3,-1,1,6),(2,1,1,6),(3,-1,-2,6),(2,1,-2,6),(3,0,0,6),(1,0,-5,6),(2,0,-4,6),
2067        (2,2,0,6),(3,0,3,6),(3,0,-3,6),(0,0,6,2),(3,-1,4,6),(2,0,5,6),(2,1,4,6),(4,-1,-1,6),(3,1,-1,6),
2068        (3,1,2,6),(4,-1,2,6),(2,2,-3,6),(1,1,-6,6),(1,1,6,6),(2,2,3,6),(2,1,-5,6),(3,-1,-5,6))],
2069    'P 3':[(4.,4.,6.,90.,90.,120.),((0,0,1,2),(1,0,0,6),(1,0,1,6),(0,0,2,2),(1,0,-1,6),(1,0,2,6),(1,0,-2,6),
2070        (1,1,0,6),(0,0,3,2),(1,1,1,6),(1,1,-1,6),(1,0,3,6),(1,0,-3,6),(2,0,0,6),(2,0,-1,6),(1,1,-2,6),
2071        (1,1,2,6),(2,0,1,6),(2,0,-2,6),(2,0,2,6),(0,0,4,2),(1,1,-3,6),(1,1,3,6),(1,0,-4,6),(1,0,4,6),
2072        (2,0,-3,6),(2,1,0,6),(2,0,3,6),(3,-1,0,6),(2,1,1,6),(3,-1,-1,6),(2,1,-1,6),(3,-1,1,6),(1,1,4,6),
2073        (3,-1,2,6),(3,-1,-2,6),(1,1,-4,6),(0,0,5,2),(2,1,2,6),(2,1,-2,6),(3,0,0,6),(3,0,1,6),(2,0,4,6),
2074        (2,0,-4,6),(3,0,-1,6),(1,0,-5,6),(1,0,5,6),(3,-1,-3,6),(2,1,-3,6),(2,1,3,6),(3,-1,3,6),(3,0,-2,6),
2075        (3,0,2,6),(1,1,5,6),(1,1,-5,6),(2,2,0,6),(3,0,3,6),(3,0,-3,6),(0,0,6,2),(2,0,-5,6),(2,1,-4,6),
2076        (2,2,-1,6),(3,-1,-4,6),(2,2,1,6),(3,-1,4,6),(2,1,4,6),(2,0,5,6),(1,0,-6,6),(1,0,6,6),(4,-1,0,6),
2077        (3,1,0,6),(3,1,-1,6),(3,1,1,6),(4,-1,-1,6),(2,2,2,6),(4,-1,1,6),(2,2,-2,6),(3,1,2,6),(3,1,-2,6),
2078        (3,0,4,6),(3,0,-4,6),(4,-1,-2,6),(4,-1,2,6),(2,2,-3,6),(1,1,6,6),(1,1,-6,6),(2,2,3,6),(3,-1,5,6),
2079        (2,1,5,6),(2,1,-5,6),(3,-1,-5,6))],
2080    'P 3 m 1':[(4.,4.,6.,90.,90.,120.),((0,0,1,2),(1,0,0,6),(1,0,-1,6),(1,0,1,6),(0,0,2,2),(1,0,-2,6),
2081        (1,0,2,6),(1,1,0,6),(0,0,3,2),(1,1,1,12),(1,0,-3,6),(1,0,3,6),(2,0,0,6),(1,1,2,12),(2,0,1,6),
2082        (2,0,-1,6),(0,0,4,2),(2,0,-2,6),(2,0,2,6),(1,1,3,12),(1,0,-4,6),(1,0,4,6),(2,0,3,6),(2,1,0,12),
2083        (2,0,-3,6),(2,1,1,12),(2,1,-1,12),(1,1,4,12),(2,1,2,12),(0,0,5,2),(2,1,-2,12),(3,0,0,6),(1,0,-5,6),
2084        (3,0,1,6),(3,0,-1,6),(1,0,5,6),(2,0,4,6),(2,0,-4,6),(2,1,3,12),(2,1,-3,12),(3,0,-2,6),(3,0,2,6),
2085        (1,1,5,12),(3,0,-3,6),(0,0,6,2),(2,2,0,6),(3,0,3,6),(2,1,4,12),(2,2,1,12),(2,0,5,6),(2,1,-4,12),
2086        (2,0,-5,6),(1,0,-6,6),(1,0,6,6),(3,1,0,12),(3,1,-1,12),(3,1,1,12),(2,2,2,12),(3,1,2,12),
2087        (3,0,4,6),(3,1,-2,12),(3,0,-4,6),(1,1,6,12),(2,2,3,12))],
2088    'P 3 1 m':[(4.,4.,6.,90.,90.,120.),((0,0,1,2),(1,0,0,6),(0,0,2,2),(1,0,1,12),(1,0,2,12),(1,1,0,6),
2089        (0,0,3,2),(1,1,-1,6),(1,1,1,6),(1,0,3,12),(2,0,0,6),(2,0,1,12),(1,1,2,6),(1,1,-2,6),(2,0,2,12),
2090        (0,0,4,2),(1,1,-3,6),(1,1,3,6),(1,0,4,12),(2,1,0,12),(2,0,3,12),(2,1,1,12),(2,1,-1,12),(1,1,-4,6),
2091        (1,1,4,6),(0,0,5,2),(2,1,-2,12),(2,1,2,12),(3,0,0,6),(1,0,5,12),(2,0,4,12),(3,0,1,12),(2,1,-3,12),
2092        (2,1,3,12),(3,0,2,12),(1,1,5,6),(1,1,-5,6),(3,0,3,12),(0,0,6,2),(2,2,0,6),(2,1,-4,12),(2,0,5,12),
2093        (2,2,-1,6),(2,2,1,6),(2,1,4,12),(3,1,0,12),(1,0,6,12),(2,2,2,6),(3,1,-1,12),(2,2,-2,6),(3,1,1,12),
2094        (3,1,-2,12),(3,0,4,12),(3,1,2,12),(1,1,-6,6),(2,2,3,6),(2,2,-3,6),(1,1,6,6))],
2095    }
2096   
2097    global FLnhTestData
2098    FLnhTestData = [{
2099    'C(4,0,0)': (0.965, 0.42760447),
2100    'C(2,0,0)': (1.0122, -0.80233610),
2101    'C(2,0,2)': (0.0061, 8.37491546E-03),
2102    'C(6,0,4)': (-0.0898, 4.37985696E-02),
2103    'C(6,0,6)': (-0.1369, -9.04081762E-02),
2104    'C(6,0,0)': (0.5935, -0.18234928),
2105    'C(4,0,4)': (0.1872, 0.16358127),
2106    'C(6,0,2)': (0.6193, 0.27573633),
2107    'C(4,0,2)': (-0.1897, 0.12530720)},[1,0,0]]
2108def test0():
2109    if NeedTestData: TestData()
2110    msg = 'test cell2Gmat, fillgmat, Gmat2cell'
2111    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2112        G, g = cell2Gmat(cell)
2113        assert np.allclose(G,tG),msg
2114        assert np.allclose(g,tg),msg
2115        tcell = Gmat2cell(g)
2116        assert np.allclose(cell,tcell),msg
2117        tcell = Gmat2cell(G)
2118        assert np.allclose(tcell,trcell),msg
2119if __name__ == '__main__': selftestlist.append(test0)
2120
2121def test1():
2122    'test cell2A and A2Gmat'
2123    _ReportTest()
2124    if NeedTestData: TestData()
2125    msg = 'test cell2A and A2Gmat'
2126    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2127        G, g = A2Gmat(cell2A(cell))
2128        assert np.allclose(G,tG),msg
2129        assert np.allclose(g,tg),msg
2130if __name__ == '__main__': selftestlist.append(test1)
2131
2132def test2():
2133    'test Gmat2A, A2cell, A2Gmat, Gmat2cell'
2134    _ReportTest()
2135    if NeedTestData: TestData()
2136    msg = 'test Gmat2A, A2cell, A2Gmat, Gmat2cell'
2137    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2138        G, g = cell2Gmat(cell)
2139        tcell = A2cell(Gmat2A(G))
2140        assert np.allclose(cell,tcell),msg
2141if __name__ == '__main__': selftestlist.append(test2)
2142
2143def test3():
2144    'test invcell2Gmat'
2145    _ReportTest()
2146    if NeedTestData: TestData()
2147    msg = 'test invcell2Gmat'
2148    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2149        G, g = invcell2Gmat(trcell)
2150        assert np.allclose(G,tG),msg
2151        assert np.allclose(g,tg),msg
2152if __name__ == '__main__': selftestlist.append(test3)
2153
2154def test4():
2155    'test calc_rVsq, calc_rV, calc_V'
2156    _ReportTest()
2157    if NeedTestData: TestData()
2158    msg = 'test calc_rVsq, calc_rV, calc_V'
2159    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2160        assert np.allclose(calc_rV(cell2A(cell)),trV), msg
2161        assert np.allclose(calc_V(cell2A(cell)),tV), msg
2162if __name__ == '__main__': selftestlist.append(test4)
2163
2164def test5():
2165    'test A2invcell'
2166    _ReportTest()
2167    if NeedTestData: TestData()
2168    msg = 'test A2invcell'
2169    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2170        rcell = A2invcell(cell2A(cell))
2171        assert np.allclose(rcell,trcell),msg
2172if __name__ == '__main__': selftestlist.append(test5)
2173
2174def test6():
2175    'test cell2AB'
2176    _ReportTest()
2177    if NeedTestData: TestData()
2178    msg = 'test cell2AB'
2179    for (cell,coordlist) in CoordTestData:
2180        A,B = cell2AB(cell)
2181        for (frac,ortho) in coordlist:
2182            to = np.inner(A,frac)
2183            tf = np.inner(B,to)
2184            assert np.allclose(ortho,to), msg
2185            assert np.allclose(frac,tf), msg
2186            to = np.sum(A*frac,axis=1)
2187            tf = np.sum(B*to,axis=1)
2188            assert np.allclose(ortho,to), msg
2189            assert np.allclose(frac,tf), msg
2190if __name__ == '__main__': selftestlist.append(test6)
2191
2192def test7():
2193    'test GetBraviasNum(...) and GenHBravais(...)'
2194    _ReportTest()
2195    import os.path
2196    import sys
2197    import GSASIIspc as spc
2198    testdir = os.path.join(os.path.split(os.path.abspath( __file__ ))[0],'testinp')
2199    if os.path.exists(testdir):
2200        if testdir not in sys.path: sys.path.insert(0,testdir)
2201    import sgtbxlattinp
2202    derror = 1e-4
2203    def indexmatch(hklin, hkllist, system):
2204        for hklref in hkllist:
2205            hklref = list(hklref)
2206            # these permutations are far from complete, but are sufficient to
2207            # allow the test to complete
2208            if system == 'cubic':
2209                permlist = [(1,2,3),(1,3,2),(2,1,3),(2,3,1),(3,1,2),(3,2,1),]
2210            elif system == 'monoclinic':
2211                permlist = [(1,2,3),(-1,2,-3)]
2212            else:
2213                permlist = [(1,2,3)]
2214
2215            for perm in permlist:
2216                hkl = [abs(i) * hklin[abs(i)-1] / i for i in perm]
2217                if hkl == hklref: return True
2218                if [-i for i in hkl] == hklref: return True
2219        else:
2220            return False
2221
2222    for key in sgtbxlattinp.sgtbx7:
2223        spdict = spc.SpcGroup(key)
2224        cell = sgtbxlattinp.sgtbx7[key][0]
2225        system = spdict[1]['SGSys']
2226        center = spdict[1]['SGLatt']
2227
2228        bravcode = GetBraviasNum(center, system)
2229
2230        g2list = GenHBravais(sgtbxlattinp.dmin, bravcode, cell2A(cell))
2231
2232        assert len(sgtbxlattinp.sgtbx7[key][1]) == len(g2list), 'Reflection lists differ for %s' % key
2233        for h,k,l,d,num in g2list:
2234            for hkllist,dref in sgtbxlattinp.sgtbx7[key][1]: 
2235                if abs(d-dref) < derror:
2236                    if indexmatch((h,k,l,), hkllist, system):
2237                        break
2238            else:
2239                assert 0,'No match for %s at %s (%s)' % ((h,k,l),d,key)
2240if __name__ == '__main__': selftestlist.append(test7)
2241
2242def test8():
2243    'test GenHLaue'
2244    _ReportTest()
2245    import GSASIIspc as spc
2246    import sgtbxlattinp
2247    derror = 1e-4
2248    dmin = sgtbxlattinp.dmin
2249
2250    def indexmatch(hklin, hklref, system, axis):
2251        # these permutations are far from complete, but are sufficient to
2252        # allow the test to complete
2253        if system == 'cubic':
2254            permlist = [(1,2,3),(1,3,2),(2,1,3),(2,3,1),(3,1,2),(3,2,1),]
2255        elif system == 'monoclinic' and axis=='b':
2256            permlist = [(1,2,3),(-1,2,-3)]
2257        elif system == 'monoclinic' and axis=='a':
2258            permlist = [(1,2,3),(1,-2,-3)]
2259        elif system == 'monoclinic' and axis=='c':
2260            permlist = [(1,2,3),(-1,-2,3)]
2261        elif system == 'trigonal':
2262            permlist = [(1,2,3),(2,1,3),(-1,-2,3),(-2,-1,3)]
2263        elif system == 'rhombohedral':
2264            permlist = [(1,2,3),(2,3,1),(3,1,2)]
2265        else:
2266            permlist = [(1,2,3)]
2267
2268        hklref = list(hklref)
2269        for perm in permlist:
2270            hkl = [abs(i) * hklin[abs(i)-1] / i for i in perm]
2271            if hkl == hklref: return True
2272            if [-i for i in hkl] == hklref: return True
2273        return False
2274
2275    for key in sgtbxlattinp.sgtbx8:
2276        spdict = spc.SpcGroup(key)[1]
2277        cell = sgtbxlattinp.sgtbx8[key][0]
2278        Axis = spdict['SGUniq']
2279        system = spdict['SGSys']
2280
2281        g2list = GenHLaue(dmin,spdict,cell2A(cell))
2282        #if len(g2list) != len(sgtbxlattinp.sgtbx8[key][1]):
2283        #    print 'failed',key,':' ,len(g2list),'vs',len(sgtbxlattinp.sgtbx8[key][1])
2284        #    print 'GSAS-II:'
2285        #    for h,k,l,d in g2list: print '  ',(h,k,l),d
2286        #    print 'SGTBX:'
2287        #    for hkllist,dref in sgtbxlattinp.sgtbx8[key][1]: print '  ',hkllist,dref
2288        assert len(g2list) == len(sgtbxlattinp.sgtbx8[key][1]), (
2289            'Reflection lists differ for %s' % key
2290            )
2291        #match = True
2292        for h,k,l,d in g2list:
2293            for hkllist,dref in sgtbxlattinp.sgtbx8[key][1]: 
2294                if abs(d-dref) < derror:
2295                    if indexmatch((h,k,l,), hkllist, system, Axis): break
2296            else:
2297                assert 0,'No match for %s at %s (%s)' % ((h,k,l),d,key)
2298                #match = False
2299        #if not match:
2300            #for hkllist,dref in sgtbxlattinp.sgtbx8[key][1]: print '  ',hkllist,dref
2301            #print center, Laue, Axis, system
2302if __name__ == '__main__': selftestlist.append(test8)
2303           
2304def test9():
2305    'test GenHLaue'
2306    _ReportTest()
2307    import GSASIIspc as G2spc
2308    if NeedTestData: TestData()
2309    for spc in LaueTestData:
2310        data = LaueTestData[spc]
2311        cell = data[0]
2312        hklm = np.array(data[1])
2313        H = hklm[-1][:3]
2314        hklO = hklm.T[:3].T
2315        A = cell2A(cell)
2316        dmin = 1./np.sqrt(calc_rDsq(H,A))
2317        SGData = G2spc.SpcGroup(spc)[1]
2318        hkls = np.array(GenHLaue(dmin,SGData,A))
2319        hklN = hkls.T[:3].T
2320        #print spc,hklO.shape,hklN.shape
2321        err = True
2322        for H in hklO:
2323            if H not in hklN:
2324                print ('%d %s'%(H,' missing from hkl from GSASII'))
2325                err = False
2326        assert(err)
2327if __name__ == '__main__': selftestlist.append(test9)
2328       
2329       
2330   
2331
2332if __name__ == '__main__':
2333    # run self-tests
2334    selftestquiet = False
2335    for test in selftestlist:
2336        test()
2337    print ("OK")
Note: See TracBrowser for help on using the repository browser.