source: trunk/GSASIIlattice.py @ 3103

Last change on this file since 3103 was 3103, checked in by vondreele, 4 years ago

show seq. result in structure drawing - steps thru histograms changing structure at each step

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