source: trunk/GSASIIlattice.py @ 5352

Last change on this file since 5352 was 5162, checked in by vondreele, 3 years ago

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