source: trunk/GSASIIlattice.py @ 3801

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

combine cell2AB & Gmat2AB in G2lat
fixes to mag. str. fctr. math. Correct sym. op. transformation of moments & fix moment to cartesian axis problem
  • 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-01-27 15:04:12 +0000 (Sun, 27 Jan 2019) $
27# $Author: vondreele $
28# $Revision: 3801 $
29# $URL: trunk/GSASIIlattice.py $
30# $Id: GSASIIlattice.py 3801 2019-01-27 15:04:12Z vondreele $
31########### SVN repository information ###################
32from __future__ import division, print_function
33import math
34import copy
35import sys
36import random as ran
37import numpy as np
38import numpy.linalg as nl
39import GSASIIpath
40import GSASIImath as G2mth
41import GSASIIspc as G2spc
42import GSASIIElem as G2elem
43GSASIIpath.SetVersionNumber("$Revision: 3801 $")
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,Bmat,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#            GSASIIpath.IPyBreak()
927            Stack.append(HX)
928    return Stack
929       
930def MaxIndex(dmin,A):
931    'needs doc string'
932    Hmax = [0,0,0]
933    try:
934        cell = A2cell(A)
935    except:
936        cell = [1.,1.,1.,90.,90.,90.]
937    for i in range(3):
938        Hmax[i] = int(round(cell[i]/dmin))
939    return Hmax
940   
941def transposeHKLF(transMat,Super,refList):
942    ''' Apply transformation matrix to hkl(m)
943    param: transmat: 3x3 or 4x4 array
944    param: Super: 0 or 1 for extra index
945    param: refList list of h,k,l,....
946    return: newRefs transformed list of h',k',l',,,
947    return: badRefs list of noninteger h',k',l'...
948    '''
949    newRefs = np.copy(refList)
950    badRefs = []
951    for H in newRefs:
952        newH = np.inner(transMat,H[:3+Super])
953        H[:3+Super] = np.rint(newH)
954        if not np.allclose(newH,H[:3+Super],atol=0.01):
955            badRefs.append(newH)
956    return newRefs,badRefs
957   
958def sortHKLd(HKLd,ifreverse,ifdup,ifSS=False):
959    '''sort reflection list on d-spacing; can sort in either order
960
961    :param HKLd: a list of [h,k,l,d,...];
962    :param ifreverse: True for largest d first
963    :param ifdup: True if duplicate d-spacings allowed
964    :return: sorted reflection list
965    '''
966    T = []
967    N = 3
968    if ifSS:
969        N = 4
970    for i,H in enumerate(HKLd):
971        if ifdup:
972            T.append((H[N],i))
973        else:
974            T.append(H[N])           
975    D = dict(zip(T,HKLd))
976    T.sort()
977    if ifreverse:
978        T.reverse()
979    X = []
980    okey = ''
981    for key in T: 
982        if key != okey: X.append(D[key])    #remove duplicate d-spacings
983        okey = key
984    return X
985   
986def SwapIndx(Axis,H):
987    'needs doc string'
988    if Axis in [1,-1]:
989        return H
990    elif Axis in [2,-3]:
991        return [H[1],H[2],H[0]]
992    else:
993        return [H[2],H[0],H[1]]
994       
995def Rh2Hx(Rh):
996    'needs doc string'
997    Hx = [0,0,0]
998    Hx[0] = Rh[0]-Rh[1]
999    Hx[1] = Rh[1]-Rh[2]
1000    Hx[2] = np.sum(Rh)
1001    return Hx
1002   
1003def Hx2Rh(Hx):
1004    'needs doc string'
1005    Rh = [0,0,0]
1006    itk = -Hx[0]+Hx[1]+Hx[2]
1007    if itk%3 != 0:
1008        return 0        #error - not rhombohedral reflection
1009    else:
1010        Rh[1] = itk//3
1011        Rh[0] = Rh[1]+Hx[0]
1012        Rh[2] = Rh[1]-Hx[1]
1013        if Rh[0] < 0:
1014            for i in range(3):
1015                Rh[i] = -Rh[i]
1016        return Rh
1017       
1018def CentCheck(Cent,H):
1019    'needs doc string'
1020    h,k,l = H
1021    if Cent == 'A' and (k+l)%2:
1022        return False
1023    elif Cent == 'B' and (h+l)%2:
1024        return False
1025    elif Cent == 'C' and (h+k)%2:
1026        return False
1027    elif Cent == 'I' and (h+k+l)%2:
1028        return False
1029    elif Cent == 'F' and ((h+k)%2 or (h+l)%2 or (k+l)%2):
1030        return False
1031    elif Cent == 'R' and (-h+k+l)%3:
1032        return False
1033    else:
1034        return True
1035                                   
1036def GetBraviasNum(center,system):
1037    """Determine the Bravais lattice number, as used in GenHBravais
1038   
1039    :param center: one of: 'P', 'C', 'I', 'F', 'R' (see SGLatt from GSASIIspc.SpcGroup)
1040    :param system: one of 'cubic', 'hexagonal', 'tetragonal', 'orthorhombic', 'trigonal' (for R)
1041      'monoclinic', 'triclinic' (see SGSys from GSASIIspc.SpcGroup)
1042    :return: a number between 0 and 13
1043      or throws a ValueError exception if the combination of center, system is not found (i.e. non-standard)
1044
1045    """
1046    if center.upper() == 'F' and system.lower() == 'cubic':
1047        return 0
1048    elif center.upper() == 'I' and system.lower() == 'cubic':
1049        return 1
1050    elif center.upper() == 'P' and system.lower() == 'cubic':
1051        return 2
1052    elif center.upper() == 'R' and system.lower() == 'trigonal':
1053        return 3
1054    elif center.upper() == 'P' and system.lower() == 'hexagonal':
1055        return 4
1056    elif center.upper() == 'I' and system.lower() == 'tetragonal':
1057        return 5
1058    elif center.upper() == 'P' and system.lower() == 'tetragonal':
1059        return 6
1060    elif center.upper() == 'F' and system.lower() == 'orthorhombic':
1061        return 7
1062    elif center.upper() == 'I' and system.lower() == 'orthorhombic':
1063        return 8
1064    elif center.upper() == 'A' and system.lower() == 'orthorhombic':
1065        return 9
1066    elif center.upper() == 'B' and system.lower() == 'orthorhombic':
1067        return 10
1068    elif center.upper() == 'C' and system.lower() == 'orthorhombic':
1069        return 11
1070    elif center.upper() == 'P' and system.lower() == 'orthorhombic':
1071        return 12
1072    elif center.upper() == 'C' and system.lower() == 'monoclinic':
1073        return 13
1074    elif center.upper() == 'P' and system.lower() == 'monoclinic':
1075        return 14
1076    elif center.upper() == 'P' and system.lower() == 'triclinic':
1077        return 15
1078    raise ValueError('non-standard Bravais lattice center=%s, cell=%s' % (center,system))
1079
1080def GenHBravais(dmin,Bravais,A):
1081    """Generate the positionally unique powder diffraction reflections
1082     
1083    :param dmin: minimum d-spacing in A
1084    :param Bravais: lattice type (see GetBraviasNum). Bravais is one of:
1085   
1086            * 0 F cubic
1087            * 1 I cubic
1088            * 2 P cubic
1089            * 3 R hexagonal (trigonal not rhombohedral)
1090            * 4 P hexagonal
1091            * 5 I tetragonal
1092            * 6 P tetragonal
1093            * 7 F orthorhombic
1094            * 8 I orthorhombic
1095            * 9 A orthorhombic
1096            * 10 B orthorhombic
1097            * 11 C orthorhombic
1098            * 12 P orthorhombic
1099            * 13 I monoclinic
1100            * 14 C monoclinic
1101            * 15 P monoclinic
1102            * 16 P triclinic
1103           
1104    :param A: reciprocal metric tensor elements as [G11,G22,G33,2*G12,2*G13,2*G23]
1105    :return: HKL unique d list of [h,k,l,d,-1] sorted with largest d first
1106           
1107    """
1108    if Bravais in [9,]:
1109        Cent = 'A'
1110    elif Bravais in [10,]:
1111        Cent = 'B'
1112    elif Bravais in [11,14]:
1113        Cent = 'C'
1114    elif Bravais in [1,5,8,13]:
1115        Cent = 'I'
1116    elif Bravais in [0,7]:
1117        Cent = 'F'
1118    elif Bravais in [3]:
1119        Cent = 'R'
1120    else:
1121        Cent = 'P'
1122    Hmax = MaxIndex(dmin,A)
1123    dminsq = 1./(dmin**2)
1124    HKL = []
1125    if Bravais == 16:                       #triclinic
1126        for l in range(-Hmax[2],Hmax[2]+1):
1127            for k in range(-Hmax[1],Hmax[1]+1):
1128                hmin = 0
1129                if (k < 0): hmin = 1
1130                if (k ==0 and l < 0): hmin = 1
1131                for h in range(hmin,Hmax[0]+1):
1132                    H=[h,k,l]
1133                    rdsq = calc_rDsq(H,A)
1134                    if 0 < rdsq <= dminsq:
1135                        HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
1136    elif Bravais in [13,14,15]:                #monoclinic - b unique
1137        Hmax = SwapIndx(2,Hmax)
1138        for h in range(Hmax[0]+1):
1139            for k in range(-Hmax[1],Hmax[1]+1):
1140                lmin = 0
1141                if k < 0:lmin = 1
1142                for l in range(lmin,Hmax[2]+1):
1143                    [h,k,l] = SwapIndx(-2,[h,k,l])
1144                    H = []
1145                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
1146                    if H:
1147                        rdsq = calc_rDsq(H,A)
1148                        if 0 < rdsq <= dminsq:
1149                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
1150                    [h,k,l] = SwapIndx(2,[h,k,l])
1151    elif Bravais in [7,8,9,10,11,12]:            #orthorhombic
1152        for h in range(Hmax[0]+1):
1153            for k in range(Hmax[1]+1):
1154                for l in range(Hmax[2]+1):
1155                    H = []
1156                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
1157                    if H:
1158                        rdsq = calc_rDsq(H,A)
1159                        if 0 < rdsq <= dminsq:
1160                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
1161    elif Bravais in [5,6]:                  #tetragonal
1162        for l in range(Hmax[2]+1):
1163            for k in range(Hmax[1]+1):
1164                for h in range(k,Hmax[0]+1):
1165                    H = []
1166                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
1167                    if H:
1168                        rdsq = calc_rDsq(H,A)
1169                        if 0 < rdsq <= dminsq:
1170                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
1171    elif Bravais in [3,4]:
1172        lmin = 0
1173        if Bravais == 3: lmin = -Hmax[2]                  #hexagonal/trigonal
1174        for l in range(lmin,Hmax[2]+1):
1175            for k in range(Hmax[1]+1):
1176                hmin = k
1177                if l < 0: hmin += 1
1178                for h in range(hmin,Hmax[0]+1):
1179                    H = []
1180                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
1181                    if H:
1182                        rdsq = calc_rDsq(H,A)
1183                        if 0 < rdsq <= dminsq:
1184                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
1185
1186    else:                                   #cubic
1187        for l in range(Hmax[2]+1):
1188            for k in range(l,Hmax[1]+1):
1189                for h in range(k,Hmax[0]+1):
1190                    H = []
1191                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
1192                    if H:
1193                        rdsq = calc_rDsq(H,A)
1194                        if 0 < rdsq <= dminsq:
1195                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
1196    return sortHKLd(HKL,True,False)
1197   
1198def getHKLmax(dmin,SGData,A):
1199    'finds maximum allowed hkl for given A within dmin'
1200    SGLaue = SGData['SGLaue']
1201    if SGLaue in ['3R','3mR']:        #Rhombohedral axes
1202        Hmax = [0,0,0]
1203        cell = A2cell(A)
1204        aHx = cell[0]*math.sqrt(2.0*(1.0-cosd(cell[3])))
1205        cHx = cell[0]*math.sqrt(3.0*(1.0+2.0*cosd(cell[3])))
1206        Hmax[0] = Hmax[1] = int(round(aHx/dmin))
1207        Hmax[2] = int(round(cHx/dmin))
1208        #print Hmax,aHx,cHx
1209    else:                           # all others
1210        Hmax = MaxIndex(dmin,A)
1211    return Hmax
1212   
1213def GenHLaue(dmin,SGData,A):
1214    """Generate the crystallographically unique powder diffraction reflections
1215    for a lattice and Bravais type
1216   
1217    :param dmin: minimum d-spacing
1218    :param SGData: space group dictionary with at least
1219   
1220        * 'SGLaue': Laue group symbol: one of '-1','2/m','mmm','4/m','6/m','4/mmm','6/mmm', '3m1', '31m', '3', '3R', '3mR', 'm3', 'm3m'
1221        * 'SGLatt': lattice centering: one of 'P','A','B','C','I','F'
1222        * 'SGUniq': code for unique monoclinic axis one of 'a','b','c' (only if 'SGLaue' is '2/m') otherwise an empty string
1223       
1224    :param A: reciprocal metric tensor elements as [G11,G22,G33,2*G12,2*G13,2*G23]
1225    :return: HKL = list of [h,k,l,d] sorted with largest d first and is unique
1226            part of reciprocal space ignoring anomalous dispersion
1227           
1228    """
1229    import math
1230    SGLaue = SGData['SGLaue']
1231    SGLatt = SGData['SGLatt']
1232    SGUniq = SGData['SGUniq']
1233    #finds maximum allowed hkl for given A within dmin
1234    Hmax = getHKLmax(dmin,SGData,A)
1235       
1236    dminsq = 1./(dmin**2)
1237    HKL = []
1238    if SGLaue == '-1':                       #triclinic
1239        for l in range(-Hmax[2],Hmax[2]+1):
1240            for k in range(-Hmax[1],Hmax[1]+1):
1241                hmin = 0
1242                if (k < 0) or (k ==0 and l < 0): hmin = 1
1243                for h in range(hmin,Hmax[0]+1):
1244                    H = []
1245                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
1246                    if H:
1247                        rdsq = calc_rDsq(H,A)
1248                        if 0 < rdsq <= dminsq:
1249                            HKL.append([h,k,l,1./math.sqrt(rdsq)])
1250    elif SGLaue == '2/m':                #monoclinic
1251        axisnum = 1 + ['a','b','c'].index(SGUniq)
1252        Hmax = SwapIndx(axisnum,Hmax)
1253        for h in range(Hmax[0]+1):
1254            for k in range(-Hmax[1],Hmax[1]+1):
1255                lmin = 0
1256                if k < 0:lmin = 1
1257                for l in range(lmin,Hmax[2]+1):
1258                    [h,k,l] = SwapIndx(-axisnum,[h,k,l])
1259                    H = []
1260                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
1261                    if H:
1262                        rdsq = calc_rDsq(H,A)
1263                        if 0 < rdsq <= dminsq:
1264                            HKL.append([h,k,l,1./math.sqrt(rdsq)])
1265                    [h,k,l] = SwapIndx(axisnum,[h,k,l])
1266    elif SGLaue in ['mmm','4/m','6/m']:            #orthorhombic
1267        for l in range(Hmax[2]+1):
1268            for h in range(Hmax[0]+1):
1269                kmin = 1
1270                if SGLaue == 'mmm' or h ==0: kmin = 0
1271                for k in range(kmin,Hmax[1]+1):
1272                    H = []
1273                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
1274                    if H:
1275                        rdsq = calc_rDsq(H,A)
1276                        if 0 < rdsq <= dminsq:
1277                            HKL.append([h,k,l,1./math.sqrt(rdsq)])
1278    elif SGLaue in ['4/mmm','6/mmm']:                  #tetragonal & hexagonal
1279        for l in range(Hmax[2]+1):
1280            for h in range(Hmax[0]+1):
1281                for k in range(h+1):
1282                    H = []
1283                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
1284                    if H:
1285                        rdsq = calc_rDsq(H,A)
1286                        if 0 < rdsq <= dminsq:
1287                            HKL.append([h,k,l,1./math.sqrt(rdsq)])
1288    elif SGLaue in ['3m1','31m','3','3R','3mR']:                  #trigonals
1289        for l in range(-Hmax[2],Hmax[2]+1):
1290            hmin = 0
1291            if l < 0: hmin = 1
1292            for h in range(hmin,Hmax[0]+1):
1293                if SGLaue in ['3R','3']:
1294                    kmax = h
1295                    kmin = -int((h-1.)/2.)
1296                else:
1297                    kmin = 0
1298                    kmax = h
1299                    if SGLaue in ['3m1','3mR'] and l < 0: kmax = h-1
1300                    if SGLaue == '31m' and l < 0: kmin = 1
1301                for k in range(kmin,kmax+1):
1302                    H = []
1303                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
1304                    if SGLaue in ['3R','3mR']:
1305                        H = Hx2Rh(H)
1306                    if H:
1307                        rdsq = calc_rDsq(H,A)
1308                        if 0 < rdsq <= dminsq:
1309                            HKL.append([H[0],H[1],H[2],1./math.sqrt(rdsq)])
1310    else:                                   #cubic
1311        for h in range(Hmax[0]+1):
1312            for k in range(h+1):
1313                lmin = 0
1314                lmax = k
1315                if SGLaue =='m3':
1316                    lmax = h-1
1317                    if h == k: lmax += 1
1318                for l in range(lmin,lmax+1):
1319                    H = []
1320                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
1321                    if H:
1322                        rdsq = calc_rDsq(H,A)
1323                        if 0 < rdsq <= dminsq:
1324                            HKL.append([h,k,l,1./math.sqrt(rdsq)])
1325    return sortHKLd(HKL,True,True)
1326   
1327def GenPfHKLs(nMax,SGData,A):   
1328    """Generate the unique pole figure reflections for a lattice and Bravais type.
1329    Min d-spacing=1.0A & no more than nMax returned
1330   
1331    :param nMax: maximum number of hkls returned
1332    :param SGData: space group dictionary with at least
1333   
1334        * 'SGLaue': Laue group symbol: one of '-1','2/m','mmm','4/m','6/m','4/mmm','6/mmm', '3m1', '31m', '3', '3R', '3mR', 'm3', 'm3m'
1335        * 'SGLatt': lattice centering: one of 'P','A','B','C','I','F'
1336        * 'SGUniq': code for unique monoclinic axis one of 'a','b','c' (only if 'SGLaue' is '2/m') otherwise an empty string
1337       
1338    :param A: reciprocal metric tensor elements as [G11,G22,G33,2*G12,2*G13,2*G23]
1339    :return: HKL = list of 'h k l' strings sorted with largest d first; no duplicate zones
1340           
1341    """
1342    HKL = np.array(GenHLaue(1.0,SGData,A)).T[:3].T     #strip d-spacings
1343    N = min(nMax,len(HKL))
1344    return ['%d %d %d'%(h[0],h[1],h[2]) for h in HKL[:N]]       
1345
1346def GenSSHLaue(dmin,SGData,SSGData,Vec,maxH,A):
1347    'needs a doc string'
1348    ifMag = False
1349    if 'MagSpGrp' in SGData:
1350        ifMag = True
1351    HKLs = []
1352    vec = np.array(Vec)
1353    vstar = np.sqrt(calc_rDsq(vec,A))     #find extra needed for -n SS reflections
1354    dvec = 1./(maxH*vstar+1./dmin)
1355    HKL = GenHLaue(dvec,SGData,A)       
1356    SSdH = [vec*h for h in range(-maxH,maxH+1)]
1357    SSdH = dict(zip(range(-maxH,maxH+1),SSdH))
1358    for h,k,l,d in HKL:
1359        ext = G2spc.GenHKLf([h,k,l],SGData)[0]  #h,k,l must be integral values here
1360        if not ext and d >= dmin:
1361            HKLs.append([h,k,l,0,d])
1362        for dH in SSdH:
1363            if dH:
1364                DH = SSdH[dH]
1365                H = [h+DH[0],k+DH[1],l+DH[2]]
1366                d = 1./np.sqrt(calc_rDsq(H,A))
1367                if d >= dmin:
1368                    HKLM = np.array([h,k,l,dH])
1369                    if (G2spc.checkSSLaue([h,k,l,dH],SGData,SSGData) and G2spc.checkSSextc(HKLM,SSGData)) or ifMag:
1370                        HKLs.append([h,k,l,dH,d])   
1371    return HKLs
1372   
1373def LaueUnique2(SGData,refList):
1374    ''' Impose Laue symmetry on hkl
1375   
1376    :param SGData: space group data from 'P '+Laue
1377    :param HKLF: np.array([[h,k,l,...]]) reflection set to be converted
1378   
1379    :return: HKLF new reflection array with imposed Laue symmetry
1380    '''
1381    for ref in refList:
1382        H = ref[:3]
1383        Uniq = G2spc.GenHKLf(H,SGData)[2]
1384        Uniq = G2mth.sortArray(G2mth.sortArray(G2mth.sortArray(Uniq,2),1),0)
1385        ref[:3] = Uniq[-1]
1386    return refList
1387   
1388def LaueUnique(Laue,HKLF):
1389    ''' Impose Laue symmetry on hkl
1390   
1391    :param str Laue: Laue symbol, as below
1392   
1393      centrosymmetric Laue groups::
1394       
1395            ['-1','2/m','112/m','2/m11','mmm','-42m','-4m2','4/mmm','-3',
1396            '-31m','-3m1','6/m','6/mmm','m3','m3m']
1397     
1398      noncentrosymmetric Laue groups::
1399     
1400           ['1','2','211','112','m','m11','11m','222','mm2','m2m','2mm',
1401           '4','-4','422','4mm','3','312','321','31m','3m1','6','-6',
1402           '622','6mm','-62m','-6m2','23','432','-43m']
1403     
1404    :param HKLF: np.array([[h,k,l,...]]) reflection set to be converted
1405   
1406    :returns: HKLF new reflection array with imposed Laue symmetry
1407    '''
1408   
1409    HKLFT = HKLF.T
1410    mat41 = np.array([[0,1,0],[-1,0,0],[0,0,1]])    #hkl -> k,-h,l
1411    mat43 = np.array([[0,-1,0],[1,0,0],[0,0,1]])    #hkl -> -k,h,l
1412    mat4bar = np.array([[0,-1,0],[1,0,0],[0,0,-1]]) #hkl -> k,-h,-l
1413    mat31 = np.array([[-1,-1,0],[1,0,0],[0,0,1]])   #hkl -> ihl = -h-k,h,l
1414    mat32 = np.array([[0,1,0],[-1,-1,0],[0,0,1]])   #hkl -> kil = k,-h-k,l
1415    matd3 = np.array([[0,1,0],[0,0,1],[1,0,0]])     #hkl -> k,l,h
1416    matd3q = np.array([[0,0,-1],[-1,0,0],[0,1,0]])  #hkl -> -l,-h,k
1417    matd3t = np.array([[0,0,-1],[1,0,0],[0,-1,0]])  #hkl -> -l,h,-k
1418    mat6 = np.array([[1,1,0],[-1,0,0],[0,0,1]])     #hkl -> h+k,-h,l really 65
1419    matdm = np.array([[0,1,0],[1,0,0],[0,0,1]])     #hkl -> k,h,l
1420    matdmp = np.array([[-1,-1,0],[0,1,0],[0,0,1]])  #hkl -> -h-k,k,l
1421    matkm = np.array([[-1,0,0],[1,1,0],[0,0,1]])    #hkl -> -h,h+k,l
1422    matd2 = np.array([[0,1,0],[1,0,0],[0,0,-1]])    #hkl -> k,h,-l
1423    matdm3 = np.array([[1,0,0],[0,0,1],[0,1,0]])    #hkl -> h,l,k
1424    mat2d43 = np.array([[0,1,0],[1,0,0],[0,0,1]])   #hkl -> k,-h,l
1425    matk2 = np.array([[-1,0,0],[1,1,0],[0,0,-1]])   #hkl -> -h,-i,-l
1426    #triclinic
1427    if Laue == '1': #ok
1428        pass
1429    elif Laue == '-1':  #ok
1430        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
1431        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]<0),HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
1432        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
1433    #monoclinic
1434    #noncentrosymmetric - all ok
1435    elif Laue == '2': 
1436        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3])
1437        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3])
1438    elif Laue == '1 1 2':
1439        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1440        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]<0),HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1441    elif Laue == '2 1 1':   
1442        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1443        HKLFT[:3] = np.where((HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1444    elif Laue == 'm':
1445        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1446    elif Laue == 'm 1 1':
1447        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1448    elif Laue == '1 1 m':
1449        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1450    #centrosymmetric - all ok
1451    elif Laue == '2/m 1 1':       
1452        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1453        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1454        HKLFT[:3] = np.where((HKLFT[2]*HKLFT[0]==0)&(HKLFT[1]<0),HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1455    elif Laue == '2/m':
1456        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3])
1457        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1458        HKLFT[:3] = np.where((HKLFT[0]*HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1459    elif Laue == '1 1 2/m':
1460        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1461        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1462        HKLFT[:3] = np.where((HKLFT[1]*HKLFT[2]==0)&(HKLFT[0]<0),HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1463    #orthorhombic
1464    #noncentrosymmetric - all OK
1465    elif Laue == '2 2 2':
1466        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1467        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1468        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1469        HKLFT[:3] = np.where((HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1470    elif Laue == 'm m 2':
1471        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1472        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1473    elif Laue == '2 m m': 
1474        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1475        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1476    elif Laue == 'm 2 m':
1477        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1478        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1479    #centrosymmetric - all ok
1480    elif Laue == 'm m m':
1481        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1482        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1483        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1484    #tetragonal
1485    #noncentrosymmetric - all ok
1486    elif Laue == '4':
1487        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1488        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1489        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]>0),np.squeeze(np.inner(HKLF[:,:3],mat41[nxs,:,:])).T,HKLFT[:3])
1490    elif Laue == '-4': 
1491        HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])     
1492        HKLFT[:3] = np.where(HKLFT[0]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1493        HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])     
1494        HKLFT[:3] = np.where(HKLFT[1]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1495        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1496    elif Laue == '4 2 2':
1497        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1498        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1499        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1500        HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]<HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1501        HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])   #in lieu od 2-fold
1502    elif Laue == '4 m m':
1503        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1504        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1505        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1506        HKLFT[:3] = np.where(HKLFT[0]<HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1507    elif Laue == '-4 2 m':
1508        HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])     
1509        HKLFT[:3] = np.where(HKLFT[0]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1510        HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])     
1511        HKLFT[:3] = np.where(HKLFT[1]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1512        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1513        HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1514        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1515    elif Laue == '-4 m 2':
1516        HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1517        HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])     
1518        HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]<=0),np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1519        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]<0),HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])     
1520        HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]==0),np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1521        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3]) 
1522        HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[0]>HKLFT[1]),np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1523    #centrosymmetric - all ok
1524    elif Laue == '4/m':
1525        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1526        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1527        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1528        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]>0),np.squeeze(np.inner(HKLF[:,:3],mat41[nxs,:,:])).T,HKLFT[:3])
1529    elif Laue == '4/m m m':
1530        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1531        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1532        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])       
1533        HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],mat41[nxs,:,:])).T,HKLFT[:3])
1534        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1535    #trigonal - all hex cell
1536    #noncentrosymmetric - all ok
1537    elif Laue == '3':
1538        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1539        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1540        HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1541    elif Laue == '3 1 2':
1542        HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matk2[nxs,:,:])).T,HKLFT[:3])
1543        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1544        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1545        HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1546        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],matk2[nxs,:,:])).T,HKLFT[:3])
1547    elif Laue == '3 2 1':
1548        HKLFT[:3] = np.where(HKLFT[0]<=-2*HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1549        HKLFT[:3] = np.where(HKLFT[1]<-2*HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1550        HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1551        HKLFT[:3] = np.where((HKLFT[2]>0)&(HKLFT[1]==HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1552        HKLFT[:3] = np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T
1553        HKLFT[:3] = np.where((HKLFT[0]!=0)&(HKLFT[2]>0)&(HKLFT[0]==-2*HKLFT[1]),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1554    elif Laue == '3 1 m':
1555        HKLFT[:3] = np.where(HKLFT[0]>=HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1556        HKLFT[:3] = np.where(2*HKLFT[1]<-HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1557        HKLFT[:3] = np.where(HKLFT[1]>-2*HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matdmp[nxs,:,:])).T,HKLFT[:3])
1558        HKLFT[:3] = np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T
1559    elif Laue == '3 m 1':
1560        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1561        HKLFT[:3] = np.where((HKLFT[1]+HKLFT[0])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1562        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],matkm[nxs,:,:])).T,HKLFT[:3])
1563    #centrosymmetric
1564    elif Laue == '-3':  #ok
1565        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
1566        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1567        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1568        HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1569        HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[0]<0),-np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1570        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],-mat31[nxs,:,:])).T,HKLFT[:3])   
1571    elif Laue == '-3 m 1':  #ok
1572        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1573        HKLFT[:3] = np.where((HKLFT[1]+HKLFT[0])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1574        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],matkm[nxs,:,:])).T,HKLFT[:3])
1575        HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1576        HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]<HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1577    elif Laue == '-3 1 m':  #ok
1578        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
1579        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1580        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1581        HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1582        HKLFT[:3] = np.where(HKLFT[0]<=0,np.squeeze(np.inner(HKLF[:,:3],-mat31[nxs,:,:])).T,HKLFT[:3])   
1583        HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1584    #hexagonal
1585    #noncentrosymmetric
1586    elif Laue == '6':   #ok
1587        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1588        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1589        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1590        HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1591    elif Laue == '-6':  #ok
1592        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1593        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1594        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1595        HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1596    elif Laue == '6 2 2':   #ok
1597        HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1598        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1599        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1600        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1601        HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1602        HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[0]>HKLFT[1]),np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1603    elif Laue == '6 m m':   #ok
1604        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1605        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1606        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1607        HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1608        HKLFT[:3] = np.where(HKLFT[0]>HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1609    elif Laue == '-6 m 2':  #ok
1610        HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matk2[nxs,:,:])).T,HKLFT[:3])
1611        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1612        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1613        HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1614        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],matk2[nxs,:,:])).T,HKLFT[:3])
1615        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1616    elif Laue == '-6 2 m':  #ok
1617        HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1618        HKLFT[:3] = np.where(HKLFT[0]<=-2*HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1619        HKLFT[:3] = np.where(HKLFT[1]<-2*HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1620        HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1621        HKLFT[:3] = np.where((HKLFT[2]>0)&(HKLFT[1]==HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1622        HKLFT[:3] = np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T
1623        HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1624        HKLFT[:3] = np.where(HKLFT[0]>HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1625    #centrosymmetric
1626    elif Laue == '6/m': #ok
1627        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1628        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1629        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1630        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1631        HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1632    elif Laue == '6/m m m': #ok
1633        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1634        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1635        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1636        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1637        HKLFT[:3] = np.where(HKLFT[0]>HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm.T[nxs,:,:])).T,HKLFT[:3])
1638    #cubic - all ok
1639    #noncentrosymmetric -
1640    elif Laue == '2 3': 
1641        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1642        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1643        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1644        HKLFT[:3] = np.where((HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1645        HKLFT[:3] = np.where((HKLFT[2]>=0)&((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
1646        HKLFT[:3] = np.where((HKLFT[2]>=0)&((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
1647        HKLFT[:3] = np.where((HKLFT[2]<0)&((HKLFT[0]>-HKLFT[2])|(HKLFT[1]>-HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3t[nxs,:,:])).T,HKLFT[:3])
1648        HKLFT[:3] = np.where((HKLFT[2]<0)&((HKLFT[0]>-HKLFT[2])|(HKLFT[1]>=-HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3t[nxs,:,:])).T,HKLFT[:3])
1649        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3])       
1650    elif Laue == '4 3 2':   
1651        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1652        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1653        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1654        HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]<HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1655        HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])   #in lieu od 2-fold
1656        HKLFT[:3] = np.where((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2]),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
1657        HKLFT[:3] = np.where((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2]),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
1658        HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat2d43[nxs,:,:])).T,HKLFT[:3])
1659    elif Laue == '-4 3 m': 
1660        HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])     
1661        HKLFT[:3] = np.where(HKLFT[0]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1662        HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])     
1663        HKLFT[:3] = np.where(HKLFT[1]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1664        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1665        HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1666        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1667        HKLFT[:3] = np.where((HKLFT[2]>=0)&((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
1668        HKLFT[:3] = np.where((HKLFT[2]>=0)&((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
1669        HKLFT[:3] = np.where((HKLFT[2]>=0)&(HKLFT[1]<HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1670        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3]) 
1671        HKLFT[:3] = np.where((HKLFT[0]<0)&(HKLFT[2]<-HKLFT[0])&(HKLFT[1]>HKLFT[2]),np.squeeze(np.inner(HKLF[:,:3],matd3q[nxs,:,:])).T,HKLFT[:3])
1672        HKLFT[:3] = np.where((HKLFT[0]<0)&(HKLFT[2]>=-HKLFT[0])&(HKLFT[1]>HKLFT[2]),np.squeeze(np.inner(HKLF[:,:3],matdm3[nxs,:,:])).T,HKLFT[:3])
1673    #centrosymmetric
1674    elif Laue == 'm 3':
1675        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1676        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1677        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])           
1678        HKLFT[:3] = np.where((HKLFT[2]>=0)&((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
1679        HKLFT[:3] = np.where((HKLFT[2]>=0)&((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
1680    elif Laue == 'm 3 m':
1681        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1682        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1683        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])           
1684        HKLFT[:3] = np.where((HKLFT[2]>=0)&((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
1685        HKLFT[:3] = np.where((HKLFT[2]>=0)&((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
1686        HKLFT[:3] = np.where(HKLFT[0]>HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1687    return HKLFT.T
1688       
1689
1690#Spherical harmonics routines
1691def OdfChk(SGLaue,L,M):
1692    'needs doc string'
1693    if not L%2 and abs(M) <= L:
1694        if SGLaue == '0':                      #cylindrical symmetry
1695            if M == 0: return True
1696        elif SGLaue == '-1':
1697            return True
1698        elif SGLaue == '2/m':
1699            if not abs(M)%2: return True
1700        elif SGLaue == 'mmm':
1701            if not abs(M)%2 and M >= 0: return True
1702        elif SGLaue == '4/m':
1703            if not abs(M)%4: return True
1704        elif SGLaue == '4/mmm':
1705            if not abs(M)%4 and M >= 0: return True
1706        elif SGLaue in ['3R','3']:
1707            if not abs(M)%3: return True
1708        elif SGLaue in ['3mR','3m1','31m']:
1709            if not abs(M)%3 and M >= 0: return True
1710        elif SGLaue == '6/m':
1711            if not abs(M)%6: return True
1712        elif SGLaue == '6/mmm':
1713            if not abs(M)%6 and M >= 0: return True
1714        elif SGLaue == 'm3':
1715            if M > 0:
1716                if L%12 == 2:
1717                    if M <= L//12: return True
1718                else:
1719                    if M <= L//12+1: return True
1720        elif SGLaue == 'm3m':
1721            if M > 0:
1722                if L%12 == 2:
1723                    if M <= L//12: return True
1724                else:
1725                    if M <= L//12+1: return True
1726    return False
1727       
1728def GenSHCoeff(SGLaue,SamSym,L,IfLMN=True):
1729    'needs doc string'
1730    coeffNames = []
1731    for iord in [2*i+2 for i in range(L//2)]:
1732        for m in [i-iord for i in range(2*iord+1)]:
1733            if OdfChk(SamSym,iord,m):
1734                for n in [i-iord for i in range(2*iord+1)]:
1735                    if OdfChk(SGLaue,iord,n):
1736                        if IfLMN:
1737                            coeffNames.append('C(%d,%d,%d)'%(iord,m,n))
1738                        else:
1739                            coeffNames.append('C(%d,%d)'%(iord,n))
1740    return coeffNames
1741   
1742def CrsAng(H,cell,SGData):
1743    'needs doc string'
1744    a,b,c,al,be,ga = cell
1745    SQ3 = 1.732050807569
1746    H1 = np.array([1,0,0])
1747    H2 = np.array([0,1,0])
1748    H3 = np.array([0,0,1])
1749    H4 = np.array([1,1,1])
1750    G,g = cell2Gmat(cell)
1751    Laue = SGData['SGLaue']
1752    Naxis = SGData['SGUniq']
1753    if len(H.shape) == 1:
1754        DH = np.inner(H,np.inner(G,H))
1755    else:
1756        DH = np.array([np.inner(h,np.inner(G,h)) for h in H])
1757    if Laue == '2/m':
1758        if Naxis == 'a':
1759            DR = np.inner(H1,np.inner(G,H1))
1760            DHR = np.inner(H,np.inner(G,H1))
1761        elif Naxis == 'b':
1762            DR = np.inner(H2,np.inner(G,H2))
1763            DHR = np.inner(H,np.inner(G,H2))
1764        else:
1765            DR = np.inner(H3,np.inner(G,H3))
1766            DHR = np.inner(H,np.inner(G,H3))
1767    elif Laue in ['R3','R3m']:
1768        DR = np.inner(H4,np.inner(G,H4))
1769        DHR = np.inner(H,np.inner(G,H4))
1770    else:
1771        DR = np.inner(H3,np.inner(G,H3))
1772        DHR = np.inner(H,np.inner(G,H3))
1773    DHR /= np.sqrt(DR*DH)
1774    phi = np.where(DHR <= 1.0,acosd(DHR),0.0)
1775    if Laue == '-1':
1776        BA = H.T[1]*a/(b-H.T[0]*cosd(ga))
1777        BB = H.T[0]*sind(ga)**2
1778    elif Laue == '2/m':
1779        if Naxis == 'a':
1780            BA = H.T[2]*b/(c-H.T[1]*cosd(al))
1781            BB = H.T[1]*sind(al)**2
1782        elif Naxis == 'b':
1783            BA = H.T[0]*c/(a-H.T[2]*cosd(be))
1784            BB = H.T[2]*sind(be)**2
1785        else:
1786            BA = H.T[1]*a/(b-H.T[0]*cosd(ga))
1787            BB = H.T[0]*sind(ga)**2
1788    elif Laue in ['mmm','4/m','4/mmm']:
1789        BA = H.T[1]*a
1790        BB = H.T[0]*b
1791    elif Laue in ['3R','3mR']:
1792        BA = H.T[0]+H.T[1]-2.0*H.T[2]
1793        BB = SQ3*(H.T[0]-H.T[1])
1794    elif Laue in ['m3','m3m']:
1795        BA = H.T[1]
1796        BB = H.T[0]
1797    else:
1798        BA = H.T[0]+2.0*H.T[1]
1799        BB = SQ3*H.T[0]
1800    beta = atan2d(BA,BB)
1801    return phi,beta
1802   
1803def SamAng(Tth,Gangls,Sangl,IFCoup):
1804    """Compute sample orientation angles vs laboratory coord. system
1805
1806    :param Tth:        Signed theta                                   
1807    :param Gangls:     Sample goniometer angles phi,chi,omega,azmuth 
1808    :param Sangl:      Sample angle zeros om-0, chi-0, phi-0         
1809    :param IFCoup:     True if omega & 2-theta coupled in CW scan
1810    :returns: 
1811        psi,gam:    Sample odf angles                             
1812        dPSdA,dGMdA:    Angle zero derivatives
1813    """                         
1814   
1815    if IFCoup:
1816        GSomeg = sind(Gangls[2]+Tth)
1817        GComeg = cosd(Gangls[2]+Tth)
1818    else:
1819        GSomeg = sind(Gangls[2])
1820        GComeg = cosd(Gangls[2])
1821    GSTth = sind(Tth)
1822    GCTth = cosd(Tth)     
1823    GSazm = sind(Gangls[3])
1824    GCazm = cosd(Gangls[3])
1825    GSchi = sind(Gangls[1])
1826    GCchi = cosd(Gangls[1])
1827    GSphi = sind(Gangls[0]+Sangl[2])
1828    GCphi = cosd(Gangls[0]+Sangl[2])
1829    SSomeg = sind(Sangl[0])
1830    SComeg = cosd(Sangl[0])
1831    SSchi = sind(Sangl[1])
1832    SCchi = cosd(Sangl[1])
1833    AT = -GSTth*GComeg+GCTth*GCazm*GSomeg
1834    BT = GSTth*GSomeg+GCTth*GCazm*GComeg
1835    CT = -GCTth*GSazm*GSchi
1836    DT = -GCTth*GSazm*GCchi
1837   
1838    BC1 = -AT*GSphi+(CT+BT*GCchi)*GCphi
1839    BC2 = DT-BT*GSchi
1840    BC3 = AT*GCphi+(CT+BT*GCchi)*GSphi
1841     
1842    BC = BC1*SComeg*SCchi+BC2*SComeg*SSchi-BC3*SSomeg     
1843    psi = acosd(BC)
1844   
1845    BD = 1.0-BC**2
1846    C = np.where(BD>1.e-6,rpd/np.sqrt(BD),0.)
1847    dPSdA = [-C*(-BC1*SSomeg*SCchi-BC2*SSomeg*SSchi-BC3*SComeg),
1848        -C*(-BC1*SComeg*SSchi+BC2*SComeg*SCchi),
1849        -C*(-BC1*SSomeg-BC3*SComeg*SCchi)]
1850     
1851    BA = -BC1*SSchi+BC2*SCchi
1852    BB = BC1*SSomeg*SCchi+BC2*SSomeg*SSchi+BC3*SComeg
1853    gam = atan2d(BB,BA)
1854
1855    BD = (BA**2+BB**2)/rpd
1856
1857    dBAdO = 0
1858    dBAdC = -BC1*SCchi-BC2*SSchi
1859    dBAdF = BC3*SSchi
1860   
1861    dBBdO = BC1*SComeg*SCchi+BC2*SComeg*SSchi-BC3*SSomeg
1862    dBBdC = -BC1*SSomeg*SSchi+BC2*SSomeg*SCchi
1863    dBBdF = BC1*SComeg-BC3*SSomeg*SCchi
1864   
1865    dGMdA = np.where(BD > 1.e-6,[(BA*dBBdO-BB*dBAdO)/BD,(BA*dBBdC-BB*dBAdC)/BD, \
1866        (BA*dBBdF-BB*dBAdF)/BD],[np.zeros_like(BD),np.zeros_like(BD),np.zeros_like(BD)])
1867       
1868    return psi,gam,dPSdA,dGMdA
1869
1870BOH = {
1871'L=2':[[],[],[]],
1872'L=4':[[0.30469720,0.36418281],[],[]],
1873'L=6':[[-0.14104740,0.52775103],[],[]],
1874'L=8':[[0.28646862,0.21545346,0.32826995],[],[]],
1875'L=10':[[-0.16413497,0.33078546,0.39371345],[],[]],
1876'L=12':[[0.26141975,0.27266871,0.03277460,0.32589402],
1877    [0.09298802,-0.23773812,0.49446631,0.0],[]],
1878'L=14':[[-0.17557309,0.25821932,0.27709173,0.33645360],[],[]],
1879'L=16':[[0.24370673,0.29873515,0.06447688,0.00377,0.32574495],
1880    [0.12039646,-0.25330128,0.23950998,0.40962508,0.0],[]],
1881'L=18':[[-0.16914245,0.17017340,0.34598142,0.07433932,0.32696037],
1882    [-0.06901768,0.16006562,-0.24743528,0.47110273,0.0],[]],
1883'L=20':[[0.23067026,0.31151832,0.09287682,0.01089683,0.00037564,0.32573563],
1884    [0.13615420,-0.25048007,0.12882081,0.28642879,0.34620433,0.0],[]],
1885'L=22':[[-0.16109560,0.10244188,0.36285175,0.13377513,0.01314399,0.32585583],
1886    [-0.09620055,0.20244115,-0.22389483,0.17928946,0.42017231,0.0],[]],
1887'L=24':[[0.22050742,0.31770654,0.11661736,0.02049853,0.00150861,0.00003426,0.32573505],
1888    [0.13651722,-0.21386648,0.00522051,0.33939435,0.10837396,0.32914497,0.0],
1889    [0.05378596,-0.11945819,0.16272298,-0.26449730,0.44923956,0.0,0.0]],
1890'L=26':[[-0.15435003,0.05261630,0.35524646,0.18578869,0.03259103,0.00186197,0.32574594],
1891    [-0.11306511,0.22072681,-0.18706142,0.05439948,0.28122966,0.35634355,0.0],[]],
1892'L=28':[[0.21225019,0.32031716,0.13604702,0.03132468,0.00362703,0.00018294,0.00000294,0.32573501],
1893    [0.13219496,-0.17206256,-0.08742608,0.32671661,0.17973107,0.02567515,0.32619598,0.0],
1894    [0.07989184,-0.16735346,0.18839770,-0.20705337,0.12926808,0.42715602,0.0,0.0]],
1895'L=30':[[-0.14878368,0.01524973,0.33628434,0.22632587,0.05790047,0.00609812,0.00022898,0.32573594],
1896    [-0.11721726,0.20915005,-0.11723436,-0.07815329,0.31318947,0.13655742,0.33241385,0.0],
1897    [-0.04297703,0.09317876,-0.11831248,0.17355132,-0.28164031,0.42719361,0.0,0.0]],
1898'L=32':[[0.20533892,0.32087437,0.15187897,0.04249238,0.00670516,0.00054977,0.00002018,0.00000024,0.32573501],
1899    [0.12775091,-0.13523423,-0.14935701,0.28227378,0.23670434,0.05661270,0.00469819,0.32578978,0.0],
1900    [0.09703829,-0.19373733,0.18610682,-0.14407046,0.00220535,0.26897090,0.36633402,0.0,0.0]],
1901'L=34':[[-0.14409234,-0.01343681,0.31248977,0.25557722,0.08571889,0.01351208,0.00095792,0.00002550,0.32573508],
1902    [-0.11527834,0.18472133,-0.04403280,-0.16908618,0.27227021,0.21086614,0.04041752,0.32688152,0.0],
1903    [-0.06773139,0.14120811,-0.15835721,0.18357456,-0.19364673,0.08377174,0.43116318,0.0,0.0]]
1904}
1905
1906Lnorm = lambda L: 4.*np.pi/(2.0*L+1.)
1907
1908def GetKcl(L,N,SGLaue,phi,beta):
1909    'needs doc string'
1910    import pytexture as ptx
1911    if SGLaue in ['m3','m3m']:
1912        if 'array' in str(type(phi)) and np.any(phi.shape):
1913            Kcl = np.zeros_like(phi)
1914        else:
1915            Kcl = 0.
1916        for j in range(0,L+1,4):
1917            im = j//4
1918            if 'array' in str(type(phi)) and np.any(phi.shape):
1919                pcrs = ptx.pyplmpsi(L,j,len(phi),phi)[0]
1920            else:
1921                pcrs = ptx.pyplmpsi(L,j,1,phi)[0]
1922            Kcl += BOH['L=%d'%(L)][N-1][im]*pcrs*cosd(j*beta)       
1923    else:
1924        if 'array' in str(type(phi)) and np.any(phi.shape):
1925            pcrs = ptx.pyplmpsi(L,N,len(phi),phi)[0]
1926        else:
1927            pcrs = ptx.pyplmpsi(L,N,1,phi)[0]
1928        pcrs *= RSQ2PI
1929        if N:
1930            pcrs *= SQ2
1931        if SGLaue in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']:
1932            if SGLaue in ['3mR','3m1','31m']: 
1933                if N%6 == 3:
1934                    Kcl = pcrs*sind(N*beta)
1935                else:
1936                    Kcl = pcrs*cosd(N*beta)
1937            else:
1938                Kcl = pcrs*cosd(N*beta)
1939        else:
1940            Kcl = pcrs*(cosd(N*beta)+sind(N*beta))
1941    return Kcl
1942   
1943def GetKsl(L,M,SamSym,psi,gam):
1944    'needs doc string'
1945    import pytexture as ptx
1946    if 'array' in str(type(psi)) and np.any(psi.shape):
1947        psrs,dpdps = ptx.pyplmpsi(L,M,len(psi),psi)
1948    else:
1949        psrs,dpdps = ptx.pyplmpsi(L,M,1,psi)
1950    psrs *= RSQ2PI
1951    dpdps *= RSQ2PI
1952    if M:
1953        psrs *= SQ2
1954        dpdps *= SQ2
1955    if SamSym in ['mmm',]:
1956        dum = cosd(M*gam)
1957        Ksl = psrs*dum
1958        dKsdp = dpdps*dum
1959        dKsdg = -psrs*M*sind(M*gam)
1960    else:
1961        dum = cosd(M*gam)+sind(M*gam)
1962        Ksl = psrs*dum
1963        dKsdp = dpdps*dum
1964        dKsdg = psrs*M*(-sind(M*gam)+cosd(M*gam))
1965    return Ksl,dKsdp,dKsdg
1966   
1967def GetKclKsl(L,N,SGLaue,psi,phi,beta):
1968    """
1969    This is used for spherical harmonics description of preferred orientation;
1970        cylindrical symmetry only (M=0) and no sample angle derivatives returned
1971    """
1972    import pytexture as ptx
1973    Ksl,x = ptx.pyplmpsi(L,0,1,psi)
1974    Ksl *= RSQ2PI
1975    if SGLaue in ['m3','m3m']:
1976        Kcl = 0.0
1977        for j in range(0,L+1,4):
1978            im = j//4
1979            pcrs,dum = ptx.pyplmpsi(L,j,1,phi)
1980            Kcl += BOH['L=%d'%(L)][N-1][im]*pcrs*cosd(j*beta)       
1981    else:
1982        pcrs,dum = ptx.pyplmpsi(L,N,1,phi)
1983        pcrs *= RSQ2PI
1984        if N:
1985            pcrs *= SQ2
1986        if SGLaue in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']:
1987            if SGLaue in ['3mR','3m1','31m']: 
1988                if N%6 == 3:
1989                    Kcl = pcrs*sind(N*beta)
1990                else:
1991                    Kcl = pcrs*cosd(N*beta)
1992            else:
1993                Kcl = pcrs*cosd(N*beta)
1994        else:
1995            Kcl = pcrs*(cosd(N*beta)+sind(N*beta))
1996    return Kcl*Ksl,Lnorm(L)
1997   
1998def Glnh(Start,SHCoef,psi,gam,SamSym):
1999    'needs doc string'
2000    import pytexture as ptx
2001
2002    if Start:
2003        ptx.pyqlmninit()
2004        Start = False
2005    Fln = np.zeros(len(SHCoef))
2006    for i,term in enumerate(SHCoef):
2007        l,m,n = eval(term.strip('C'))
2008        pcrs,dum = ptx.pyplmpsi(l,m,1,psi)
2009        pcrs *= RSQPI
2010        if m == 0:
2011            pcrs /= SQ2
2012        if SamSym in ['mmm',]:
2013            Ksl = pcrs*cosd(m*gam)
2014        else:
2015            Ksl = pcrs*(cosd(m*gam)+sind(m*gam))
2016        Fln[i] = SHCoef[term]*Ksl*Lnorm(l)
2017    ODFln = dict(zip(SHCoef.keys(),list(zip(SHCoef.values(),Fln))))
2018    return ODFln
2019
2020def Flnh(Start,SHCoef,phi,beta,SGData):
2021    'needs doc string'
2022    import pytexture as ptx
2023   
2024    if Start:
2025        ptx.pyqlmninit()
2026        Start = False
2027    Fln = np.zeros(len(SHCoef))
2028    for i,term in enumerate(SHCoef):
2029        l,m,n = eval(term.strip('C'))
2030        if SGData['SGLaue'] in ['m3','m3m']:
2031            Kcl = 0.0
2032            for j in range(0,l+1,4):
2033                im = j//4
2034                pcrs,dum = ptx.pyplmpsi(l,j,1,phi)
2035                Kcl += BOH['L='+str(l)][n-1][im]*pcrs*cosd(j*beta)       
2036        else:                #all but cubic
2037            pcrs,dum = ptx.pyplmpsi(l,n,1,phi)
2038            pcrs *= RSQPI
2039            if n == 0:
2040                pcrs /= SQ2
2041            if SGData['SGLaue'] in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']:
2042               if SGData['SGLaue'] in ['3mR','3m1','31m']: 
2043                   if n%6 == 3:
2044                       Kcl = pcrs*sind(n*beta)
2045                   else:
2046                       Kcl = pcrs*cosd(n*beta)
2047               else:
2048                   Kcl = pcrs*cosd(n*beta)
2049            else:
2050                Kcl = pcrs*(cosd(n*beta)+sind(n*beta))
2051        Fln[i] = SHCoef[term]*Kcl*Lnorm(l)
2052    ODFln = dict(zip(SHCoef.keys(),list(zip(SHCoef.values(),Fln))))
2053    return ODFln
2054   
2055def polfcal(ODFln,SamSym,psi,gam):
2056    '''Perform a pole figure computation.
2057    Note that the the number of gam values must either be 1 or must
2058    match psi. Updated for numpy 1.8.0
2059    '''
2060    import pytexture as ptx
2061    PolVal = np.ones_like(psi)
2062    for term in ODFln:
2063        if abs(ODFln[term][1]) > 1.e-3:
2064            l,m,n = eval(term.strip('C'))
2065            psrs,dum = ptx.pyplmpsi(l,m,len(psi),psi)
2066            if SamSym in ['-1','2/m']:
2067                if m:
2068                    Ksl = RSQPI*psrs*(cosd(m*gam)+sind(m*gam))
2069                else:
2070                    Ksl = RSQPI*psrs/SQ2
2071            else:
2072                if m:
2073                    Ksl = RSQPI*psrs*cosd(m*gam)
2074                else:
2075                    Ksl = RSQPI*psrs/SQ2
2076            PolVal += ODFln[term][1]*Ksl
2077    return PolVal
2078   
2079def invpolfcal(ODFln,SGData,phi,beta):
2080    'needs doc string'
2081    import pytexture as ptx
2082   
2083    invPolVal = np.ones_like(beta)
2084    for term in ODFln:
2085        if abs(ODFln[term][1]) > 1.e-3:
2086            l,m,n = eval(term.strip('C'))
2087            if SGData['SGLaue'] in ['m3','m3m']:
2088                Kcl = 0.0
2089                for j in range(0,l+1,4):
2090                    im = j//4
2091                    pcrs,dum = ptx.pyplmpsi(l,j,len(beta),phi)
2092                    Kcl += BOH['L=%d'%(l)][n-1][im]*pcrs*cosd(j*beta)       
2093            else:                #all but cubic
2094                pcrs,dum = ptx.pyplmpsi(l,n,len(beta),phi)
2095                pcrs *= RSQPI
2096                if n == 0:
2097                    pcrs /= SQ2
2098                if SGData['SGLaue'] in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']:
2099                   if SGData['SGLaue'] in ['3mR','3m1','31m']: 
2100                       if n%6 == 3:
2101                           Kcl = pcrs*sind(n*beta)
2102                       else:
2103                           Kcl = pcrs*cosd(n*beta)
2104                   else:
2105                       Kcl = pcrs*cosd(n*beta)
2106                else:
2107                    Kcl = pcrs*(cosd(n*beta)+sind(n*beta))
2108            invPolVal += ODFln[term][1]*Kcl
2109    return invPolVal
2110   
2111   
2112def textureIndex(SHCoef):
2113    'needs doc string'
2114    Tindx = 1.0
2115    for term in SHCoef:
2116        l = eval(term.strip('C'))[0]
2117        Tindx += SHCoef[term]**2/(2.0*l+1.)
2118    return Tindx
2119   
2120# self-test materials follow.
2121selftestlist = []
2122'''Defines a list of self-tests'''
2123selftestquiet = True
2124def _ReportTest():
2125    'Report name and doc string of current routine when ``selftestquiet`` is False'
2126    if not selftestquiet:
2127        import inspect
2128        caller = inspect.stack()[1][3]
2129        doc = eval(caller).__doc__
2130        if doc is not None:
2131            print('testing '+__file__+' with '+caller+' ('+doc+')')
2132        else:
2133            print('testing '+__file__()+" with "+caller)
2134NeedTestData = True
2135def TestData():
2136    array = np.array
2137    global NeedTestData
2138    NeedTestData = False
2139    global CellTestData
2140    # output from uctbx computed on platform darwin on 2010-05-28
2141    CellTestData = [
2142# cell, g, G, cell*, V, V*
2143  [(4, 4, 4, 90, 90, 90), 
2144   array([[  1.60000000e+01,   9.79717439e-16,   9.79717439e-16],
2145       [  9.79717439e-16,   1.60000000e+01,   9.79717439e-16],
2146       [  9.79717439e-16,   9.79717439e-16,   1.60000000e+01]]), array([[  6.25000000e-02,   3.82702125e-18,   3.82702125e-18],
2147       [  3.82702125e-18,   6.25000000e-02,   3.82702125e-18],
2148       [  3.82702125e-18,   3.82702125e-18,   6.25000000e-02]]), (0.25, 0.25, 0.25, 90.0, 90.0, 90.0), 64.0, 0.015625],
2149# cell, g, G, cell*, V, V*
2150  [(4.0999999999999996, 5.2000000000000002, 6.2999999999999998, 100, 80, 130), 
2151   array([[ 16.81      , -13.70423184,   4.48533243],
2152       [-13.70423184,  27.04      ,  -5.6887143 ],
2153       [  4.48533243,  -5.6887143 ,  39.69      ]]), array([[ 0.10206349,  0.05083339, -0.00424823],
2154       [ 0.05083339,  0.06344997,  0.00334956],
2155       [-0.00424823,  0.00334956,  0.02615544]]), (0.31947376387537696, 0.25189277536327803, 0.16172643497798223, 85.283666420376008, 94.716333579624006, 50.825714168082683), 100.98576357983838, 0.0099023858863968445],
2156# cell, g, G, cell*, V, V*
2157  [(3.5, 3.5, 6, 90, 90, 120), 
2158   array([[  1.22500000e+01,  -6.12500000e+00,   1.28587914e-15],
2159       [ -6.12500000e+00,   1.22500000e+01,   1.28587914e-15],
2160       [  1.28587914e-15,   1.28587914e-15,   3.60000000e+01]]), array([[  1.08843537e-01,   5.44217687e-02,   3.36690552e-18],
2161       [  5.44217687e-02,   1.08843537e-01,   3.36690552e-18],
2162       [  3.36690552e-18,   3.36690552e-18,   2.77777778e-02]]), (0.32991443953692895, 0.32991443953692895, 0.16666666666666669, 90.0, 90.0, 60.000000000000021), 63.652867178156257, 0.015710211406520427],
2163  ]
2164    global CoordTestData
2165    CoordTestData = [
2166# cell, ((frac, ortho),...)
2167  ((4,4,4,90,90,90,), [
2168 ((0.10000000000000001, 0.0, 0.0),(0.40000000000000002, 0.0, 0.0)),
2169 ((0.0, 0.10000000000000001, 0.0),(2.4492935982947065e-17, 0.40000000000000002, 0.0)),
2170 ((0.0, 0.0, 0.10000000000000001),(2.4492935982947065e-17, -2.4492935982947065e-17, 0.40000000000000002)),
2171 ((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(0.40000000000000013, 0.79999999999999993, 1.2)),
2172 ((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(0.80000000000000016, 1.2, 0.40000000000000002)),
2173 ((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(1.2, 0.80000000000000004, 0.40000000000000002)),
2174 ((0.5, 0.5, 0.5),(2.0, 1.9999999999999998, 2.0)),
2175]),
2176# cell, ((frac, ortho),...)
2177  ((4.1,5.2,6.3,100,80,130,), [
2178 ((0.10000000000000001, 0.0, 0.0),(0.40999999999999998, 0.0, 0.0)),
2179 ((0.0, 0.10000000000000001, 0.0),(-0.33424955703700043, 0.39834311042186865, 0.0)),
2180 ((0.0, 0.0, 0.10000000000000001),(0.10939835193016617, -0.051013289294572106, 0.6183281045774256)),
2181 ((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(0.069695941716497567, 0.64364635296002093, 1.8549843137322766)),
2182 ((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(-0.073350319180835066, 1.1440160419710339, 0.6183281045774256)),
2183 ((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(0.67089923785616512, 0.74567293154916525, 0.6183281045774256)),
2184 ((0.5, 0.5, 0.5),(0.92574397446582857, 1.7366491056364828, 3.0916405228871278)),
2185]),
2186# cell, ((frac, ortho),...)
2187  ((3.5,3.5,6,90,90,120,), [
2188 ((0.10000000000000001, 0.0, 0.0),(0.35000000000000003, 0.0, 0.0)),
2189 ((0.0, 0.10000000000000001, 0.0),(-0.17499999999999993, 0.3031088913245536, 0.0)),
2190 ((0.0, 0.0, 0.10000000000000001),(3.6739403974420595e-17, -3.6739403974420595e-17, 0.60000000000000009)),
2191 ((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(2.7675166561703527e-16, 0.60621778264910708, 1.7999999999999998)),
2192 ((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(0.17500000000000041, 0.90932667397366063, 0.60000000000000009)),
2193 ((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(0.70000000000000018, 0.6062177826491072, 0.60000000000000009)),
2194 ((0.5, 0.5, 0.5),(0.87500000000000067, 1.5155444566227676, 3.0)),
2195]),
2196]
2197    global LaueTestData             #generated by GSAS
2198    LaueTestData = {
2199    'R 3 m':[(4.,4.,6.,90.,90.,120.),((1,0,1,6),(1,0,-2,6),(0,0,3,2),(1,1,0,6),(2,0,-1,6),(2,0,2,6),
2200        (1,1,3,12),(1,0,4,6),(2,1,1,12),(2,1,-2,12),(3,0,0,6),(1,0,-5,6),(2,0,-4,6),(3,0,-3,6),(3,0,3,6),
2201        (0,0,6,2),(2,2,0,6),(2,1,4,12),(2,0,5,6),(3,1,-1,12),(3,1,2,12),(1,1,6,12),(2,2,3,12),(2,1,-5,12))],
2202    'R 3':[(4.,4.,6.,90.,90.,120.),((1,0,1,6),(1,0,-2,6),(0,0,3,2),(1,1,0,6),(2,0,-1,6),(2,0,2,6),(1,1,3,6),
2203        (1,1,-3,6),(1,0,4,6),(3,-1,1,6),(2,1,1,6),(3,-1,-2,6),(2,1,-2,6),(3,0,0,6),(1,0,-5,6),(2,0,-4,6),
2204        (2,2,0,6),(3,0,3,6),(3,0,-3,6),(0,0,6,2),(3,-1,4,6),(2,0,5,6),(2,1,4,6),(4,-1,-1,6),(3,1,-1,6),
2205        (3,1,2,6),(4,-1,2,6),(2,2,-3,6),(1,1,-6,6),(1,1,6,6),(2,2,3,6),(2,1,-5,6),(3,-1,-5,6))],
2206    'P 3':[(4.,4.,6.,90.,90.,120.),((0,0,1,2),(1,0,0,6),(1,0,1,6),(0,0,2,2),(1,0,-1,6),(1,0,2,6),(1,0,-2,6),
2207        (1,1,0,6),(0,0,3,2),(1,1,1,6),(1,1,-1,6),(1,0,3,6),(1,0,-3,6),(2,0,0,6),(2,0,-1,6),(1,1,-2,6),
2208        (1,1,2,6),(2,0,1,6),(2,0,-2,6),(2,0,2,6),(0,0,4,2),(1,1,-3,6),(1,1,3,6),(1,0,-4,6),(1,0,4,6),
2209        (2,0,-3,6),(2,1,0,6),(2,0,3,6),(3,-1,0,6),(2,1,1,6),(3,-1,-1,6),(2,1,-1,6),(3,-1,1,6),(1,1,4,6),
2210        (3,-1,2,6),(3,-1,-2,6),(1,1,-4,6),(0,0,5,2),(2,1,2,6),(2,1,-2,6),(3,0,0,6),(3,0,1,6),(2,0,4,6),
2211        (2,0,-4,6),(3,0,-1,6),(1,0,-5,6),(1,0,5,6),(3,-1,-3,6),(2,1,-3,6),(2,1,3,6),(3,-1,3,6),(3,0,-2,6),
2212        (3,0,2,6),(1,1,5,6),(1,1,-5,6),(2,2,0,6),(3,0,3,6),(3,0,-3,6),(0,0,6,2),(2,0,-5,6),(2,1,-4,6),
2213        (2,2,-1,6),(3,-1,-4,6),(2,2,1,6),(3,-1,4,6),(2,1,4,6),(2,0,5,6),(1,0,-6,6),(1,0,6,6),(4,-1,0,6),
2214        (3,1,0,6),(3,1,-1,6),(3,1,1,6),(4,-1,-1,6),(2,2,2,6),(4,-1,1,6),(2,2,-2,6),(3,1,2,6),(3,1,-2,6),
2215        (3,0,4,6),(3,0,-4,6),(4,-1,-2,6),(4,-1,2,6),(2,2,-3,6),(1,1,6,6),(1,1,-6,6),(2,2,3,6),(3,-1,5,6),
2216        (2,1,5,6),(2,1,-5,6),(3,-1,-5,6))],
2217    'P 3 m 1':[(4.,4.,6.,90.,90.,120.),((0,0,1,2),(1,0,0,6),(1,0,-1,6),(1,0,1,6),(0,0,2,2),(1,0,-2,6),
2218        (1,0,2,6),(1,1,0,6),(0,0,3,2),(1,1,1,12),(1,0,-3,6),(1,0,3,6),(2,0,0,6),(1,1,2,12),(2,0,1,6),
2219        (2,0,-1,6),(0,0,4,2),(2,0,-2,6),(2,0,2,6),(1,1,3,12),(1,0,-4,6),(1,0,4,6),(2,0,3,6),(2,1,0,12),
2220        (2,0,-3,6),(2,1,1,12),(2,1,-1,12),(1,1,4,12),(2,1,2,12),(0,0,5,2),(2,1,-2,12),(3,0,0,6),(1,0,-5,6),
2221        (3,0,1,6),(3,0,-1,6),(1,0,5,6),(2,0,4,6),(2,0,-4,6),(2,1,3,12),(2,1,-3,12),(3,0,-2,6),(3,0,2,6),
2222        (1,1,5,12),(3,0,-3,6),(0,0,6,2),(2,2,0,6),(3,0,3,6),(2,1,4,12),(2,2,1,12),(2,0,5,6),(2,1,-4,12),
2223        (2,0,-5,6),(1,0,-6,6),(1,0,6,6),(3,1,0,12),(3,1,-1,12),(3,1,1,12),(2,2,2,12),(3,1,2,12),
2224        (3,0,4,6),(3,1,-2,12),(3,0,-4,6),(1,1,6,12),(2,2,3,12))],
2225    'P 3 1 m':[(4.,4.,6.,90.,90.,120.),((0,0,1,2),(1,0,0,6),(0,0,2,2),(1,0,1,12),(1,0,2,12),(1,1,0,6),
2226        (0,0,3,2),(1,1,-1,6),(1,1,1,6),(1,0,3,12),(2,0,0,6),(2,0,1,12),(1,1,2,6),(1,1,-2,6),(2,0,2,12),
2227        (0,0,4,2),(1,1,-3,6),(1,1,3,6),(1,0,4,12),(2,1,0,12),(2,0,3,12),(2,1,1,12),(2,1,-1,12),(1,1,-4,6),
2228        (1,1,4,6),(0,0,5,2),(2,1,-2,12),(2,1,2,12),(3,0,0,6),(1,0,5,12),(2,0,4,12),(3,0,1,12),(2,1,-3,12),
2229        (2,1,3,12),(3,0,2,12),(1,1,5,6),(1,1,-5,6),(3,0,3,12),(0,0,6,2),(2,2,0,6),(2,1,-4,12),(2,0,5,12),
2230        (2,2,-1,6),(2,2,1,6),(2,1,4,12),(3,1,0,12),(1,0,6,12),(2,2,2,6),(3,1,-1,12),(2,2,-2,6),(3,1,1,12),
2231        (3,1,-2,12),(3,0,4,12),(3,1,2,12),(1,1,-6,6),(2,2,3,6),(2,2,-3,6),(1,1,6,6))],
2232    }
2233   
2234    global FLnhTestData
2235    FLnhTestData = [{
2236    'C(4,0,0)': (0.965, 0.42760447),
2237    'C(2,0,0)': (1.0122, -0.80233610),
2238    'C(2,0,2)': (0.0061, 8.37491546E-03),
2239    'C(6,0,4)': (-0.0898, 4.37985696E-02),
2240    'C(6,0,6)': (-0.1369, -9.04081762E-02),
2241    'C(6,0,0)': (0.5935, -0.18234928),
2242    'C(4,0,4)': (0.1872, 0.16358127),
2243    'C(6,0,2)': (0.6193, 0.27573633),
2244    'C(4,0,2)': (-0.1897, 0.12530720)},[1,0,0]]
2245def test0():
2246    if NeedTestData: TestData()
2247    msg = 'test cell2Gmat, fillgmat, Gmat2cell'
2248    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2249        G, g = cell2Gmat(cell)
2250        assert np.allclose(G,tG),msg
2251        assert np.allclose(g,tg),msg
2252        tcell = Gmat2cell(g)
2253        assert np.allclose(cell,tcell),msg
2254        tcell = Gmat2cell(G)
2255        assert np.allclose(tcell,trcell),msg
2256if __name__ == '__main__': selftestlist.append(test0)
2257
2258def test1():
2259    'test cell2A and A2Gmat'
2260    _ReportTest()
2261    if NeedTestData: TestData()
2262    msg = 'test cell2A and A2Gmat'
2263    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2264        G, g = A2Gmat(cell2A(cell))
2265        assert np.allclose(G,tG),msg
2266        assert np.allclose(g,tg),msg
2267if __name__ == '__main__': selftestlist.append(test1)
2268
2269def test2():
2270    'test Gmat2A, A2cell, A2Gmat, Gmat2cell'
2271    _ReportTest()
2272    if NeedTestData: TestData()
2273    msg = 'test Gmat2A, A2cell, A2Gmat, Gmat2cell'
2274    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2275        G, g = cell2Gmat(cell)
2276        tcell = A2cell(Gmat2A(G))
2277        assert np.allclose(cell,tcell),msg
2278if __name__ == '__main__': selftestlist.append(test2)
2279
2280def test3():
2281    'test invcell2Gmat'
2282    _ReportTest()
2283    if NeedTestData: TestData()
2284    msg = 'test invcell2Gmat'
2285    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2286        G, g = invcell2Gmat(trcell)
2287        assert np.allclose(G,tG),msg
2288        assert np.allclose(g,tg),msg
2289if __name__ == '__main__': selftestlist.append(test3)
2290
2291def test4():
2292    'test calc_rVsq, calc_rV, calc_V'
2293    _ReportTest()
2294    if NeedTestData: TestData()
2295    msg = 'test calc_rVsq, calc_rV, calc_V'
2296    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2297        assert np.allclose(calc_rV(cell2A(cell)),trV), msg
2298        assert np.allclose(calc_V(cell2A(cell)),tV), msg
2299if __name__ == '__main__': selftestlist.append(test4)
2300
2301def test5():
2302    'test A2invcell'
2303    _ReportTest()
2304    if NeedTestData: TestData()
2305    msg = 'test A2invcell'
2306    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2307        rcell = A2invcell(cell2A(cell))
2308        assert np.allclose(rcell,trcell),msg
2309if __name__ == '__main__': selftestlist.append(test5)
2310
2311def test6():
2312    'test cell2AB'
2313    _ReportTest()
2314    if NeedTestData: TestData()
2315    msg = 'test cell2AB'
2316    for (cell,coordlist) in CoordTestData:
2317        A,B = cell2AB(cell)
2318        for (frac,ortho) in coordlist:
2319            to = np.inner(A,frac)
2320            tf = np.inner(B,to)
2321            assert np.allclose(ortho,to), msg
2322            assert np.allclose(frac,tf), msg
2323            to = np.sum(A*frac,axis=1)
2324            tf = np.sum(B*to,axis=1)
2325            assert np.allclose(ortho,to), msg
2326            assert np.allclose(frac,tf), msg
2327if __name__ == '__main__': selftestlist.append(test6)
2328
2329def test7():
2330    'test GetBraviasNum(...) and GenHBravais(...)'
2331    _ReportTest()
2332    import os.path
2333    import sys
2334    import GSASIIspc as spc
2335    testdir = os.path.join(os.path.split(os.path.abspath( __file__ ))[0],'testinp')
2336    if os.path.exists(testdir):
2337        if testdir not in sys.path: sys.path.insert(0,testdir)
2338    import sgtbxlattinp
2339    derror = 1e-4
2340    def indexmatch(hklin, hkllist, system):
2341        for hklref in hkllist:
2342            hklref = list(hklref)
2343            # these permutations are far from complete, but are sufficient to
2344            # allow the test to complete
2345            if system == 'cubic':
2346                permlist = [(1,2,3),(1,3,2),(2,1,3),(2,3,1),(3,1,2),(3,2,1),]
2347            elif system == 'monoclinic':
2348                permlist = [(1,2,3),(-1,2,-3)]
2349            else:
2350                permlist = [(1,2,3)]
2351
2352            for perm in permlist:
2353                hkl = [abs(i) * hklin[abs(i)-1] / i for i in perm]
2354                if hkl == hklref: return True
2355                if [-i for i in hkl] == hklref: return True
2356        else:
2357            return False
2358
2359    for key in sgtbxlattinp.sgtbx7:
2360        spdict = spc.SpcGroup(key)
2361        cell = sgtbxlattinp.sgtbx7[key][0]
2362        system = spdict[1]['SGSys']
2363        center = spdict[1]['SGLatt']
2364
2365        bravcode = GetBraviasNum(center, system)
2366
2367        g2list = GenHBravais(sgtbxlattinp.dmin, bravcode, cell2A(cell))
2368
2369        assert len(sgtbxlattinp.sgtbx7[key][1]) == len(g2list), 'Reflection lists differ for %s' % key
2370        for h,k,l,d,num in g2list:
2371            for hkllist,dref in sgtbxlattinp.sgtbx7[key][1]: 
2372                if abs(d-dref) < derror:
2373                    if indexmatch((h,k,l,), hkllist, system):
2374                        break
2375            else:
2376                assert 0,'No match for %s at %s (%s)' % ((h,k,l),d,key)
2377if __name__ == '__main__': selftestlist.append(test7)
2378
2379def test8():
2380    'test GenHLaue'
2381    _ReportTest()
2382    import GSASIIspc as spc
2383    import sgtbxlattinp
2384    derror = 1e-4
2385    dmin = sgtbxlattinp.dmin
2386
2387    def indexmatch(hklin, hklref, system, axis):
2388        # these permutations are far from complete, but are sufficient to
2389        # allow the test to complete
2390        if system == 'cubic':
2391            permlist = [(1,2,3),(1,3,2),(2,1,3),(2,3,1),(3,1,2),(3,2,1),]
2392        elif system == 'monoclinic' and axis=='b':
2393            permlist = [(1,2,3),(-1,2,-3)]
2394        elif system == 'monoclinic' and axis=='a':
2395            permlist = [(1,2,3),(1,-2,-3)]
2396        elif system == 'monoclinic' and axis=='c':
2397            permlist = [(1,2,3),(-1,-2,3)]
2398        elif system == 'trigonal':
2399            permlist = [(1,2,3),(2,1,3),(-1,-2,3),(-2,-1,3)]
2400        elif system == 'rhombohedral':
2401            permlist = [(1,2,3),(2,3,1),(3,1,2)]
2402        else:
2403            permlist = [(1,2,3)]
2404
2405        hklref = list(hklref)
2406        for perm in permlist:
2407            hkl = [abs(i) * hklin[abs(i)-1] / i for i in perm]
2408            if hkl == hklref: return True
2409            if [-i for i in hkl] == hklref: return True
2410        return False
2411
2412    for key in sgtbxlattinp.sgtbx8:
2413        spdict = spc.SpcGroup(key)[1]
2414        cell = sgtbxlattinp.sgtbx8[key][0]
2415        Axis = spdict['SGUniq']
2416        system = spdict['SGSys']
2417
2418        g2list = GenHLaue(dmin,spdict,cell2A(cell))
2419        #if len(g2list) != len(sgtbxlattinp.sgtbx8[key][1]):
2420        #    print 'failed',key,':' ,len(g2list),'vs',len(sgtbxlattinp.sgtbx8[key][1])
2421        #    print 'GSAS-II:'
2422        #    for h,k,l,d in g2list: print '  ',(h,k,l),d
2423        #    print 'SGTBX:'
2424        #    for hkllist,dref in sgtbxlattinp.sgtbx8[key][1]: print '  ',hkllist,dref
2425        assert len(g2list) == len(sgtbxlattinp.sgtbx8[key][1]), (
2426            'Reflection lists differ for %s' % key
2427            )
2428        #match = True
2429        for h,k,l,d in g2list:
2430            for hkllist,dref in sgtbxlattinp.sgtbx8[key][1]: 
2431                if abs(d-dref) < derror:
2432                    if indexmatch((h,k,l,), hkllist, system, Axis): break
2433            else:
2434                assert 0,'No match for %s at %s (%s)' % ((h,k,l),d,key)
2435                #match = False
2436        #if not match:
2437            #for hkllist,dref in sgtbxlattinp.sgtbx8[key][1]: print '  ',hkllist,dref
2438            #print center, Laue, Axis, system
2439if __name__ == '__main__': selftestlist.append(test8)
2440           
2441def test9():
2442    'test GenHLaue'
2443    _ReportTest()
2444    import GSASIIspc as G2spc
2445    if NeedTestData: TestData()
2446    for spc in LaueTestData:
2447        data = LaueTestData[spc]
2448        cell = data[0]
2449        hklm = np.array(data[1])
2450        H = hklm[-1][:3]
2451        hklO = hklm.T[:3].T
2452        A = cell2A(cell)
2453        dmin = 1./np.sqrt(calc_rDsq(H,A))
2454        SGData = G2spc.SpcGroup(spc)[1]
2455        hkls = np.array(GenHLaue(dmin,SGData,A))
2456        hklN = hkls.T[:3].T
2457        #print spc,hklO.shape,hklN.shape
2458        err = True
2459        for H in hklO:
2460            if H not in hklN:
2461                print ('%d %s'%(H,' missing from hkl from GSASII'))
2462                err = False
2463        assert(err)
2464if __name__ == '__main__': selftestlist.append(test9)
2465       
2466       
2467   
2468
2469if __name__ == '__main__':
2470    # run self-tests
2471    selftestquiet = False
2472    for test in selftestlist:
2473        test()
2474    print ("OK")
Note: See TracBrowser for help on using the repository browser.