source: trunk/GSASIIlattice.py @ 4898

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

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