source: trunk/GSASIIlattice.py @ 3646

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

modify UseMagAtomDialog? to have ifOK option - only OK button G2phsGUI)
change Unit Cells List mag phase tables to respond to 1)space group - give operator table (no longer row number) & 2)Uniq - give unique mag atom positions (just OK button)
make TestMagAtoms? work over entire mag cell volume to get all unique mag atoms - Uniq in table now correct
G2lat - new routine ExpandCell? - fills mag unit cell (before transformation)

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