source: trunk/GSASIIlattice.py @ 5116

Last change on this file since 5116 was 5116, checked in by vondreele, 10 months ago

fixes to EDX calcs for LeBail?; start on Pawley for EDX.

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