source: trunk/GSASIIlattice.py @ 4791

Last change on this file since 4791 was 4791, checked in by toby, 2 years ago

docs chapter numbering

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