source: trunk/GSASIIlattice.py @ 3406

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

fix Transform to create correct(?) constraints between lattices, atom positions & anisotropic thermal parameters.
Checked for axes permutations in orthorhombic.

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