source: trunk/GSASIIlattice.py @ 3621

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

Modify Transform tool to allow setting spins & displaying operators
Transform operation matches mag space group selection from Bilbao
select nonstandard orthorhombic works (except P 2221 & P 2'2'21) - get same set of mag atom site symmetries

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