source: trunk/GSASIIlattice.py @ 3736

Last change on this file since 3736 was 3736, checked in by vondreele, 5 years ago

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