source: trunk/GSASIIlattice.py @ 4434

Last change on this file since 4434 was 4434, checked in by vondreele, 17 months ago

Change Force=True for anisotropic atms in call to GenAtom? in FillUnitCell? - otherwise breaks RMCProfile data prep
fix a 'Va' atom check in RMCProfile setup

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