source: trunk/GSASIIlattice.py @ 3586

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

rename MiltiFloatDialog? to MultiDataDialog? - better now that it does more than floats
add choice as an option to MultiDataDialog? - makes a ComboBox? for selection of an item
add optional choice of magnetic atom to test for possible moment as part of making table of mag space groups

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