source: trunk/GSASIIlattice.py @ 3596

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

fix OnTransform? for magnetic phases to give correct atom sites/symmetry

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