source: trunk/GSASIIlattice.py @ 5147

Last change on this file since 5147 was 5147, checked in by vondreele, 11 months ago

Add new routine to G2lat - AplusDij? & use it in a couple of places
add a new column to the ISODISTORT mode display - shows values obtained from ISODISTORT
remove calls to GetHstrainShift? for all powder data cases - wrong thing to do.

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