source: trunk/GSASIIlattice.py @ 3583

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

fixes so orthorhombic mag space groups from Bilbao conserve axes - nonstd space groups produced

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