source: trunk/GSASIIlattice.py @ 3654

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

fix to bug if no atom types selected for mag atom check
add code to transform monoclinics to nicer unit cells - currently disabled

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