source: trunk/GSASIIlattice.py @ 4998

Last change on this file since 4998 was 4998, checked in by toby, 3 months ago

partial fix of magnetic constraints (still need to deal w/unvaried items; fullrmc minor stuff

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