source: trunk/GSASIIlattice.py @ 3587

Last change on this file since 3587 was 3587, checked in by vondreele, 5 years ago

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