source: trunk/GSASIIlattice.py @ 4175

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

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