source: trunk/GSASIIlattice.py @ 5105

Last change on this file since 5105 was 5105, checked in by vondreele, 7 months ago

fix monoclinic b form for CellDijCorr? line 275; missing last zero.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 121.7 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIlattice: Unit cells*
4---------------------------
5
6Perform lattice-related computations
7
8Note that *G* is the reciprocal lattice tensor, and *g* is its inverse,
9:math:`G = g^{-1}`, where
10
11  .. math::
12
13   g = \\left( \\begin{matrix}
14   a^2 & a b\\cos\\gamma & a c\\cos\\beta \\\\
15   a b\\cos\\gamma & b^2 & b c \\cos\\alpha \\\\
16   a c\\cos\\beta &  b c \\cos\\alpha & c^2
17   \\end{matrix}\\right)
18
19The "*A* tensor" terms are defined as
20:math:`A = (\\begin{matrix} G_{11} & G_{22} & G_{33} & 2G_{12} & 2G_{13} & 2G_{23}\\end{matrix})` and *A* can be used in this fashion:
21:math:`d^* = \\sqrt {A_0 h^2 + A_1 k^2 + A_2 l^2 + A_3 hk + A_4 hl + A_5 kl}`, where
22*d* is the d-spacing, and :math:`d^*` is the reciprocal lattice spacing,
23:math:`Q = 2 \\pi d^* = 2 \\pi / d`.
24Note that GSAS-II variables ``p::Ai`` (i = 0, 1,... 5) and ``p`` is a phase number are
25used for the *Ai* values. See :func:`A2cell`, :func:`cell2A` for interconversion between A and
26unit cell parameters; :func:`cell2Gmat` :func:`Gmat2cell` for G and cell parameters.
27
28When the hydrostatic/elastic strain coefficients (*Dij*, :math:`D_{ij}`) are used, they are added to the
29*A* tensor terms (Ai, :math:`A_{i}`) so that A is redefined
30:math:`A = (\\begin{matrix} A_{0} + D_{11} & A_{1} + D_{22} & A_{2} + D_{33} & A_{3} + D_{12} & A_{4} + D_{13} & A_{5} + D_{23}\\end{matrix})`. See :func:`cellDijFill`.
31Note that GSAS-II variables ``p:h:Dij`` (i,j = 1, 2, 3) and ``p`` is a phase number
32and ``h`` a histogram number are used for the *Dij* values.
33'''
34########### SVN repository information ###################
35# $Date: 2021-11-27 15:55:47 +0000 (Sat, 27 Nov 2021) $
36# $Author: vondreele $
37# $Revision: 5105 $
38# $URL: trunk/GSASIIlattice.py $
39# $Id: GSASIIlattice.py 5105 2021-11-27 15:55:47Z 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: 5105 $")
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    '''
1144    if 'T' in Inst['Type'][0]:
1145        return TOF2dsp(Inst,pos)
1146    else:   #'PKS', 'C' or 'B'
1147        wave = G2mth.getWave(Inst)
1148        return wave/(2.0*sind((pos-Inst.get('Zero',[0,0])[1])/2.0))
1149       
1150def TOF2dsp(Inst,Pos):
1151    ''' convert powder pattern TOF, musec to d-spacing by successive approximation
1152    Pos can be numpy array
1153    '''
1154    def func(d,pos,Inst):       
1155        return (pos-Inst['difA'][1]*d**2-Inst['Zero'][1]-Inst['difB'][1]/d)/Inst['difC'][1]
1156    dsp0 = Pos/Inst['difC'][1]
1157    N = 0
1158    while True:      #successive approximations
1159        dsp = func(dsp0,Pos,Inst)
1160        if np.allclose(dsp,dsp0,atol=0.000001):
1161            return dsp
1162        dsp0 = dsp
1163        N += 1
1164        if N > 10:
1165            return dsp
1166   
1167def Dsp2pos(Inst,dsp):
1168    ''' convert d-spacing to powder pattern position (2-theta or TOF, musec)
1169    '''
1170    if 'T' in Inst['Type'][0]:
1171        pos = Inst['difC'][1]*dsp+Inst['Zero'][1]+Inst['difA'][1]*dsp**2+Inst.get('difB',[0,0,False])[1]/dsp
1172    else:   #'C' or 'B'
1173        wave = G2mth.getWave(Inst)
1174        val = min(0.995,wave/(2.*dsp))  #set max at 168deg
1175        pos = 2.0*asind(val)+Inst.get('Zero',[0,0])[1]             
1176    return pos
1177   
1178def getPeakPos(dataType,parmdict,dsp):
1179    ''' convert d-spacing to powder pattern position (2-theta or TOF, musec)
1180    '''
1181    if 'T' in dataType:
1182        pos = parmdict['difC']*dsp+parmdict['difA']*dsp**2+parmdict['difB']/dsp+parmdict['Zero']
1183    else:   #'C' or 'B'
1184        pos = 2.0*asind(parmdict['Lam']/(2.*dsp))+parmdict['Zero']
1185    return pos
1186                   
1187def calc_rDsq(H,A):
1188    'needs doc string'
1189    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]
1190    return rdsq
1191   
1192def calc_rDsq2(H,G):
1193    'needs doc string'
1194    return np.inner(H,np.inner(G,H))
1195   
1196def calc_rDsqSS(H,A,vec):
1197    'needs doc string'
1198    rdsq = calc_rDsq(H[:3]+(H[3]*vec).T,A)
1199    return rdsq
1200       
1201def calc_rDsqZ(H,A,Z,tth,lam):
1202    'needs doc string'
1203    rdsq = calc_rDsq(H,A)+Z*sind(tth)*2.0*rpd/lam**2
1204    return rdsq
1205       
1206def calc_rDsqZSS(H,A,vec,Z,tth,lam):
1207    'needs doc string'
1208    rdsq = calc_rDsq(H[:3]+(H[3][:,np.newaxis]*vec).T,A)+Z*sind(tth)*2.0*rpd/lam**2
1209    return rdsq
1210       
1211def calc_rDsqT(H,A,Z,tof,difC):
1212    'needs doc string'
1213    rdsq = calc_rDsq(H,A)+Z/difC
1214    return rdsq
1215       
1216def calc_rDsqTSS(H,A,vec,Z,tof,difC):
1217    'needs doc string'
1218    rdsq = calc_rDsq(H[:3]+(H[3][:,np.newaxis]*vec).T,A)+Z/difC
1219    return rdsq
1220   
1221def PlaneIntercepts(Amat,H,phase,stack):
1222    ''' find unit cell intercepts for a stack of hkl planes
1223    '''
1224    Steps = range(-1,2,2)
1225    if stack:
1226        Steps = range(-10,10,1)
1227    Stack = []
1228    Ux = np.array([[0,0],[1,0],[1,1],[0,1]])
1229    for step in Steps:
1230        HX = []
1231        for i in [0,1,2]:
1232            if H[i]:
1233               h,k,l = [(i+1)%3,(i+2)%3,(i+3)%3]
1234               for j in [0,1,2,3]:
1235                    hx = [0,0,0]
1236                    intcpt = ((phase)/360.+step-H[h]*Ux[j,0]-H[k]*Ux[j,1])/H[l]
1237                    if 0. <= intcpt <= 1.:                       
1238                        hx[h] = Ux[j,0]
1239                        hx[k] = Ux[j,1]
1240                        hx[l] = intcpt
1241                        HX.append(hx)
1242        if len(HX)> 2:
1243            HX = np.array(HX)
1244            DX = np.inner(HX-HX[0],Amat)
1245            D = np.sqrt(np.sum(DX**2,axis=1))
1246            Dsort = np.argsort(D)
1247            HX = HX[Dsort]
1248            DX = DX[Dsort]
1249            D = D[Dsort]
1250            DX[1:,:] = DX[1:,:]/D[1:,nxs]
1251            A = 2.*np.ones(HX.shape[0])
1252            A[1:] = [np.dot(DX[1],dx) for dx in DX[1:]]
1253            HX = HX[np.argsort(A)]
1254            Stack.append(HX)
1255    return Stack
1256       
1257def MaxIndex(dmin,A):
1258    'needs doc string'
1259    Hmax = [0,0,0]
1260    try:
1261        cell = A2cell(A)
1262    except:
1263        cell = [1.,1.,1.,90.,90.,90.]
1264    for i in range(3):
1265        Hmax[i] = int(np.round(cell[i]/dmin))
1266    return Hmax
1267   
1268def transposeHKLF(transMat,Super,refList):
1269    ''' Apply transformation matrix to hkl(m)
1270    param: transmat: 3x3 or 4x4 array
1271    param: Super: 0 or 1 for extra index
1272    param: refList list of h,k,l,....
1273    return: newRefs transformed list of h',k',l',,,
1274    return: badRefs list of noninteger h',k',l'...
1275    '''
1276    newRefs = np.copy(refList)
1277    badRefs = []
1278    for H in newRefs:
1279        newH = np.inner(transMat,H[:3+Super])
1280        H[:3+Super] = np.rint(newH)
1281        if not np.allclose(newH,H[:3+Super],atol=0.01):
1282            badRefs.append(newH)
1283    return newRefs,badRefs
1284   
1285def sortHKLd(HKLd,ifreverse,ifdup,ifSS=False):
1286    '''sort reflection list on d-spacing; can sort in either order
1287
1288    :param HKLd: a list of [h,k,l,d,...];
1289    :param ifreverse: True for largest d first
1290    :param ifdup: True if duplicate d-spacings allowed
1291    :return: sorted reflection list
1292    '''
1293    T = []
1294    N = 3
1295    if ifSS:
1296        N = 4
1297    for i,H in enumerate(HKLd):
1298        if ifdup:
1299            T.append((H[N],i))
1300        else:
1301            T.append(H[N])           
1302    D = dict(zip(T,HKLd))
1303    T.sort()
1304    if ifreverse:
1305        T.reverse()
1306    X = []
1307    okey = ''
1308    for key in T: 
1309        if key != okey: X.append(D[key])    #remove duplicate d-spacings
1310        okey = key
1311    return X
1312   
1313def SwapIndx(Axis,H):
1314    'needs doc string'
1315    if Axis in [1,-1]:
1316        return H
1317    elif Axis in [2,-3]:
1318        return [H[1],H[2],H[0]]
1319    else:
1320        return [H[2],H[0],H[1]]
1321   
1322def SwapItems(Alist,pos1,pos2):
1323    'exchange 2 items in a list'
1324    try:
1325        get = Alist[pos1],Alist[pos2]
1326        Alist[pos2],Alist[pos1] = get
1327    except IndexError:
1328        pass
1329    return Alist
1330       
1331def Rh2Hx(Rh):
1332    'needs doc string'
1333    Hx = [0,0,0]
1334    Hx[0] = Rh[0]-Rh[1]
1335    Hx[1] = Rh[1]-Rh[2]
1336    Hx[2] = np.sum(Rh)
1337    return Hx
1338   
1339def Hx2Rh(Hx):
1340    'needs doc string'
1341    Rh = [0,0,0]
1342    itk = -Hx[0]+Hx[1]+Hx[2]
1343    if itk%3 != 0:
1344        return 0        #error - not rhombohedral reflection
1345    else:
1346        Rh[1] = itk//3
1347        Rh[0] = Rh[1]+Hx[0]
1348        Rh[2] = Rh[1]-Hx[1]
1349        if Rh[0] < 0:
1350            for i in range(3):
1351                Rh[i] = -Rh[i]
1352        return Rh
1353       
1354def CentCheck(Cent,H):
1355    'needs doc string'
1356    h,k,l = H
1357    if Cent == 'A' and (k+l)%2:
1358        return False
1359    elif Cent == 'B' and (h+l)%2:
1360        return False
1361    elif Cent == 'C' and (h+k)%2:
1362        return False
1363    elif Cent == 'I' and (h+k+l)%2:
1364        return False
1365    elif Cent == 'F' and ((h+k)%2 or (h+l)%2 or (k+l)%2):
1366        return False
1367    elif Cent == 'R' and (-h+k+l)%3:
1368        return False
1369    else:
1370        return True
1371   
1372def RBsymCheck(Atoms,ct,cx,cs,AtLookUp,Amat,RBObjIds,SGData):
1373    """ Checks members of a rigid body to see if one is a symmetry equivalent of another.
1374    If so the atom site frac is set to zero.
1375
1376    :param Atoms: atom array as defined in GSAS-II; modified here
1377    :param ct: int location of atom type in Atoms item
1378    :param cx: int location of x,y,z,frac in Atoms item
1379    :param dict AtLookUp: atom lookup by Id table
1380    :param np.array Amat: crystal-to-Cartesian transformation matrix
1381    :param list RBObjIds: atom Id belonging to rigid body being tested
1382    :param dict SGData: GSAS-II space group info.
1383    :returns: Atoms with modified atom frac entries
1384   
1385    """
1386    for i,Id in enumerate(RBObjIds):
1387        XYZo = np.array(Atoms[AtLookUp[Id]][cx:cx+3])%1.
1388        typo = Atoms[AtLookUp[Id]][ct]
1389        for Jd in RBObjIds[i+1:]:
1390            if Atoms[AtLookUp[Jd]][ct] == typo:
1391                XYZt = Atoms[AtLookUp[Jd]][cx:cx+3]
1392                Xeqv = list(G2spc.GenAtom(np.array(XYZt)%1.,SGData,True))
1393                close = [np.allclose(np.inner(Amat,XYZo),np.inner(Amat,eqv[0]),atol=0.01) for eqv in Xeqv]
1394                if True in close:
1395                    Atoms[AtLookUp[Jd]][cx+3] = 0.0
1396        Sytsym,Mult = G2spc.SytSym(Atoms[AtLookUp[Id]][cx:cx+3],SGData)[:2]
1397        Atoms[AtLookUp[Id]][cs] = Sytsym
1398        Atoms[AtLookUp[Id]][cs+1] = Mult           
1399    return Atoms
1400                                   
1401def GetBraviasNum(center,system):
1402    """Determine the Bravais lattice number, as used in GenHBravais
1403   
1404    :param center: one of: 'P', 'C', 'I', 'F', 'R' (see SGLatt from GSASIIspc.SpcGroup)
1405    :param system: one of 'cubic', 'hexagonal', 'tetragonal', 'orthorhombic', 'trigonal' (for R)
1406      'monoclinic', 'triclinic' (see SGSys from GSASIIspc.SpcGroup)
1407    :return: a number between 0 and 13
1408      or throws a ValueError exception if the combination of center, system is not found (i.e. non-standard)
1409
1410    """
1411    if center.upper() == 'F' and system.lower() == 'cubic':
1412        return 0
1413    elif center.upper() == 'I' and system.lower() == 'cubic':
1414        return 1
1415    elif center.upper() == 'P' and system.lower() == 'cubic':
1416        return 2
1417    elif center.upper() == 'R' and system.lower() == 'trigonal':
1418        return 3
1419    elif center.upper() == 'P' and system.lower() == 'hexagonal':
1420        return 4
1421    elif center.upper() == 'I' and system.lower() == 'tetragonal':
1422        return 5
1423    elif center.upper() == 'P' and system.lower() == 'tetragonal':
1424        return 6
1425    elif center.upper() == 'F' and system.lower() == 'orthorhombic':
1426        return 7
1427    elif center.upper() == 'I' and system.lower() == 'orthorhombic':
1428        return 8
1429    elif center.upper() == 'A' and system.lower() == 'orthorhombic':
1430        return 9
1431    elif center.upper() == 'B' and system.lower() == 'orthorhombic':
1432        return 10
1433    elif center.upper() == 'C' and system.lower() == 'orthorhombic':
1434        return 11
1435    elif center.upper() == 'P' and system.lower() == 'orthorhombic':
1436        return 12
1437    elif center.upper() == 'C' and system.lower() == 'monoclinic':
1438        return 13
1439    elif center.upper() == 'P' and system.lower() == 'monoclinic':
1440        return 14
1441    elif center.upper() == 'P' and system.lower() == 'triclinic':
1442        return 15
1443    raise ValueError('non-standard Bravais lattice center=%s, cell=%s' % (center,system))
1444
1445def _GenHBravais_cctbx(dmin, Bravais, A, sg_type, uctbx_unit_cell, miller_index_generator):
1446    '''Alternate form of :func:`GenHBravais` that uses CCTBX internals
1447    '''
1448    g_inv = np.array([[A[0],   A[3]/2, A[4]/2],
1449                      [A[3]/2, A[1],   A[5]/2], 
1450                      [A[4]/2, A[5]/2, A[2]]])
1451    g = np.linalg.inv(g_inv)
1452    g_elems = (g[0][0], g[1][1], g[2][2], g[0][1], g[0][2], g[1][2])
1453    try:
1454        uc = uctbx_unit_cell(metrical_matrix=g_elems)
1455    except ValueError: # this function sometimes receives an A matrix that gives
1456                           # numbers <0 in the diagonal elems of g. Not sure why.
1457        return []
1458    #if sg_type is None:
1459    #    sg_type = make_sgtype(Bravais)
1460    mig = miller_index_generator(uc, sg_type, 0, dmin)
1461    result = []
1462    for h,k,l in mig:
1463        d = uc.d((h,k,l))
1464        result.append([h, k, l, d, -1])
1465    result.sort(key=lambda l: l[3], reverse=True)
1466    return result
1467
1468def GenHBravais(dmin, Bravais, A, cctbx_args=None):
1469    """Generate the positionally unique powder diffraction reflections
1470     
1471    :param dmin: minimum d-spacing in A
1472    :param Bravais: lattice type (see GetBraviasNum). Bravais is one of:
1473   
1474            * 0 F cubic
1475            * 1 I cubic
1476            * 2 P cubic
1477            * 3 R hexagonal (trigonal not rhombohedral)
1478            * 4 P hexagonal
1479            * 5 I tetragonal
1480            * 6 P tetragonal
1481            * 7 F orthorhombic
1482            * 8 I orthorhombic
1483            * 9 A orthorhombic
1484            * 10 B orthorhombic
1485            * 11 C orthorhombic
1486            * 12 P orthorhombic
1487            * 13 I monoclinic
1488            * 14 A monoclinic
1489            * 15 C monoclinic
1490            * 16 P monoclinic
1491            * 17 P triclinic
1492           
1493    :param A: reciprocal metric tensor elements as [G11,G22,G33,2*G12,2*G13,2*G23]
1494    :param dict cctbx_args: items defined in CCTBX:
1495
1496         * 'sg_type': value from cctbx.sgtbx.space_group_type(symmorphic_sgs[ibrav])
1497         * 'uctbx_unit_cell': pointer to :meth:`cctbx.uctbx.unit_cell`
1498         * 'miller_index_generator':  pointer to :meth:`cctbx.miller.index_generator`
1499
1500    :returns: HKL unique d list of [h,k,l,d,-1] sorted with largest d first
1501           
1502    """
1503    if cctbx_args:
1504        return _GenHBravais_cctbx(dmin, Bravais, A,
1505                    cctbx_args['sg_type'], cctbx_args['uctbx_unit_cell'], cctbx_args['miller_index_generator'])
1506   
1507    if Bravais in [9,14]:
1508        Cent = 'A'
1509    elif Bravais in [10,]:
1510        Cent = 'B'
1511    elif Bravais in [11,15]:
1512        Cent = 'C'
1513    elif Bravais in [1,5,8,13]:
1514        Cent = 'I'
1515    elif Bravais in [0,7]:
1516        Cent = 'F'
1517    elif Bravais in [3]:
1518        Cent = 'R'
1519    else:
1520        Cent = 'P'
1521    Hmax = MaxIndex(dmin,A)
1522    dminsq = 1./(dmin**2)
1523    HKL = []
1524    if Bravais == 17:                       #triclinic
1525        for l in range(-Hmax[2],Hmax[2]+1):
1526            for k in range(-Hmax[1],Hmax[1]+1):
1527                hmin = 0
1528                if (k < 0): hmin = 1
1529                if (k ==0 and l < 0): hmin = 1
1530                for h in range(hmin,Hmax[0]+1):
1531                    H=[h,k,l]
1532                    rdsq = calc_rDsq(H,A)
1533                    if 0 < rdsq <= dminsq:
1534                        HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
1535    elif Bravais in [13,14,15,16]:                #monoclinic - b unique
1536        Hmax = SwapIndx(2,Hmax)
1537        for h in range(Hmax[0]+1):
1538            for k in range(-Hmax[1],Hmax[1]+1):
1539                lmin = 0
1540                if k < 0:lmin = 1
1541                for l in range(lmin,Hmax[2]+1):
1542                    [h,k,l] = SwapIndx(-2,[h,k,l])
1543                    H = []
1544                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
1545                    if H:
1546                        rdsq = calc_rDsq(H,A)
1547                        if 0 < rdsq <= dminsq:
1548                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
1549                    [h,k,l] = SwapIndx(2,[h,k,l])
1550    elif Bravais in [7,8,9,10,11,12]:            #orthorhombic
1551        for h in range(Hmax[0]+1):
1552            for k in range(Hmax[1]+1):
1553                for l in range(Hmax[2]+1):
1554                    H = []
1555                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
1556                    if H:
1557                        rdsq = calc_rDsq(H,A)
1558                        if 0 < rdsq <= dminsq:
1559                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
1560    elif Bravais in [5,6]:                  #tetragonal
1561        for l in range(Hmax[2]+1):
1562            for k in range(Hmax[1]+1):
1563                for h in range(k,Hmax[0]+1):
1564                    H = []
1565                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
1566                    if H:
1567                        rdsq = calc_rDsq(H,A)
1568                        if 0 < rdsq <= dminsq:
1569                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
1570    elif Bravais in [3,4]:
1571        lmin = 0
1572        if Bravais == 3: lmin = -Hmax[2]                  #hexagonal/trigonal
1573        for l in range(lmin,Hmax[2]+1):
1574            for k in range(Hmax[1]+1):
1575                hmin = k
1576                if l < 0: hmin += 1
1577                for h in range(hmin,Hmax[0]+1):
1578                    H = []
1579                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
1580                    if H:
1581                        rdsq = calc_rDsq(H,A)
1582                        if 0 < rdsq <= dminsq:
1583                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
1584
1585    else:                                   #cubic
1586        for l in range(Hmax[2]+1):
1587            for k in range(l,Hmax[1]+1):
1588                for h in range(k,Hmax[0]+1):
1589                    H = []
1590                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
1591                    if H:
1592                        rdsq = calc_rDsq(H,A)
1593                        if 0 < rdsq <= dminsq:
1594                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
1595    return sortHKLd(HKL,True,False)
1596   
1597def getHKLmax(dmin,SGData,A):
1598    'finds maximum allowed hkl for given A within dmin'
1599    SGLaue = SGData['SGLaue']
1600    if SGLaue in ['3R','3mR']:        #Rhombohedral axes
1601        Hmax = [0,0,0]
1602        cell = A2cell(A)
1603        aHx = cell[0]*math.sqrt(2.0*(1.0-cosd(cell[3])))
1604        cHx = cell[0]*math.sqrt(3.0*(1.0+2.0*cosd(cell[3])))
1605        Hmax[0] = Hmax[1] = int(round(aHx/dmin))
1606        Hmax[2] = int(round(cHx/dmin))
1607        #print Hmax,aHx,cHx
1608    else:                           # all others
1609        Hmax = MaxIndex(dmin,A)
1610    return Hmax
1611   
1612def GenHLaue(dmin,SGData,A):
1613    """Generate the crystallographically unique powder diffraction reflections
1614    for a lattice and Bravais type
1615   
1616    :param dmin: minimum d-spacing
1617    :param SGData: space group dictionary with at least
1618   
1619        * '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'
1620        * 'SGLatt': lattice centering: one of 'P','A','B','C','I','F'
1621        * 'SGUniq': code for unique monoclinic axis one of 'a','b','c' (only if 'SGLaue' is '2/m') otherwise an empty string
1622       
1623    :param A: reciprocal metric tensor elements as [G11,G22,G33,2*G12,2*G13,2*G23]
1624    :return: HKL = list of [h,k,l,d] sorted with largest d first and is unique
1625            part of reciprocal space ignoring anomalous dispersion
1626           
1627    """
1628    import math
1629    SGLaue = SGData['SGLaue']
1630    SGLatt = SGData['SGLatt']
1631    SGUniq = SGData['SGUniq']
1632    #finds maximum allowed hkl for given A within dmin
1633    Hmax = getHKLmax(dmin,SGData,A)
1634       
1635    dminsq = 1./(dmin**2)
1636    HKL = []
1637    if SGLaue == '-1':                       #triclinic
1638        for l in range(-Hmax[2],Hmax[2]+1):
1639            for k in range(-Hmax[1],Hmax[1]+1):
1640                hmin = 0
1641                if (k < 0) or (k ==0 and l < 0): hmin = 1
1642                for h in range(hmin,Hmax[0]+1):
1643                    H = []
1644                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
1645                    if H:
1646                        rdsq = calc_rDsq(H,A)
1647                        if 0 < rdsq <= dminsq:
1648                            HKL.append([h,k,l,1./math.sqrt(rdsq)])
1649    elif SGLaue == '2/m':                #monoclinic
1650        axisnum = 1 + ['a','b','c'].index(SGUniq)
1651        Hmax = SwapIndx(axisnum,Hmax)
1652        for h in range(Hmax[0]+1):
1653            for k in range(-Hmax[1],Hmax[1]+1):
1654                lmin = 0
1655                if k < 0:lmin = 1
1656                for l in range(lmin,Hmax[2]+1):
1657                    [h,k,l] = SwapIndx(-axisnum,[h,k,l])
1658                    H = []
1659                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
1660                    if H:
1661                        rdsq = calc_rDsq(H,A)
1662                        if 0 < rdsq <= dminsq:
1663                            HKL.append([h,k,l,1./math.sqrt(rdsq)])
1664                    [h,k,l] = SwapIndx(axisnum,[h,k,l])
1665    elif SGLaue in ['mmm','4/m','6/m']:            #orthorhombic
1666        for l in range(Hmax[2]+1):
1667            for h in range(Hmax[0]+1):
1668                kmin = 1
1669                if SGLaue == 'mmm' or h ==0: kmin = 0
1670                for k in range(kmin,Hmax[1]+1):
1671                    H = []
1672                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
1673                    if H:
1674                        rdsq = calc_rDsq(H,A)
1675                        if 0 < rdsq <= dminsq:
1676                            HKL.append([h,k,l,1./math.sqrt(rdsq)])
1677    elif SGLaue in ['4/mmm','6/mmm']:                  #tetragonal & hexagonal
1678        for l in range(Hmax[2]+1):
1679            for h in range(Hmax[0]+1):
1680                for k in range(h+1):
1681                    H = []
1682                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
1683                    if H:
1684                        rdsq = calc_rDsq(H,A)
1685                        if 0 < rdsq <= dminsq:
1686                            HKL.append([h,k,l,1./math.sqrt(rdsq)])
1687    elif SGLaue in ['3m1','31m','3','3R','3mR']:                  #trigonals
1688        for l in range(-Hmax[2],Hmax[2]+1):
1689            hmin = 0
1690            if l < 0: hmin = 1
1691            for h in range(hmin,Hmax[0]+1):
1692                if SGLaue in ['3R','3']:
1693                    kmax = h
1694                    kmin = -int((h-1.)/2.)
1695                else:
1696                    kmin = 0
1697                    kmax = h
1698                    if SGLaue in ['3m1','3mR'] and l < 0: kmax = h-1
1699                    if SGLaue == '31m' and l < 0: kmin = 1
1700                for k in range(kmin,kmax+1):
1701                    H = []
1702                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
1703                    if SGLaue in ['3R','3mR']:
1704                        H = Hx2Rh(H)
1705                    if H:
1706                        rdsq = calc_rDsq(H,A)
1707                        if 0 < rdsq <= dminsq:
1708                            HKL.append([H[0],H[1],H[2],1./math.sqrt(rdsq)])
1709    else:                                   #cubic
1710        for h in range(Hmax[0]+1):
1711            for k in range(h+1):
1712                lmin = 0
1713                lmax = k
1714                if SGLaue =='m3':
1715                    lmax = h-1
1716                    if h == k: lmax += 1
1717                for l in range(lmin,lmax+1):
1718                    H = []
1719                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
1720                    if H:
1721                        rdsq = calc_rDsq(H,A)
1722                        if 0 < rdsq <= dminsq:
1723                            HKL.append([h,k,l,1./math.sqrt(rdsq)])
1724    return sortHKLd(HKL,True,True)
1725   
1726def GenPfHKLs(nMax,SGData,A):   
1727    """Generate the unique pole figure reflections for a lattice and Bravais type.
1728    Min d-spacing=1.0A & no more than nMax returned
1729   
1730    :param nMax: maximum number of hkls returned
1731    :param SGData: space group dictionary with at least
1732   
1733        * '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'
1734        * 'SGLatt': lattice centering: one of 'P','A','B','C','I','F'
1735        * 'SGUniq': code for unique monoclinic axis one of 'a','b','c' (only if 'SGLaue' is '2/m') otherwise an empty string
1736       
1737    :param A: reciprocal metric tensor elements as [G11,G22,G33,2*G12,2*G13,2*G23]
1738    :return: HKL = list of 'h k l' strings sorted with largest d first; no duplicate zones
1739           
1740    """
1741    HKL = np.array(GenHLaue(1.0,SGData,A)).T[:3].T     #strip d-spacings
1742    N = min(nMax,len(HKL))
1743    return ['%d %d %d'%(h[0],h[1],h[2]) for h in HKL[:N]]       
1744
1745def GenSSHLaue(dmin,SGData,SSGData,Vec,maxH,A):
1746    'needs a doc string'
1747    ifMag = False
1748    if 'MagSpGrp' in SGData:
1749        ifMag = True
1750    HKLs = []
1751    vec = np.array(Vec)
1752    vstar = np.sqrt(calc_rDsq(vec,A))     #find extra needed for -n SS reflections
1753    dvec = 1./(maxH*vstar+1./dmin)
1754    HKL = GenHLaue(dvec,SGData,A)       
1755    SSdH = [vec*h for h in range(-maxH,maxH+1)]
1756    SSdH = dict(zip(range(-maxH,maxH+1),SSdH))
1757    for h,k,l,d in HKL:
1758        ext = G2spc.GenHKLf([h,k,l],SGData)[0]  #h,k,l must be integral values here
1759        if not ext and d >= dmin:
1760            HKLs.append([h,k,l,0,d])
1761        for dH in SSdH:
1762            if dH:
1763                DH = SSdH[dH]
1764                H = [h+DH[0],k+DH[1],l+DH[2]]
1765                d = 1./np.sqrt(calc_rDsq(H,A))
1766                if d >= dmin:
1767                    HKLM = np.array([h,k,l,dH])
1768                    if (G2spc.checkSSLaue([h,k,l,dH],SGData,SSGData) and G2spc.checkSSextc(HKLM,SSGData)) or ifMag:
1769                        HKLs.append([h,k,l,dH,d])   
1770    return HKLs
1771   
1772def LaueUnique2(SGData,refList):
1773    ''' Impose Laue symmetry on hkl
1774   
1775    :param SGData: space group data from 'P '+Laue
1776    :param HKLF: np.array([[h,k,l,...]]) reflection set to be converted
1777   
1778    :return: HKLF new reflection array with imposed Laue symmetry
1779    '''
1780    for ref in refList:
1781        H = ref[:3]
1782        Uniq = G2spc.GenHKLf(H,SGData)[2]
1783        Uniq = G2mth.sortArray(G2mth.sortArray(G2mth.sortArray(Uniq,2),1),0)
1784        ref[:3] = Uniq[-1]
1785    return refList
1786   
1787def LaueUnique(Laue,HKLF):
1788    ''' Impose Laue symmetry on hkl
1789   
1790    :param str Laue: Laue symbol, as below
1791   
1792      centrosymmetric Laue groups::
1793       
1794            ['-1','2/m','112/m','2/m11','mmm','-42m','-4m2','4/mmm','-3','-3m',
1795            '-31m','-3m1','6/m','6/mmm','m3','m3m']
1796     
1797      noncentrosymmetric Laue groups::
1798     
1799           ['1','2','211','112','m','m11','11m','222','mm2','m2m','2mm',
1800           '4','-4','422','4mm','3','312','321','3m','31m','3m1','6','-6',
1801           '622','6mm','-62m','-6m2','23','432','-43m']
1802     
1803    :param HKLF: np.array([[h,k,l,...]]) reflection set to be converted
1804   
1805    :returns: HKLF new reflection array with imposed Laue symmetry
1806    '''
1807   
1808    HKLFT = HKLF.T
1809    mat41 = np.array([[0,1,0],[-1,0,0],[0,0,1]])    #hkl -> k,-h,l
1810    mat43 = np.array([[0,-1,0],[1,0,0],[0,0,1]])    #hkl -> -k,h,l
1811    mat4bar = np.array([[0,-1,0],[1,0,0],[0,0,-1]]) #hkl -> k,-h,-l
1812    mat31 = np.array([[-1,-1,0],[1,0,0],[0,0,1]])   #hkl -> ihl = -h-k,h,l
1813    mat32 = np.array([[0,1,0],[-1,-1,0],[0,0,1]])   #hkl -> kil = k,-h-k,l
1814    matd3 = np.array([[0,1,0],[0,0,1],[1,0,0]])     #hkl -> k,l,h
1815    matd3q = np.array([[0,0,-1],[-1,0,0],[0,1,0]])  #hkl -> -l,-h,k
1816    matd3t = np.array([[0,0,-1],[1,0,0],[0,-1,0]])  #hkl -> -l,h,-k
1817    mat6 = np.array([[1,1,0],[-1,0,0],[0,0,1]])     #hkl -> h+k,-h,l really 65
1818    matdm = np.array([[0,1,0],[1,0,0],[0,0,1]])     #hkl -> k,h,l
1819    matdmp = np.array([[-1,-1,0],[0,1,0],[0,0,1]])  #hkl -> -h-k,k,l
1820    matkm = np.array([[-1,0,0],[1,1,0],[0,0,1]])    #hkl -> -h,h+k,l
1821    matd2 = np.array([[0,1,0],[1,0,0],[0,0,-1]])    #hkl -> k,h,-l
1822    matdm3 = np.array([[1,0,0],[0,0,1],[0,1,0]])    #hkl -> h,l,k
1823    mat2d43 = np.array([[0,1,0],[1,0,0],[0,0,1]])   #hkl -> k,-h,l
1824    matk2 = np.array([[-1,0,0],[1,1,0],[0,0,-1]])   #hkl -> -h,-i,-l
1825    #triclinic
1826    if Laue == '1': #ok
1827        pass
1828    elif Laue == '-1':  #ok
1829        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
1830        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]<0),HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
1831        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
1832    #monoclinic
1833    #noncentrosymmetric - all ok
1834    elif Laue == '2': 
1835        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3])
1836        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3])
1837    elif Laue == '1 1 2':
1838        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1839        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]<0),HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1840    elif Laue == '2 1 1':   
1841        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1842        HKLFT[:3] = np.where((HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1843    elif Laue == 'm':
1844        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1845    elif Laue == 'm 1 1':
1846        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1847    elif Laue == '1 1 m':
1848        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1849    #centrosymmetric - all ok
1850    elif Laue == '2/m 1 1':       
1851        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1852        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1853        HKLFT[:3] = np.where((HKLFT[2]*HKLFT[0]==0)&(HKLFT[1]<0),HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1854    elif Laue == '2/m':
1855        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3])
1856        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1857        HKLFT[:3] = np.where((HKLFT[0]*HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1858    elif Laue == '1 1 2/m':
1859        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1860        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1861        HKLFT[:3] = np.where((HKLFT[1]*HKLFT[2]==0)&(HKLFT[0]<0),HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1862    #orthorhombic
1863    #noncentrosymmetric - all OK
1864    elif Laue == '2 2 2':
1865        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1866        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1867        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1868        HKLFT[:3] = np.where((HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1869    elif Laue == 'm m 2':
1870        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1871        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1872    elif Laue == '2 m m': 
1873        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1874        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1875    elif Laue == 'm 2 m':
1876        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1877        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1878    #centrosymmetric - all ok
1879    elif Laue == 'm m m':
1880        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1881        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1882        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1883    #tetragonal
1884    #noncentrosymmetric - all ok
1885    elif Laue == '4':
1886        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1887        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1888        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]>0),np.squeeze(np.inner(HKLF[:,:3],mat41[nxs,:,:])).T,HKLFT[:3])
1889    elif Laue == '-4': 
1890        HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])     
1891        HKLFT[:3] = np.where(HKLFT[0]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1892        HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])     
1893        HKLFT[:3] = np.where(HKLFT[1]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1894        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1895    elif Laue == '4 2 2':
1896        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1897        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1898        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1899        HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]<HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1900        HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])   #in lieu od 2-fold
1901    elif Laue == '4 m m':
1902        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1903        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1904        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1905        HKLFT[:3] = np.where(HKLFT[0]<HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1906    elif Laue == '-4 2 m':
1907        HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])     
1908        HKLFT[:3] = np.where(HKLFT[0]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1909        HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])     
1910        HKLFT[:3] = np.where(HKLFT[1]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1911        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1912        HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1913        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1914    elif Laue == '-4 m 2':
1915        HKLFT[:3] = np.where(HKLFT[2]<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[2]==0)&(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[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])     
1919        HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]==0),np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1920        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3]) 
1921        HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[0]>HKLFT[1]),np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1922    #centrosymmetric - all ok
1923    elif Laue == '4/m':
1924        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1925        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1926        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1927        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]>0),np.squeeze(np.inner(HKLF[:,:3],mat41[nxs,:,:])).T,HKLFT[:3])
1928    elif Laue == '4/m m m':
1929        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1930        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1931        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])       
1932        HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],mat41[nxs,:,:])).T,HKLFT[:3])
1933        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1934    #trigonal - all hex cell
1935    #noncentrosymmetric - all ok
1936    elif Laue == '3':
1937        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1938        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1939        HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1940    elif Laue == '3 1 2':
1941        HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matk2[nxs,:,:])).T,HKLFT[:3])
1942        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1943        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1944        HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1945        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],matk2[nxs,:,:])).T,HKLFT[:3])
1946    elif Laue == '3 2 1':
1947        HKLFT[:3] = np.where(HKLFT[0]<=-2*HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1948        HKLFT[:3] = np.where(HKLFT[1]<-2*HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1949        HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1950        HKLFT[:3] = np.where((HKLFT[2]>0)&(HKLFT[1]==HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1951        HKLFT[:3] = np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T
1952        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])
1953    elif Laue == '3 1 m':
1954        HKLFT[:3] = np.where(HKLFT[0]>=HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1955        HKLFT[:3] = np.where(2*HKLFT[1]<-HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1956        HKLFT[:3] = np.where(HKLFT[1]>-2*HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matdmp[nxs,:,:])).T,HKLFT[:3])
1957        HKLFT[:3] = np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T
1958    elif (Laue == '3 m 1' or Laue == '3 m'):
1959        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1960        HKLFT[:3] = np.where((HKLFT[1]+HKLFT[0])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1961        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],matkm[nxs,:,:])).T,HKLFT[:3])
1962    #centrosymmetric
1963    elif Laue == '-3':  #ok
1964        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
1965        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1966        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1967        HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1968        HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[0]<0),-np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1969        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],-mat31[nxs,:,:])).T,HKLFT[:3])   
1970    elif (Laue == '-3 m 1' or Laue == '-3 m'):  #ok
1971        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1972        HKLFT[:3] = np.where((HKLFT[1]+HKLFT[0])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1973        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],matkm[nxs,:,:])).T,HKLFT[:3])
1974        HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1975        HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]<HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1976    elif Laue == '-3 1 m':  #ok
1977        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
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[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1980        HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1981        HKLFT[:3] = np.where(HKLFT[0]<=0,np.squeeze(np.inner(HKLF[:,:3],-mat31[nxs,:,:])).T,HKLFT[:3])   
1982        HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1983    #hexagonal
1984    #noncentrosymmetric
1985    elif Laue == '6':   #ok
1986        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1987        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1988        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1989        HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1990    elif Laue == '-6':  #ok
1991        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1992        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1993        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1994        HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1995    elif Laue == '6 2 2':   #ok
1996        HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1997        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1998        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1999        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
2000        HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
2001        HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[0]>HKLFT[1]),np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
2002    elif Laue == '6 m m':   #ok
2003        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
2004        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
2005        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[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[0]>HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
2008    elif Laue == '-6 m 2':  #ok
2009        HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matk2[nxs,:,:])).T,HKLFT[:3])
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[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
2013        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],matk2[nxs,:,:])).T,HKLFT[:3])
2014        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
2015    elif Laue == '-6 2 m':  #ok
2016        HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
2017        HKLFT[:3] = np.where(HKLFT[0]<=-2*HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
2018        HKLFT[:3] = np.where(HKLFT[1]<-2*HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
2019        HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
2020        HKLFT[:3] = np.where((HKLFT[2]>0)&(HKLFT[1]==HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
2021        HKLFT[:3] = np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T
2022        HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
2023        HKLFT[:3] = np.where(HKLFT[0]>HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
2024    #centrosymmetric
2025    elif Laue == '6/m': #ok
2026        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
2027        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
2028        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
2029        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
2030        HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
2031    elif Laue == '6/m m m': #ok
2032        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
2033        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
2034        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
2035        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
2036        HKLFT[:3] = np.where(HKLFT[0]>HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm.T[nxs,:,:])).T,HKLFT[:3])
2037    #cubic - all ok
2038    #noncentrosymmetric -
2039    elif Laue == '2 3': 
2040        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
2041        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
2042        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
2043        HKLFT[:3] = np.where((HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
2044        HKLFT[:3] = np.where((HKLFT[2]>=0)&((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
2045        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])
2046        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])
2047        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])
2048        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3])       
2049    elif Laue == '4 3 2':   
2050        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
2051        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
2052        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
2053        HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]<HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
2054        HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])   #in lieu od 2-fold
2055        HKLFT[:3] = np.where((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2]),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
2056        HKLFT[:3] = np.where((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2]),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
2057        HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat2d43[nxs,:,:])).T,HKLFT[:3])
2058    elif Laue == '-4 3 m': 
2059        HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])     
2060        HKLFT[:3] = np.where(HKLFT[0]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
2061        HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])     
2062        HKLFT[:3] = np.where(HKLFT[1]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
2063        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
2064        HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
2065        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
2066        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])
2067        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])
2068        HKLFT[:3] = np.where((HKLFT[2]>=0)&(HKLFT[1]<HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
2069        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3]) 
2070        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])
2071        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])
2072    #centrosymmetric
2073    elif Laue == 'm 3':
2074        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
2075        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],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[2]>=0)&((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
2078        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])
2079    elif Laue == 'm 3 m':
2080        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
2081        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
2082        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])           
2083        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])
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[0]>HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
2086    return HKLFT.T
2087       
2088
2089#Spherical harmonics routines
2090def OdfChk(SGLaue,L,M):
2091    'needs doc string'
2092    if not L%2 and abs(M) <= L:
2093        if SGLaue == '0':                      #cylindrical symmetry
2094            if M == 0: return True
2095        elif SGLaue == '-1':
2096            return True
2097        elif SGLaue == '2/m':
2098            if not abs(M)%2: return True
2099        elif SGLaue == 'mmm':
2100            if not abs(M)%2 and M >= 0: return True
2101        elif SGLaue == '4/m':
2102            if not abs(M)%4: return True
2103        elif SGLaue == '4/mmm':
2104            if not abs(M)%4 and M >= 0: return True
2105        elif SGLaue in ['3R','3']:
2106            if not abs(M)%3: return True
2107        elif SGLaue in ['3mR','3m1','31m']:
2108            if not abs(M)%3 and M >= 0: return True
2109        elif SGLaue == '6/m':
2110            if not abs(M)%6: return True
2111        elif SGLaue == '6/mmm':
2112            if not abs(M)%6 and M >= 0: return True
2113        elif SGLaue == 'm3':
2114            if M > 0:
2115                if L%12 == 2:
2116                    if M <= L//12: return True
2117                else:
2118                    if M <= L//12+1: return True
2119        elif SGLaue == 'm3m':
2120            if M > 0:
2121                if L%12 == 2:
2122                    if M <= L//12: return True
2123                else:
2124                    if M <= L//12+1: return True
2125    return False
2126       
2127def GenSHCoeff(SGLaue,SamSym,L,IfLMN=True):
2128    'needs doc string'
2129    coeffNames = []
2130    for iord in [2*i+2 for i in range(L//2)]:
2131        for m in [i-iord for i in range(2*iord+1)]:
2132            if OdfChk(SamSym,iord,m):
2133                for n in [i-iord for i in range(2*iord+1)]:
2134                    if OdfChk(SGLaue,iord,n):
2135                        if IfLMN:
2136                            coeffNames.append('C(%d,%d,%d)'%(iord,m,n))
2137                        else:
2138                            coeffNames.append('C(%d,%d)'%(iord,n))
2139    return coeffNames
2140   
2141def CrsAng(H,cell,SGData):
2142    'needs doc string'
2143    a,b,c,al,be,ga = cell
2144    SQ3 = 1.732050807569
2145    H1 = np.array([1,0,0])
2146    H2 = np.array([0,1,0])
2147    H3 = np.array([0,0,1])
2148    H4 = np.array([1,1,1])
2149    G,g = cell2Gmat(cell)
2150    Laue = SGData['SGLaue']
2151    Naxis = SGData['SGUniq']
2152    if len(H.shape) == 1:
2153        DH = np.inner(H,np.inner(G,H))
2154    else:
2155        DH = np.array([np.inner(h,np.inner(G,h)) for h in H])
2156    if Laue == '2/m':
2157        if Naxis == 'a':
2158            DR = np.inner(H1,np.inner(G,H1))
2159            DHR = np.inner(H,np.inner(G,H1))
2160        elif Naxis == 'b':
2161            DR = np.inner(H2,np.inner(G,H2))
2162            DHR = np.inner(H,np.inner(G,H2))
2163        else:
2164            DR = np.inner(H3,np.inner(G,H3))
2165            DHR = np.inner(H,np.inner(G,H3))
2166    elif Laue in ['R3','R3m']:
2167        DR = np.inner(H4,np.inner(G,H4))
2168        DHR = np.inner(H,np.inner(G,H4))
2169    else:
2170        DR = np.inner(H3,np.inner(G,H3))
2171        DHR = np.inner(H,np.inner(G,H3))
2172    DHR /= np.sqrt(DR*DH)
2173    phi = np.where(DHR <= 1.0,acosd(DHR),0.0)
2174    if Laue == '-1':
2175        BA = H.T[1]*a/(b-H.T[0]*cosd(ga))
2176        BB = H.T[0]*sind(ga)**2
2177    elif Laue == '2/m':
2178        if Naxis == 'a':
2179            BA = H.T[2]*b/(c-H.T[1]*cosd(al))
2180            BB = H.T[1]*sind(al)**2
2181        elif Naxis == 'b':
2182            BA = H.T[0]*c/(a-H.T[2]*cosd(be))
2183            BB = H.T[2]*sind(be)**2
2184        else:
2185            BA = H.T[1]*a/(b-H.T[0]*cosd(ga))
2186            BB = H.T[0]*sind(ga)**2
2187    elif Laue in ['mmm','4/m','4/mmm']:
2188        BA = H.T[1]*a
2189        BB = H.T[0]*b
2190    elif Laue in ['3R','3mR']:
2191        BA = H.T[0]+H.T[1]-2.0*H.T[2]
2192        BB = SQ3*(H.T[0]-H.T[1])
2193    elif Laue in ['m3','m3m']:
2194        BA = H.T[1]
2195        BB = H.T[0]
2196    else:
2197        BA = H.T[0]+2.0*H.T[1]
2198        BB = SQ3*H.T[0]
2199    beta = atan2d(BA,BB)
2200    return phi,beta
2201   
2202def SamAng(Tth,Gangls,Sangl,IFCoup):
2203    """Compute sample orientation angles vs laboratory coord. system
2204
2205    :param Tth:        Signed theta                                   
2206    :param Gangls:     Sample goniometer angles phi,chi,omega,azmuth 
2207    :param Sangl:      Sample angle zeros om-0, chi-0, phi-0         
2208    :param IFCoup:     True if omega & 2-theta coupled in CW scan
2209    :returns: 
2210        psi,gam:    Sample odf angles                             
2211        dPSdA,dGMdA:    Angle zero derivatives
2212    """                         
2213   
2214    if IFCoup:
2215        GSomeg = sind(Gangls[2]+Tth)
2216        GComeg = cosd(Gangls[2]+Tth)
2217    else:
2218        GSomeg = sind(Gangls[2])
2219        GComeg = cosd(Gangls[2])
2220    GSTth = sind(Tth)
2221    GCTth = cosd(Tth)     
2222    GSazm = sind(Gangls[3])
2223    GCazm = cosd(Gangls[3])
2224    GSchi = sind(Gangls[1])
2225    GCchi = cosd(Gangls[1])
2226    GSphi = sind(Gangls[0]+Sangl[2])
2227    GCphi = cosd(Gangls[0]+Sangl[2])
2228    SSomeg = sind(Sangl[0])
2229    SComeg = cosd(Sangl[0])
2230    SSchi = sind(Sangl[1])
2231    SCchi = cosd(Sangl[1])
2232    AT = -GSTth*GComeg+GCTth*GCazm*GSomeg
2233    BT = GSTth*GSomeg+GCTth*GCazm*GComeg
2234    CT = -GCTth*GSazm*GSchi
2235    DT = -GCTth*GSazm*GCchi
2236   
2237    BC1 = -AT*GSphi+(CT+BT*GCchi)*GCphi
2238    BC2 = DT-BT*GSchi
2239    BC3 = AT*GCphi+(CT+BT*GCchi)*GSphi
2240     
2241    BC = BC1*SComeg*SCchi+BC2*SComeg*SSchi-BC3*SSomeg     
2242    psi = acosd(BC)
2243   
2244    BD = 1.0-BC**2
2245    C = np.where(BD>1.e-6,rpd/np.sqrt(BD),0.)
2246    dPSdA = [-C*(-BC1*SSomeg*SCchi-BC2*SSomeg*SSchi-BC3*SComeg),
2247        -C*(-BC1*SComeg*SSchi+BC2*SComeg*SCchi),
2248        -C*(-BC1*SSomeg-BC3*SComeg*SCchi)]
2249     
2250    BA = -BC1*SSchi+BC2*SCchi
2251    BB = BC1*SSomeg*SCchi+BC2*SSomeg*SSchi+BC3*SComeg
2252    gam = atan2d(BB,BA)
2253
2254    BD = (BA**2+BB**2)/rpd
2255
2256    dBAdO = 0
2257    dBAdC = -BC1*SCchi-BC2*SSchi
2258    dBAdF = BC3*SSchi
2259   
2260    dBBdO = BC1*SComeg*SCchi+BC2*SComeg*SSchi-BC3*SSomeg
2261    dBBdC = -BC1*SSomeg*SSchi+BC2*SSomeg*SCchi
2262    dBBdF = BC1*SComeg-BC3*SSomeg*SCchi
2263   
2264    dGMdA = np.where(BD > 1.e-6,[(BA*dBBdO-BB*dBAdO)/BD,(BA*dBBdC-BB*dBAdC)/BD, \
2265        (BA*dBBdF-BB*dBAdF)/BD],[np.zeros_like(BD),np.zeros_like(BD),np.zeros_like(BD)])
2266       
2267    return psi,gam,dPSdA,dGMdA
2268
2269BOH = {
2270'L=2':[[],[],[]],
2271'L=4':[[0.30469720,0.36418281],[],[]],
2272'L=6':[[-0.14104740,0.52775103],[],[]],
2273'L=8':[[0.28646862,0.21545346,0.32826995],[],[]],
2274'L=10':[[-0.16413497,0.33078546,0.39371345],[],[]],
2275'L=12':[[0.26141975,0.27266871,0.03277460,0.32589402],
2276    [0.09298802,-0.23773812,0.49446631,0.0],[]],
2277'L=14':[[-0.17557309,0.25821932,0.27709173,0.33645360],[],[]],
2278'L=16':[[0.24370673,0.29873515,0.06447688,0.00377,0.32574495],
2279    [0.12039646,-0.25330128,0.23950998,0.40962508,0.0],[]],
2280'L=18':[[-0.16914245,0.17017340,0.34598142,0.07433932,0.32696037],
2281    [-0.06901768,0.16006562,-0.24743528,0.47110273,0.0],[]],
2282'L=20':[[0.23067026,0.31151832,0.09287682,0.01089683,0.00037564,0.32573563],
2283    [0.13615420,-0.25048007,0.12882081,0.28642879,0.34620433,0.0],[]],
2284'L=22':[[-0.16109560,0.10244188,0.36285175,0.13377513,0.01314399,0.32585583],
2285    [-0.09620055,0.20244115,-0.22389483,0.17928946,0.42017231,0.0],[]],
2286'L=24':[[0.22050742,0.31770654,0.11661736,0.02049853,0.00150861,0.00003426,0.32573505],
2287    [0.13651722,-0.21386648,0.00522051,0.33939435,0.10837396,0.32914497,0.0],
2288    [0.05378596,-0.11945819,0.16272298,-0.26449730,0.44923956,0.0,0.0]],
2289'L=26':[[-0.15435003,0.05261630,0.35524646,0.18578869,0.03259103,0.00186197,0.32574594],
2290    [-0.11306511,0.22072681,-0.18706142,0.05439948,0.28122966,0.35634355,0.0],[]],
2291'L=28':[[0.21225019,0.32031716,0.13604702,0.03132468,0.00362703,0.00018294,0.00000294,0.32573501],
2292    [0.13219496,-0.17206256,-0.08742608,0.32671661,0.17973107,0.02567515,0.32619598,0.0],
2293    [0.07989184,-0.16735346,0.18839770,-0.20705337,0.12926808,0.42715602,0.0,0.0]],
2294'L=30':[[-0.14878368,0.01524973,0.33628434,0.22632587,0.05790047,0.00609812,0.00022898,0.32573594],
2295    [-0.11721726,0.20915005,-0.11723436,-0.07815329,0.31318947,0.13655742,0.33241385,0.0],
2296    [-0.04297703,0.09317876,-0.11831248,0.17355132,-0.28164031,0.42719361,0.0,0.0]],
2297'L=32':[[0.20533892,0.32087437,0.15187897,0.04249238,0.00670516,0.00054977,0.00002018,0.00000024,0.32573501],
2298    [0.12775091,-0.13523423,-0.14935701,0.28227378,0.23670434,0.05661270,0.00469819,0.32578978,0.0],
2299    [0.09703829,-0.19373733,0.18610682,-0.14407046,0.00220535,0.26897090,0.36633402,0.0,0.0]],
2300'L=34':[[-0.14409234,-0.01343681,0.31248977,0.25557722,0.08571889,0.01351208,0.00095792,0.00002550,0.32573508],
2301    [-0.11527834,0.18472133,-0.04403280,-0.16908618,0.27227021,0.21086614,0.04041752,0.32688152,0.0],
2302    [-0.06773139,0.14120811,-0.15835721,0.18357456,-0.19364673,0.08377174,0.43116318,0.0,0.0]]
2303}
2304
2305Lnorm = lambda L: 4.*np.pi/(2.0*L+1.)
2306
2307def GetKcl(L,N,SGLaue,phi,beta):
2308    'needs doc string'
2309    import pytexture as ptx
2310    if SGLaue in ['m3','m3m']:
2311        if 'array' in str(type(phi)) and np.any(phi.shape):
2312            Kcl = np.zeros_like(phi)
2313        else:
2314            Kcl = 0.
2315        for j in range(0,L+1,4):
2316            im = j//4
2317            if 'array' in str(type(phi)) and np.any(phi.shape):
2318                pcrs = ptx.pyplmpsi(L,j,len(phi),phi)[0]
2319            else:
2320                pcrs = ptx.pyplmpsi(L,j,1,phi)[0]
2321            Kcl += BOH['L=%d'%(L)][N-1][im]*pcrs*cosd(j*beta)       
2322    else:
2323        if 'array' in str(type(phi)) and np.any(phi.shape):
2324            pcrs = ptx.pyplmpsi(L,N,len(phi),phi)[0]
2325        else:
2326            pcrs = ptx.pyplmpsi(L,N,1,phi)[0]
2327        pcrs *= RSQ2PI
2328        if N:
2329            pcrs *= SQ2
2330        if SGLaue in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']:
2331            if SGLaue in ['3mR','3m1','31m']: 
2332                if N%6 == 3:
2333                    Kcl = pcrs*sind(N*beta)
2334                else:
2335                    Kcl = pcrs*cosd(N*beta)
2336            else:
2337                Kcl = pcrs*cosd(N*beta)
2338        else:
2339            Kcl = pcrs*(cosd(N*beta)+sind(N*beta))
2340    return Kcl
2341   
2342def GetKsl(L,M,SamSym,psi,gam):
2343    'needs doc string'
2344    import pytexture as ptx
2345    if 'array' in str(type(psi)) and np.any(psi.shape):
2346        psrs,dpdps = ptx.pyplmpsi(L,M,len(psi),psi)
2347    else:
2348        psrs,dpdps = ptx.pyplmpsi(L,M,1,psi)
2349    psrs *= RSQ2PI
2350    dpdps *= RSQ2PI
2351    if M:
2352        psrs *= SQ2
2353        dpdps *= SQ2
2354    if SamSym in ['mmm',]:
2355        dum = cosd(M*gam)
2356        Ksl = psrs*dum
2357        dKsdp = dpdps*dum
2358        dKsdg = -psrs*M*sind(M*gam)
2359    else:
2360        dum = cosd(M*gam)+sind(M*gam)
2361        Ksl = psrs*dum
2362        dKsdp = dpdps*dum
2363        dKsdg = psrs*M*(-sind(M*gam)+cosd(M*gam))
2364    return Ksl,dKsdp,dKsdg
2365   
2366def GetKclKsl(L,N,SGLaue,psi,phi,beta):
2367    """
2368    This is used for spherical harmonics description of preferred orientation;
2369        cylindrical symmetry only (M=0) and no sample angle derivatives returned
2370    """
2371    import pytexture as ptx
2372    Ksl,x = ptx.pyplmpsi(L,0,1,psi)
2373    Ksl *= RSQ2PI
2374    if SGLaue in ['m3','m3m']:
2375        Kcl = 0.0
2376        for j in range(0,L+1,4):
2377            im = j//4
2378            pcrs,dum = ptx.pyplmpsi(L,j,1,phi)
2379            Kcl += BOH['L=%d'%(L)][N-1][im]*pcrs*cosd(j*beta)       
2380    else:
2381        pcrs,dum = ptx.pyplmpsi(L,N,1,phi)
2382        pcrs *= RSQ2PI
2383        if N:
2384            pcrs *= SQ2
2385        if SGLaue in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']:
2386            if SGLaue in ['3mR','3m1','31m']: 
2387                if N%6 == 3:
2388                    Kcl = pcrs*sind(N*beta)
2389                else:
2390                    Kcl = pcrs*cosd(N*beta)
2391            else:
2392                Kcl = pcrs*cosd(N*beta)
2393        else:
2394            Kcl = pcrs*(cosd(N*beta)+sind(N*beta))
2395    return Kcl*Ksl,Lnorm(L)
2396   
2397def Glnh(Start,SHCoef,psi,gam,SamSym):
2398    'needs doc string'
2399    import pytexture as ptx
2400
2401    if Start:
2402        ptx.pyqlmninit()
2403        Start = False
2404    Fln = np.zeros(len(SHCoef))
2405    for i,term in enumerate(SHCoef):
2406        l,m,n = eval(term.strip('C'))
2407        pcrs,dum = ptx.pyplmpsi(l,m,1,psi)
2408        pcrs *= RSQPI
2409        if m == 0:
2410            pcrs /= SQ2
2411        if SamSym in ['mmm',]:
2412            Ksl = pcrs*cosd(m*gam)
2413        else:
2414            Ksl = pcrs*(cosd(m*gam)+sind(m*gam))
2415        Fln[i] = SHCoef[term]*Ksl*Lnorm(l)
2416    ODFln = dict(zip(SHCoef.keys(),list(zip(SHCoef.values(),Fln))))
2417    return ODFln
2418
2419def Flnh(Start,SHCoef,phi,beta,SGData):
2420    'needs doc string'
2421    import pytexture as ptx
2422   
2423    if Start:
2424        ptx.pyqlmninit()
2425        Start = False
2426    Fln = np.zeros(len(SHCoef))
2427    for i,term in enumerate(SHCoef):
2428        l,m,n = eval(term.strip('C'))
2429        if SGData['SGLaue'] in ['m3','m3m']:
2430            Kcl = 0.0
2431            for j in range(0,l+1,4):
2432                im = j//4
2433                pcrs,dum = ptx.pyplmpsi(l,j,1,phi)
2434                Kcl += BOH['L='+str(l)][n-1][im]*pcrs*cosd(j*beta)       
2435        else:                #all but cubic
2436            pcrs,dum = ptx.pyplmpsi(l,n,1,phi)
2437            pcrs *= RSQPI
2438            if n == 0:
2439                pcrs /= SQ2
2440            if SGData['SGLaue'] in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']:
2441               if SGData['SGLaue'] in ['3mR','3m1','31m']: 
2442                   if n%6 == 3:
2443                       Kcl = pcrs*sind(n*beta)
2444                   else:
2445                       Kcl = pcrs*cosd(n*beta)
2446               else:
2447                   Kcl = pcrs*cosd(n*beta)
2448            else:
2449                Kcl = pcrs*(cosd(n*beta)+sind(n*beta))
2450        Fln[i] = SHCoef[term]*Kcl*Lnorm(l)
2451    ODFln = dict(zip(SHCoef.keys(),list(zip(SHCoef.values(),Fln))))
2452    return ODFln
2453   
2454def polfcal(ODFln,SamSym,psi,gam):
2455    '''Perform a pole figure computation.
2456    Note that the the number of gam values must either be 1 or must
2457    match psi. Updated for numpy 1.8.0
2458    '''
2459    import pytexture as ptx
2460    PolVal = np.ones_like(psi)
2461    for term in ODFln:
2462        if abs(ODFln[term][1]) > 1.e-3:
2463            l,m,n = eval(term.strip('C'))
2464            psrs,dum = ptx.pyplmpsi(l,m,len(psi),psi)
2465            if SamSym in ['-1','2/m']:
2466                if m:
2467                    Ksl = RSQPI*psrs*(cosd(m*gam)+sind(m*gam))
2468                else:
2469                    Ksl = RSQPI*psrs/SQ2
2470            else:
2471                if m:
2472                    Ksl = RSQPI*psrs*cosd(m*gam)
2473                else:
2474                    Ksl = RSQPI*psrs/SQ2
2475            PolVal += ODFln[term][1]*Ksl
2476    return PolVal
2477   
2478def invpolfcal(ODFln,SGData,phi,beta):
2479    'needs doc string'
2480    import pytexture as ptx
2481   
2482    invPolVal = np.ones_like(beta)
2483    for term in ODFln:
2484        if abs(ODFln[term][1]) > 1.e-3:
2485            l,m,n = eval(term.strip('C'))
2486            if SGData['SGLaue'] in ['m3','m3m']:
2487                Kcl = 0.0
2488                for j in range(0,l+1,4):
2489                    im = j//4
2490                    pcrs,dum = ptx.pyplmpsi(l,j,len(beta),phi)
2491                    Kcl += BOH['L=%d'%(l)][n-1][im]*pcrs*cosd(j*beta)       
2492            else:                #all but cubic
2493                pcrs,dum = ptx.pyplmpsi(l,n,len(beta),phi)
2494                pcrs *= RSQPI
2495                if n == 0:
2496                    pcrs /= SQ2
2497                if SGData['SGLaue'] in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']:
2498                   if SGData['SGLaue'] in ['3mR','3m1','31m']: 
2499                       if n%6 == 3:
2500                           Kcl = pcrs*sind(n*beta)
2501                       else:
2502                           Kcl = pcrs*cosd(n*beta)
2503                   else:
2504                       Kcl = pcrs*cosd(n*beta)
2505                else:
2506                    Kcl = pcrs*(cosd(n*beta)+sind(n*beta))
2507            invPolVal += ODFln[term][1]*Kcl
2508    return invPolVal
2509   
2510   
2511def textureIndex(SHCoef):
2512    'needs doc string'
2513    Tindx = 1.0
2514    for term in SHCoef:
2515        l = eval(term.strip('C'))[0]
2516        Tindx += SHCoef[term]**2/(2.0*l+1.)
2517    return Tindx
2518
2519UniqueCellByLaue = [
2520        [['m3','m3m'],(0,)],
2521        [['3R','3mR'],(0,3)],
2522        [['3','3m1','31m','6/m','6/mmm','4/m','4/mmm'],(0,2)],
2523        [['mmm'],(0,1,2)],
2524        [['2/m'+'a'],(0,1,2,3)],
2525        [['2/m'+'b'],(0,1,2,4)],
2526        [['2/m'+'c'],(0,1,2,5)],
2527        [['-1'],(0,1,2,3,4,5)],
2528    ]
2529'''List the unique cell terms by index for each Laue class'''
2530
2531cellAlbl = ('a','b','c', 'alpha', 'beta', 'gamma')
2532'ASCII labels for a, b, c, alpha, beta, gamma'
2533
2534cellUlbl = ('a','b','c',u'\u03B1',u'\u03B2',u'\u03B3')
2535'unicode labels for a, b, c, alpha, beta, gamma'
2536 
2537# self-test materials follow.
2538selftestlist = []
2539'''Defines a list of self-tests'''
2540selftestquiet = True
2541def _ReportTest():
2542    'Report name and doc string of current routine when ``selftestquiet`` is False'
2543    if not selftestquiet:
2544        import inspect
2545        caller = inspect.stack()[1][3]
2546        doc = eval(caller).__doc__
2547        if doc is not None:
2548            print('testing '+__file__+' with '+caller+' ('+doc+')')
2549        else:
2550            print('testing '+__file__()+" with "+caller)
2551NeedTestData = True
2552def TestData():
2553    array = np.array
2554    global NeedTestData
2555    NeedTestData = False
2556    global CellTestData
2557    # output from uctbx computed on platform darwin on 2010-05-28
2558    CellTestData = [
2559# cell, g, G, cell*, V, V*
2560  [(4, 4, 4, 90, 90, 90), 
2561   array([[  1.60000000e+01,   9.79717439e-16,   9.79717439e-16],
2562       [  9.79717439e-16,   1.60000000e+01,   9.79717439e-16],
2563       [  9.79717439e-16,   9.79717439e-16,   1.60000000e+01]]), array([[  6.25000000e-02,   3.82702125e-18,   3.82702125e-18],
2564       [  3.82702125e-18,   6.25000000e-02,   3.82702125e-18],
2565       [  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],
2566# cell, g, G, cell*, V, V*
2567  [(4.0999999999999996, 5.2000000000000002, 6.2999999999999998, 100, 80, 130), 
2568   array([[ 16.81      , -13.70423184,   4.48533243],
2569       [-13.70423184,  27.04      ,  -5.6887143 ],
2570       [  4.48533243,  -5.6887143 ,  39.69      ]]), array([[ 0.10206349,  0.05083339, -0.00424823],
2571       [ 0.05083339,  0.06344997,  0.00334956],
2572       [-0.00424823,  0.00334956,  0.02615544]]), (0.31947376387537696, 0.25189277536327803, 0.16172643497798223, 85.283666420376008, 94.716333579624006, 50.825714168082683), 100.98576357983838, 0.0099023858863968445],
2573# cell, g, G, cell*, V, V*
2574  [(3.5, 3.5, 6, 90, 90, 120), 
2575   array([[  1.22500000e+01,  -6.12500000e+00,   1.28587914e-15],
2576       [ -6.12500000e+00,   1.22500000e+01,   1.28587914e-15],
2577       [  1.28587914e-15,   1.28587914e-15,   3.60000000e+01]]), array([[  1.08843537e-01,   5.44217687e-02,   3.36690552e-18],
2578       [  5.44217687e-02,   1.08843537e-01,   3.36690552e-18],
2579       [  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],
2580  ]
2581    global CoordTestData
2582    CoordTestData = [
2583# cell, ((frac, ortho),...)
2584  ((4,4,4,90,90,90,), [
2585 ((0.10000000000000001, 0.0, 0.0),(0.40000000000000002, 0.0, 0.0)),
2586 ((0.0, 0.10000000000000001, 0.0),(2.4492935982947065e-17, 0.40000000000000002, 0.0)),
2587 ((0.0, 0.0, 0.10000000000000001),(2.4492935982947065e-17, -2.4492935982947065e-17, 0.40000000000000002)),
2588 ((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(0.40000000000000013, 0.79999999999999993, 1.2)),
2589 ((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(0.80000000000000016, 1.2, 0.40000000000000002)),
2590 ((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(1.2, 0.80000000000000004, 0.40000000000000002)),
2591 ((0.5, 0.5, 0.5),(2.0, 1.9999999999999998, 2.0)),
2592]),
2593# cell, ((frac, ortho),...)
2594  ((4.1,5.2,6.3,100,80,130,), [
2595 ((0.10000000000000001, 0.0, 0.0),(0.40999999999999998, 0.0, 0.0)),
2596 ((0.0, 0.10000000000000001, 0.0),(-0.33424955703700043, 0.39834311042186865, 0.0)),
2597 ((0.0, 0.0, 0.10000000000000001),(0.10939835193016617, -0.051013289294572106, 0.6183281045774256)),
2598 ((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(0.069695941716497567, 0.64364635296002093, 1.8549843137322766)),
2599 ((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(-0.073350319180835066, 1.1440160419710339, 0.6183281045774256)),
2600 ((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(0.67089923785616512, 0.74567293154916525, 0.6183281045774256)),
2601 ((0.5, 0.5, 0.5),(0.92574397446582857, 1.7366491056364828, 3.0916405228871278)),
2602]),
2603# cell, ((frac, ortho),...)
2604  ((3.5,3.5,6,90,90,120,), [
2605 ((0.10000000000000001, 0.0, 0.0),(0.35000000000000003, 0.0, 0.0)),
2606 ((0.0, 0.10000000000000001, 0.0),(-0.17499999999999993, 0.3031088913245536, 0.0)),
2607 ((0.0, 0.0, 0.10000000000000001),(3.6739403974420595e-17, -3.6739403974420595e-17, 0.60000000000000009)),
2608 ((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(2.7675166561703527e-16, 0.60621778264910708, 1.7999999999999998)),
2609 ((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(0.17500000000000041, 0.90932667397366063, 0.60000000000000009)),
2610 ((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(0.70000000000000018, 0.6062177826491072, 0.60000000000000009)),
2611 ((0.5, 0.5, 0.5),(0.87500000000000067, 1.5155444566227676, 3.0)),
2612]),
2613]
2614    global LaueTestData             #generated by GSAS
2615    LaueTestData = {
2616    '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),
2617        (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),
2618        (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))],
2619    '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),
2620        (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),
2621        (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),
2622        (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))],
2623    '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),
2624        (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),
2625        (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),
2626        (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),
2627        (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),
2628        (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),
2629        (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),
2630        (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),
2631        (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),
2632        (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),
2633        (2,1,5,6),(2,1,-5,6),(3,-1,-5,6))],
2634    '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),
2635        (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),
2636        (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),
2637        (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),
2638        (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),
2639        (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),
2640        (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),
2641        (3,0,4,6),(3,1,-2,12),(3,0,-4,6),(1,1,6,12),(2,2,3,12))],
2642    '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),
2643        (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),
2644        (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),
2645        (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),
2646        (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),
2647        (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),
2648        (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))],
2649    }
2650   
2651    global FLnhTestData
2652    FLnhTestData = [{
2653    'C(4,0,0)': (0.965, 0.42760447),
2654    'C(2,0,0)': (1.0122, -0.80233610),
2655    'C(2,0,2)': (0.0061, 8.37491546E-03),
2656    'C(6,0,4)': (-0.0898, 4.37985696E-02),
2657    'C(6,0,6)': (-0.1369, -9.04081762E-02),
2658    'C(6,0,0)': (0.5935, -0.18234928),
2659    'C(4,0,4)': (0.1872, 0.16358127),
2660    'C(6,0,2)': (0.6193, 0.27573633),
2661    'C(4,0,2)': (-0.1897, 0.12530720)},[1,0,0]]
2662def test0():
2663    if NeedTestData: TestData()
2664    msg = 'test cell2Gmat, fillgmat, Gmat2cell'
2665    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2666        G, g = cell2Gmat(cell)
2667        assert np.allclose(G,tG),msg
2668        assert np.allclose(g,tg),msg
2669        tcell = Gmat2cell(g)
2670        assert np.allclose(cell,tcell),msg
2671        tcell = Gmat2cell(G)
2672        assert np.allclose(tcell,trcell),msg
2673if __name__ == '__main__': selftestlist.append(test0)
2674
2675def test1():
2676    'test cell2A and A2Gmat'
2677    _ReportTest()
2678    if NeedTestData: TestData()
2679    msg = 'test cell2A and A2Gmat'
2680    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2681        G, g = A2Gmat(cell2A(cell))
2682        assert np.allclose(G,tG),msg
2683        assert np.allclose(g,tg),msg
2684if __name__ == '__main__': selftestlist.append(test1)
2685
2686def test2():
2687    'test Gmat2A, A2cell, A2Gmat, Gmat2cell'
2688    _ReportTest()
2689    if NeedTestData: TestData()
2690    msg = 'test Gmat2A, A2cell, A2Gmat, Gmat2cell'
2691    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2692        G, g = cell2Gmat(cell)
2693        tcell = A2cell(Gmat2A(G))
2694        assert np.allclose(cell,tcell),msg
2695if __name__ == '__main__': selftestlist.append(test2)
2696
2697def test3():
2698    'test invcell2Gmat'
2699    _ReportTest()
2700    if NeedTestData: TestData()
2701    msg = 'test invcell2Gmat'
2702    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2703        G, g = invcell2Gmat(trcell)
2704        assert np.allclose(G,tG),msg
2705        assert np.allclose(g,tg),msg
2706if __name__ == '__main__': selftestlist.append(test3)
2707
2708def test4():
2709    'test calc_rVsq, calc_rV, calc_V'
2710    _ReportTest()
2711    if NeedTestData: TestData()
2712    msg = 'test calc_rVsq, calc_rV, calc_V'
2713    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2714        assert np.allclose(calc_rV(cell2A(cell)),trV), msg
2715        assert np.allclose(calc_V(cell2A(cell)),tV), msg
2716if __name__ == '__main__': selftestlist.append(test4)
2717
2718def test5():
2719    'test A2invcell'
2720    _ReportTest()
2721    if NeedTestData: TestData()
2722    msg = 'test A2invcell'
2723    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2724        rcell = A2invcell(cell2A(cell))
2725        assert np.allclose(rcell,trcell),msg
2726if __name__ == '__main__': selftestlist.append(test5)
2727
2728def test6():
2729    'test cell2AB'
2730    _ReportTest()
2731    if NeedTestData: TestData()
2732    msg = 'test cell2AB'
2733    for (cell,coordlist) in CoordTestData:
2734        A,B = cell2AB(cell)
2735        for (frac,ortho) in coordlist:
2736            to = np.inner(A,frac)
2737            tf = np.inner(B,to)
2738            assert np.allclose(ortho,to), msg
2739            assert np.allclose(frac,tf), msg
2740            to = np.sum(A*frac,axis=1)
2741            tf = np.sum(B*to,axis=1)
2742            assert np.allclose(ortho,to), msg
2743            assert np.allclose(frac,tf), msg
2744if __name__ == '__main__': selftestlist.append(test6)
2745
2746def test7():
2747    'test GetBraviasNum(...) and GenHBravais(...)'
2748    _ReportTest()
2749    import os.path
2750    import sys
2751    import GSASIIspc as spc
2752    testdir = os.path.join(os.path.split(os.path.abspath( __file__ ))[0],'testinp')
2753    if os.path.exists(testdir):
2754        if testdir not in sys.path: sys.path.insert(0,testdir)
2755    import sgtbxlattinp
2756    derror = 1e-4
2757    def indexmatch(hklin, hkllist, system):
2758        for hklref in hkllist:
2759            hklref = list(hklref)
2760            # these permutations are far from complete, but are sufficient to
2761            # allow the test to complete
2762            if system == 'cubic':
2763                permlist = [(1,2,3),(1,3,2),(2,1,3),(2,3,1),(3,1,2),(3,2,1),]
2764            elif system == 'monoclinic':
2765                permlist = [(1,2,3),(-1,2,-3)]
2766            else:
2767                permlist = [(1,2,3)]
2768
2769            for perm in permlist:
2770                hkl = [abs(i) * hklin[abs(i)-1] / i for i in perm]
2771                if hkl == hklref: return True
2772                if [-i for i in hkl] == hklref: return True
2773        else:
2774            return False
2775
2776    for key in sgtbxlattinp.sgtbx7:
2777        spdict = spc.SpcGroup(key)
2778        cell = sgtbxlattinp.sgtbx7[key][0]
2779        system = spdict[1]['SGSys']
2780        center = spdict[1]['SGLatt']
2781
2782        bravcode = GetBraviasNum(center, system)
2783
2784        g2list = GenHBravais(sgtbxlattinp.dmin, bravcode, cell2A(cell))
2785
2786        assert len(sgtbxlattinp.sgtbx7[key][1]) == len(g2list), 'Reflection lists differ for %s' % key
2787        for h,k,l,d,num in g2list:
2788            for hkllist,dref in sgtbxlattinp.sgtbx7[key][1]: 
2789                if abs(d-dref) < derror:
2790                    if indexmatch((h,k,l,), hkllist, system):
2791                        break
2792            else:
2793                assert 0,'No match for %s at %s (%s)' % ((h,k,l),d,key)
2794if __name__ == '__main__': selftestlist.append(test7)
2795
2796def test8():
2797    'test GenHLaue'
2798    _ReportTest()
2799    import GSASIIspc as spc
2800    import sgtbxlattinp
2801    derror = 1e-4
2802    dmin = sgtbxlattinp.dmin
2803
2804    def indexmatch(hklin, hklref, system, axis):
2805        # these permutations are far from complete, but are sufficient to
2806        # allow the test to complete
2807        if system == 'cubic':
2808            permlist = [(1,2,3),(1,3,2),(2,1,3),(2,3,1),(3,1,2),(3,2,1),]
2809        elif system == 'monoclinic' and axis=='b':
2810            permlist = [(1,2,3),(-1,2,-3)]
2811        elif system == 'monoclinic' and axis=='a':
2812            permlist = [(1,2,3),(1,-2,-3)]
2813        elif system == 'monoclinic' and axis=='c':
2814            permlist = [(1,2,3),(-1,-2,3)]
2815        elif system == 'trigonal':
2816            permlist = [(1,2,3),(2,1,3),(-1,-2,3),(-2,-1,3)]
2817        elif system == 'rhombohedral':
2818            permlist = [(1,2,3),(2,3,1),(3,1,2)]
2819        else:
2820            permlist = [(1,2,3)]
2821
2822        hklref = list(hklref)
2823        for perm in permlist:
2824            hkl = [abs(i) * hklin[abs(i)-1] / i for i in perm]
2825            if hkl == hklref: return True
2826            if [-i for i in hkl] == hklref: return True
2827        return False
2828
2829    for key in sgtbxlattinp.sgtbx8:
2830        spdict = spc.SpcGroup(key)[1]
2831        cell = sgtbxlattinp.sgtbx8[key][0]
2832        Axis = spdict['SGUniq']
2833        system = spdict['SGSys']
2834
2835        g2list = GenHLaue(dmin,spdict,cell2A(cell))
2836        #if len(g2list) != len(sgtbxlattinp.sgtbx8[key][1]):
2837        #    print 'failed',key,':' ,len(g2list),'vs',len(sgtbxlattinp.sgtbx8[key][1])
2838        #    print 'GSAS-II:'
2839        #    for h,k,l,d in g2list: print '  ',(h,k,l),d
2840        #    print 'SGTBX:'
2841        #    for hkllist,dref in sgtbxlattinp.sgtbx8[key][1]: print '  ',hkllist,dref
2842        assert len(g2list) == len(sgtbxlattinp.sgtbx8[key][1]), (
2843            'Reflection lists differ for %s' % key
2844            )
2845        #match = True
2846        for h,k,l,d in g2list:
2847            for hkllist,dref in sgtbxlattinp.sgtbx8[key][1]: 
2848                if abs(d-dref) < derror:
2849                    if indexmatch((h,k,l,), hkllist, system, Axis): break
2850            else:
2851                assert 0,'No match for %s at %s (%s)' % ((h,k,l),d,key)
2852                #match = False
2853        #if not match:
2854            #for hkllist,dref in sgtbxlattinp.sgtbx8[key][1]: print '  ',hkllist,dref
2855            #print center, Laue, Axis, system
2856if __name__ == '__main__': selftestlist.append(test8)
2857           
2858def test9():
2859    'test GenHLaue'
2860    _ReportTest()
2861    import GSASIIspc as G2spc
2862    if NeedTestData: TestData()
2863    for spc in LaueTestData:
2864        data = LaueTestData[spc]
2865        cell = data[0]
2866        hklm = np.array(data[1])
2867        H = hklm[-1][:3]
2868        hklO = hklm.T[:3].T
2869        A = cell2A(cell)
2870        dmin = 1./np.sqrt(calc_rDsq(H,A))
2871        SGData = G2spc.SpcGroup(spc)[1]
2872        hkls = np.array(GenHLaue(dmin,SGData,A))
2873        hklN = hkls.T[:3].T
2874        #print spc,hklO.shape,hklN.shape
2875        err = True
2876        for H in hklO:
2877            if H not in hklN:
2878                print ('%d %s'%(H,' missing from hkl from GSASII'))
2879                err = False
2880        assert(err)
2881if __name__ == '__main__': selftestlist.append(test9)
2882       
2883       
2884   
2885
2886if __name__ == '__main__':
2887    # run self-tests
2888    selftestquiet = False
2889    for test in selftestlist:
2890        test()
2891    print ("OK")
Note: See TracBrowser for help on using the repository browser.