source: trunk/GSASIIlattice.py @ 3633

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

extend Unit Cells List Bravais lattice choices to include Ammm and Bmmm. This allows test indexing to find 100, 010 or 001 propagation vectors. Also allows peak indexing in these Bravais lattices in addition to Cmmm.

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