source: trunk/GSASIIlattice.py @ 4837

Last change on this file since 4837 was 4826, checked in by vondreele, 8 months ago

fix PDF GUI stuff transmission as percent, plot limits for g(r)
fix problem with 'PKS' unit cell indexing

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