source: trunk/GSASIIlattice.py @ 3638

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

implement a count of unique magnetic atoms for each result from k-SUBGROUPSMAG, display it & add filter on how many are allowed
Modify RefreshKeep? to allow change of filters for magnetic space groups
change map Zstep size limit to 4A from 1A

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