source: trunk/GSASIIlattice.py @ 5486

Last change on this file since 5486 was 5478, checked in by vondreele, 2 years ago

fixes to RMCProfile dat file prep
more updates to sp. harm stuff

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