source: trunk/GSASIIlattice.py @ 3591

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

work on FindNonstandard? for mag space groups
enhance UseMagAtomDialog? to show mag site symmetry

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