source: trunk/GSASIIlattice.py @ 3685

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

Add progress bar for processing k-SUBGROUPSMAG results - it does take a while

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