source: trunk/GSASIIlattice.py @ 3918

Last change on this file since 3918 was 3888, checked in by vondreele, 6 years ago

provide contour plot (not finished yet) on structure drawing

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