source: trunk/GSASIIlattice.py

Last change on this file was 5038, checked in by toby, 28 hours ago

Major revision to constraints including an update to the Sequential Refinement tutorial

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