source: trunk/GSASIIlattice.py @ 1367

Last change on this file since 1367 was 1367, checked in by vondreele, 8 years ago

add Dsp2pos & Pos2dsp for conversions to G2lattice
remove wave from DoIndexPeaks? - wasn't needed
use Pos2dsp & Dsp2pos in various places in G2pwdGUI
set sort order for TOF peaks - makes indexing work!

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 58.0 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_1 h^2 + A_2 k^2 + A_3 l^2 + A_4 hk + A_5 hl + A_6 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`
24'''
25########### SVN repository information ###################
26# $Date: 2014-05-29 14:52:18 +0000 (Thu, 29 May 2014) $
27# $Author: vondreele $
28# $Revision: 1367 $
29# $URL: trunk/GSASIIlattice.py $
30# $Id: GSASIIlattice.py 1367 2014-05-29 14:52:18Z vondreele $
31########### SVN repository information ###################
32import math
33import numpy as np
34import numpy.linalg as nl
35import GSASIIpath
36GSASIIpath.SetVersionNumber("$Revision: 1367 $")
37# trig functions in degrees
38sind = lambda x: np.sin(x*np.pi/180.)
39asind = lambda x: 180.*np.arcsin(x)/np.pi
40tand = lambda x: np.tan(x*np.pi/180.)
41atand = lambda x: 180.*np.arctan(x)/np.pi
42atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
43cosd = lambda x: np.cos(x*np.pi/180.)
44acosd = lambda x: 180.*np.arccos(x)/np.pi
45rdsq2d = lambda x,p: round(1.0/np.sqrt(x),p)
46
47def sec2HMS(sec):
48    """Convert time in sec to H:M:S string
49   
50    :param sec: time in seconds
51    :return: H:M:S string (to nearest 100th second)
52   
53    """
54    H = int(sec/3600)
55    M = int(sec/60-H*60)
56    S = sec-3600*H-60*M
57    return '%d:%2d:%.2f'%(H,M,S)
58   
59def rotdMat(angle,axis=0):
60    """Prepare rotation matrix for angle in degrees about axis(=0,1,2)
61
62    :param angle: angle in degrees
63    :param axis:  axis (0,1,2 = x,y,z) about which for the rotation
64    :return: rotation matrix - 3x3 numpy array
65
66    """
67    if axis == 2:
68        return np.array([[cosd(angle),-sind(angle),0],[sind(angle),cosd(angle),0],[0,0,1]])
69    elif axis == 1:
70        return np.array([[cosd(angle),0,-sind(angle)],[0,1,0],[sind(angle),0,cosd(angle)]])
71    else:
72        return np.array([[1,0,0],[0,cosd(angle),-sind(angle)],[0,sind(angle),cosd(angle)]])
73       
74def rotdMat4(angle,axis=0):
75    """Prepare rotation matrix for angle in degrees about axis(=0,1,2) with scaling for OpenGL
76
77    :param angle: angle in degrees
78    :param axis:  axis (0,1,2 = x,y,z) about which for the rotation
79    :return: rotation matrix - 4x4 numpy array (last row/column for openGL scaling)
80
81    """
82    Mat = rotdMat(angle,axis)
83    return np.concatenate((np.concatenate((Mat,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0)
84   
85def fillgmat(cell):
86    """Compute lattice metric tensor from unit cell constants
87
88    :param cell: tuple with a,b,c,alpha, beta, gamma (degrees)
89    :return: 3x3 numpy array
90
91    """
92    a,b,c,alp,bet,gam = cell
93    g = np.array([
94        [a*a,  a*b*cosd(gam),  a*c*cosd(bet)],
95        [a*b*cosd(gam),  b*b,  b*c*cosd(alp)],
96        [a*c*cosd(bet) ,b*c*cosd(alp),   c*c]])
97    return g
98           
99def cell2Gmat(cell):
100    """Compute real and reciprocal lattice metric tensor from unit cell constants
101
102    :param cell: tuple with a,b,c,alpha, beta, gamma (degrees)
103    :return: reciprocal (G) & real (g) metric tensors (list of two numpy 3x3 arrays)
104
105    """
106    g = fillgmat(cell)
107    G = nl.inv(g)       
108    return G,g
109
110def A2Gmat(A,inverse=True):
111    """Fill real & reciprocal metric tensor (G) from A.
112
113    :param A: reciprocal metric tensor elements as [G11,G22,G33,2*G12,2*G13,2*G23]
114    :param bool inverse: if True return both G and g; else just G
115    :return: reciprocal (G) & real (g) metric tensors (list of two numpy 3x3 arrays)
116
117    """
118    G = np.zeros(shape=(3,3))
119    G = [
120        [A[0],  A[3]/2.,  A[4]/2.], 
121        [A[3]/2.,A[1],    A[5]/2.], 
122        [A[4]/2.,A[5]/2.,    A[2]]]
123    if inverse:
124        g = nl.inv(G)
125        return G,g
126    else:
127        return G
128
129def Gmat2A(G):
130    """Extract A from reciprocal metric tensor (G)
131
132    :param G: reciprocal maetric tensor (3x3 numpy array
133    :return: A = [G11,G22,G33,2*G12,2*G13,2*G23]
134
135    """
136    return [G[0][0],G[1][1],G[2][2],2.*G[0][1],2.*G[0][2],2.*G[1][2]]
137   
138def cell2A(cell):
139    """Obtain A = [G11,G22,G33,2*G12,2*G13,2*G23] from lattice parameters
140
141    :param cell: [a,b,c,alpha,beta,gamma] (degrees)
142    :return: G reciprocal metric tensor as 3x3 numpy array
143
144    """
145    G,g = cell2Gmat(cell)
146    return Gmat2A(G)
147
148def A2cell(A):
149    """Compute unit cell constants from A
150
151    :param A: [G11,G22,G33,2*G12,2*G13,2*G23] G - reciprocal metric tensor
152    :return: a,b,c,alpha, beta, gamma (degrees) - lattice parameters
153
154    """
155    G,g = A2Gmat(A)
156    return Gmat2cell(g)
157
158def Gmat2cell(g):
159    """Compute real/reciprocal lattice parameters from real/reciprocal metric tensor (g/G)
160    The math works the same either way.
161
162    :param g (or G): real (or reciprocal) metric tensor 3x3 array
163    :return: a,b,c,alpha, beta, gamma (degrees) (or a*,b*,c*,alpha*,beta*,gamma* degrees)
164
165    """
166    oldset = np.seterr('raise')
167    a = np.sqrt(max(0,g[0][0]))
168    b = np.sqrt(max(0,g[1][1]))
169    c = np.sqrt(max(0,g[2][2]))
170    alp = acosd(g[2][1]/(b*c))
171    bet = acosd(g[2][0]/(a*c))
172    gam = acosd(g[0][1]/(a*b))
173    np.seterr(**oldset)
174    return a,b,c,alp,bet,gam
175
176def invcell2Gmat(invcell):
177    """Compute real and reciprocal lattice metric tensor from reciprocal
178       unit cell constants
179       
180    :param invcell: [a*,b*,c*,alpha*, beta*, gamma*] (degrees)
181    :return: reciprocal (G) & real (g) metric tensors (list of two 3x3 arrays)
182
183    """
184    G = fillgmat(invcell)
185    g = nl.inv(G)
186    return G,g
187       
188def calc_rVsq(A):
189    """Compute the square of the reciprocal lattice volume (1/V**2) from A'
190
191    """
192    G,g = A2Gmat(A)
193    rVsq = nl.det(G)
194    if rVsq < 0:
195        return 1
196    return rVsq
197   
198def calc_rV(A):
199    """Compute the reciprocal lattice volume (V*) from A
200    """
201    return np.sqrt(calc_rVsq(A))
202   
203def calc_V(A):
204    """Compute the real lattice volume (V) from A
205    """
206    return 1./calc_rV(A)
207
208def A2invcell(A):
209    """Compute reciprocal unit cell constants from A
210    returns tuple with a*,b*,c*,alpha*, beta*, gamma* (degrees)
211    """
212    G,g = A2Gmat(A)
213    return Gmat2cell(G)
214   
215def Gmat2AB(G):
216    """Computes orthogonalization matrix from reciprocal metric tensor G
217
218    :returns: tuple of two 3x3 numpy arrays (A,B)
219
220       * A for crystal to Cartesian transformations A*x = np.inner(A,x) = X
221       * B (= inverse of A) for Cartesian to crystal transformation B*X = np.inner(B,X) = x
222
223    """
224    cellstar = Gmat2cell(G)
225    g = nl.inv(G)
226    cell = Gmat2cell(g)
227    A = np.zeros(shape=(3,3))
228    # from Giacovazzo (Fundamentals 2nd Ed.) p.75
229    A[0][0] = cell[0]                # a
230    A[0][1] = cell[1]*cosd(cell[5])  # b cos(gamma)
231    A[0][2] = cell[2]*cosd(cell[4])  # c cos(beta)
232    A[1][1] = cell[1]*sind(cell[5])  # b sin(gamma)
233    A[1][2] = -cell[2]*cosd(cellstar[3])*sind(cell[4]) # - c cos(alpha*) sin(beta)
234    A[2][2] = 1/cellstar[2]         # 1/c*
235    B = nl.inv(A)
236    return A,B
237   
238
239def cell2AB(cell):
240    """Computes orthogonalization matrix from unit cell constants
241
242    :param tuple cell: a,b,c, alpha, beta, gamma (degrees)
243    :returns: tuple of two 3x3 numpy arrays (A,B)
244       A for crystal to Cartesian transformations A*x = np.inner(A,x) = X
245       B (= inverse of A) for Cartesian to crystal transformation B*X = np.inner(B,X) = x
246    """
247    G,g = cell2Gmat(cell) 
248    cellstar = Gmat2cell(G)
249    A = np.zeros(shape=(3,3))
250    # from Giacovazzo (Fundamentals 2nd Ed.) p.75
251    A[0][0] = cell[0]                # a
252    A[0][1] = cell[1]*cosd(cell[5])  # b cos(gamma)
253    A[0][2] = cell[2]*cosd(cell[4])  # c cos(beta)
254    A[1][1] = cell[1]*sind(cell[5])  # b sin(gamma)
255    A[1][2] = -cell[2]*cosd(cellstar[3])*sind(cell[4]) # - c cos(alpha*) sin(beta)
256    A[2][2] = 1/cellstar[2]         # 1/c*
257    B = nl.inv(A)
258    return A,B
259   
260def U6toUij(U6):
261    """Fill matrix (Uij) from U6 = [U11,U22,U33,U12,U13,U23]
262    NB: there is a non numpy version in GSASIIspc: U2Uij
263
264    :param list U6: 6 terms of u11,u22,...
265    :returns:
266        Uij - numpy [3][3] array of uij
267    """
268    U = np.array([
269        [U6[0],  U6[3],  U6[4]], 
270        [U6[3],  U6[1],  U6[5]], 
271        [U6[4],  U6[5],  U6[2]]])
272    return U
273
274def UijtoU6(U):
275    """Fill vector [U11,U22,U33,U12,U13,U23] from Uij
276    NB: there is a non numpy version in GSASIIspc: Uij2U
277    """
278    U6 = np.array([U[0][0],U[1][1],U[2][2],U[0][1],U[0][2],U[1][2]])
279    return U6
280
281def Uij2betaij(Uij,G):
282    """
283    Convert Uij to beta-ij tensors -- stub for eventual completion
284   
285    :param Uij: numpy array [Uij]
286    :param G: reciprocal metric tensor
287    :returns: beta-ij - numpy array [beta-ij]
288    """
289    pass
290   
291def cell2GS(cell):
292    ''' returns Uij to betaij conversion matrix'''
293    G,g = cell2Gmat(cell)
294    GS = G
295    GS[0][1] = GS[1][0] = math.sqrt(GS[0][0]*GS[1][1])
296    GS[0][2] = GS[2][0] = math.sqrt(GS[0][0]*GS[2][2])
297    GS[1][2] = GS[2][1] = math.sqrt(GS[1][1]*GS[2][2])
298    return GS   
299   
300def Uij2Ueqv(Uij,GS,Amat):
301    ''' returns 1/3 trace of diagonalized U matrix'''
302    U = np.multiply(U6toUij(Uij),GS)
303    U = np.inner(Amat,np.inner(U,Amat).T)
304    E,R = nl.eigh(U)
305    return np.sum(E)/3.
306       
307def CosAngle(U,V,G):
308    """ calculate cos of angle between U & V in generalized coordinates
309    defined by metric tensor G
310
311    :param U: 3-vectors assume numpy arrays, can be multiple reflections as (N,3) array
312    :param V: 3-vectors assume numpy arrays, only as (3) vector
313    :param G: metric tensor for U & V defined space assume numpy array
314    :returns:
315        cos(phi)
316    """
317    u = (U.T/np.sqrt(np.sum(np.inner(U,G)*U,axis=1))).T
318    v = V/np.sqrt(np.inner(V,np.inner(G,V)))
319    cosP = np.inner(u,np.inner(G,v))
320    return cosP
321   
322def CosSinAngle(U,V,G):
323    """ calculate sin & cos of angle between U & V in generalized coordinates
324    defined by metric tensor G
325
326    :param U: 3-vectors assume numpy arrays
327    :param V: 3-vectors assume numpy arrays
328    :param G: metric tensor for U & V defined space assume numpy array
329    :returns:
330        cos(phi) & sin(phi)
331    """
332    u = U/np.sqrt(np.inner(U,np.inner(G,U)))
333    v = V/np.sqrt(np.inner(V,np.inner(G,V)))
334    cosP = np.inner(u,np.inner(G,v))
335    sinP = np.sqrt(max(0.0,1.0-cosP**2))
336    return cosP,sinP
337   
338def criticalEllipse(prob):
339    """
340    Calculate critical values for probability ellipsoids from probability
341    """
342    if not ( 0.01 <= prob < 1.0):
343        return 1.54 
344    coeff = np.array([6.44988E-09,4.16479E-07,1.11172E-05,1.58767E-04,0.00130554,
345        0.00604091,0.0114921,-0.040301,-0.6337203,1.311582])
346    llpr = math.log(-math.log(prob))
347    return np.polyval(coeff,llpr)
348   
349def CellBlock(nCells):
350    """
351    Generate block of unit cells n*n*n on a side; [0,0,0] centered, n = 2*nCells+1
352    currently only works for nCells = 0 or 1 (not >1)
353    """
354    if nCells:
355        N = 2*nCells+1
356        N2 = N*N
357        N3 = N*N*N
358        cellArray = []
359        A = np.array(range(N3))
360        cellGen = np.array([A/N2-1,A/N%N-1,A%N-1]).T
361        for cell in cellGen:
362            cellArray.append(cell)
363        return cellArray
364    else:
365        return [0,0,0]
366       
367def CellAbsorption(ElList,Volume):
368    '''Compute unit cell absorption
369
370    :param dict ElList: dictionary of element contents including mu and
371      number of atoms be cell
372    :param float Volume: unit cell volume
373    :returns: mu-total/Volume
374    '''
375    muT = 0
376    for El in ElList:
377        muT += ElList[El]['mu']*ElList[El]['FormulaNo']
378    return muT/Volume
379   
380#Permutations and Combinations
381# Four routines: combinations,uniqueCombinations, selections & permutations
382#These taken from Python Cookbook, 2nd Edition. 19.15 p724-726
383#   
384def _combinators(_handle, items, n):
385    """ factored-out common structure of all following combinators """
386    if n==0:
387        yield [ ]
388        return
389    for i, item in enumerate(items):
390        this_one = [ item ]
391        for cc in _combinators(_handle, _handle(items, i), n-1):
392            yield this_one + cc
393def combinations(items, n):
394    """ take n distinct items, order matters """
395    def skipIthItem(items, i):
396        return items[:i] + items[i+1:]
397    return _combinators(skipIthItem, items, n)
398def uniqueCombinations(items, n):
399    """ take n distinct items, order is irrelevant """
400    def afterIthItem(items, i):
401        return items[i+1:]
402    return _combinators(afterIthItem, items, n)
403def selections(items, n):
404    """ take n (not necessarily distinct) items, order matters """
405    def keepAllItems(items, i):
406        return items
407    return _combinators(keepAllItems, items, n)
408def permutations(items):
409    """ take all items, order matters """
410    return combinations(items, len(items))
411
412#reflection generation routines
413#for these: H = [h,k,l]; A is as used in calc_rDsq; G - inv metric tensor, g - metric tensor;
414#           cell - a,b,c,alp,bet,gam in A & deg
415   
416def Pos2dsp(Inst,pos):
417    ''' convert powder pattern position (2-theta or TOF, musec) to d-spacing
418    ignores secondary effects (e.g. difA in TOF)
419    '''
420    if 'C' in Inst['Type'][0]:
421        wave = G2mth.getWave(Inst)
422        dsp = wave/(2.0*sind((pos-Inst['Zero'][1])/2.0))
423    else:   #'T'OF - ignore difA
424        dsp = (pos-Inst['Zero'][1])/difC
425    return dsp
426   
427def Dsp2pos(Inst,dsp):
428    ''' convert d-spacing to powder pattern position (2-theta or TOF, musec)
429    ignores secondary effects (e.g. difA in TOF) - maybe later?
430    '''
431    if 'C' in Inst['Type'][0]:
432        wave = G2mth.getWave(Inst)
433        pos = 2.0*asind(wave/(2.*dsp))+Inst['Zero'][1]             
434    else:   #'T'OF - ignore difA
435        pos = Inst['difC'][1]*dsp+Inst['Zero'][1]
436    return pos             
437   
438def calc_rDsq(H,A):
439    'needs doc string'
440    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]
441    return rdsq
442   
443def calc_rDsq2(H,G):
444    'needs doc string'
445    return np.inner(H,np.inner(G,H))
446   
447def calc_rDsqZ(H,A,Z,tth,lam):
448    'needs doc string'
449    rpd = np.pi/180.
450    rdsq = calc_rDsq(H,A)+Z*sind(tth)*2.0*rpd/lam**2
451    return rdsq
452       
453def MaxIndex(dmin,A):
454    'needs doc string'
455    Hmax = [0,0,0]
456    try:
457        cell = A2cell(A)
458    except:
459        cell = [1,1,1,90,90,90]
460    for i in range(3):
461        Hmax[i] = int(round(cell[i]/dmin))
462    return Hmax
463   
464def sortHKLd(HKLd,ifreverse,ifdup):
465    '''needs doc string
466
467    :param HKLd: a list of [h,k,l,d,...];
468    :param ifreverse: True for largest d first
469    :param ifdup: True if duplicate d-spacings allowed
470    '''
471    T = []
472    for i,H in enumerate(HKLd):
473        if ifdup:
474            T.append((H[3],i))
475        else:
476            T.append(H[3])           
477    D = dict(zip(T,HKLd))
478    T.sort()
479    if ifreverse:
480        T.reverse()
481    X = []
482    okey = ''
483    for key in T: 
484        if key != okey: X.append(D[key])    #remove duplicate d-spacings
485        okey = key
486    return X
487   
488def SwapIndx(Axis,H):
489    'needs doc string'
490    if Axis in [1,-1]:
491        return H
492    elif Axis in [2,-3]:
493        return [H[1],H[2],H[0]]
494    else:
495        return [H[2],H[0],H[1]]
496       
497def Rh2Hx(Rh):
498    'needs doc string'
499    Hx = [0,0,0]
500    Hx[0] = Rh[0]-Rh[1]
501    Hx[1] = Rh[1]-Rh[2]
502    Hx[2] = np.sum(Rh)
503    return Hx
504   
505def Hx2Rh(Hx):
506    'needs doc string'
507    Rh = [0,0,0]
508    itk = -Hx[0]+Hx[1]+Hx[2]
509    if itk%3 != 0:
510        return 0        #error - not rhombohedral reflection
511    else:
512        Rh[1] = itk/3
513        Rh[0] = Rh[1]+Hx[0]
514        Rh[2] = Rh[1]-Hx[1]
515        if Rh[0] < 0:
516            for i in range(3):
517                Rh[i] = -Rh[i]
518        return Rh
519       
520def CentCheck(Cent,H):
521    'needs doc string'
522    h,k,l = H
523    if Cent == 'A' and (k+l)%2:
524        return False
525    elif Cent == 'B' and (h+l)%2:
526        return False
527    elif Cent == 'C' and (h+k)%2:
528        return False
529    elif Cent == 'I' and (h+k+l)%2:
530        return False
531    elif Cent == 'F' and ((h+k)%2 or (h+l)%2 or (k+l)%2):
532        return False
533    elif Cent == 'R' and (-h+k+l)%3:
534        return False
535    else:
536        return True
537                                   
538def GetBraviasNum(center,system):
539    """Determine the Bravais lattice number, as used in GenHBravais
540   
541    :param center: one of: 'P', 'C', 'I', 'F', 'R' (see SGLatt from GSASIIspc.SpcGroup)
542    :param system: one of 'cubic', 'hexagonal', 'tetragonal', 'orthorhombic', 'trigonal' (for R)
543      'monoclinic', 'triclinic' (see SGSys from GSASIIspc.SpcGroup)
544    :return: a number between 0 and 13
545      or throws a ValueError exception if the combination of center, system is not found (i.e. non-standard)
546
547    """
548    if center.upper() == 'F' and system.lower() == 'cubic':
549        return 0
550    elif center.upper() == 'I' and system.lower() == 'cubic':
551        return 1
552    elif center.upper() == 'P' and system.lower() == 'cubic':
553        return 2
554    elif center.upper() == 'R' and system.lower() == 'trigonal':
555        return 3
556    elif center.upper() == 'P' and system.lower() == 'hexagonal':
557        return 4
558    elif center.upper() == 'I' and system.lower() == 'tetragonal':
559        return 5
560    elif center.upper() == 'P' and system.lower() == 'tetragonal':
561        return 6
562    elif center.upper() == 'F' and system.lower() == 'orthorhombic':
563        return 7
564    elif center.upper() == 'I' and system.lower() == 'orthorhombic':
565        return 8
566    elif center.upper() == 'C' and system.lower() == 'orthorhombic':
567        return 9
568    elif center.upper() == 'P' and system.lower() == 'orthorhombic':
569        return 10
570    elif center.upper() == 'C' and system.lower() == 'monoclinic':
571        return 11
572    elif center.upper() == 'P' and system.lower() == 'monoclinic':
573        return 12
574    elif center.upper() == 'P' and system.lower() == 'triclinic':
575        return 13
576    raise ValueError,'non-standard Bravais lattice center=%s, cell=%s' % (center,system)
577
578def GenHBravais(dmin,Bravais,A):
579    """Generate the positionally unique powder diffraction reflections
580     
581    :param dmin: minimum d-spacing in A
582    :param Bravais: lattice type (see GetBraviasNum). Bravais is one of::
583             0 F cubic
584             1 I cubic
585             2 P cubic
586             3 R hexagonal (trigonal not rhombohedral)
587             4 P hexagonal
588             5 I tetragonal
589             6 P tetragonal
590             7 F orthorhombic
591             8 I orthorhombic
592             9 C orthorhombic
593             10 P orthorhombic
594             11 C monoclinic
595             12 P monoclinic
596             13 P triclinic
597           
598    :param A: reciprocal metric tensor elements as [G11,G22,G33,2*G12,2*G13,2*G23]
599    :return: HKL unique d list of [h,k,l,d,-1] sorted with largest d first
600           
601    """
602    import math
603    if Bravais in [9,11]:
604        Cent = 'C'
605    elif Bravais in [1,5,8]:
606        Cent = 'I'
607    elif Bravais in [0,7]:
608        Cent = 'F'
609    elif Bravais in [3]:
610        Cent = 'R'
611    else:
612        Cent = 'P'
613    Hmax = MaxIndex(dmin,A)
614    dminsq = 1./(dmin**2)
615    HKL = []
616    if Bravais == 13:                       #triclinic
617        for l in range(-Hmax[2],Hmax[2]+1):
618            for k in range(-Hmax[1],Hmax[1]+1):
619                hmin = 0
620                if (k < 0): hmin = 1
621                if (k ==0 and l < 0): hmin = 1
622                for h in range(hmin,Hmax[0]+1):
623                    H=[h,k,l]
624                    rdsq = calc_rDsq(H,A)
625                    if 0 < rdsq <= dminsq:
626                        HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
627    elif Bravais in [11,12]:                #monoclinic - b unique
628        Hmax = SwapIndx(2,Hmax)
629        for h in range(Hmax[0]+1):
630            for k in range(-Hmax[1],Hmax[1]+1):
631                lmin = 0
632                if k < 0:lmin = 1
633                for l in range(lmin,Hmax[2]+1):
634                    [h,k,l] = SwapIndx(-2,[h,k,l])
635                    H = []
636                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
637                    if H:
638                        rdsq = calc_rDsq(H,A)
639                        if 0 < rdsq <= dminsq:
640                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
641                    [h,k,l] = SwapIndx(2,[h,k,l])
642    elif Bravais in [7,8,9,10]:            #orthorhombic
643        for h in range(Hmax[0]+1):
644            for k in range(Hmax[1]+1):
645                for l in range(Hmax[2]+1):
646                    H = []
647                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
648                    if H:
649                        rdsq = calc_rDsq(H,A)
650                        if 0 < rdsq <= dminsq:
651                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
652    elif Bravais in [5,6]:                  #tetragonal
653        for l in range(Hmax[2]+1):
654            for k in range(Hmax[1]+1):
655                for h in range(k,Hmax[0]+1):
656                    H = []
657                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
658                    if H:
659                        rdsq = calc_rDsq(H,A)
660                        if 0 < rdsq <= dminsq:
661                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
662    elif Bravais in [3,4]:
663        lmin = 0
664        if Bravais == 3: lmin = -Hmax[2]                  #hexagonal/trigonal
665        for l in range(lmin,Hmax[2]+1):
666            for k in range(Hmax[1]+1):
667                hmin = k
668                if l < 0: hmin += 1
669                for h in range(hmin,Hmax[0]+1):
670                    H = []
671                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
672                    if H:
673                        rdsq = calc_rDsq(H,A)
674                        if 0 < rdsq <= dminsq:
675                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
676
677    else:                                   #cubic
678        for l in range(Hmax[2]+1):
679            for k in range(l,Hmax[1]+1):
680                for h in range(k,Hmax[0]+1):
681                    H = []
682                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
683                    if H:
684                        rdsq = calc_rDsq(H,A)
685                        if 0 < rdsq <= dminsq:
686                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
687    return sortHKLd(HKL,True,False)
688   
689def getHKLmax(dmin,SGData,A):
690    'finds maximum allowed hkl for given A within dmin'
691    SGLaue = SGData['SGLaue']
692    if SGLaue in ['3R','3mR']:        #Rhombohedral axes
693        Hmax = [0,0,0]
694        cell = A2cell(A)
695        aHx = cell[0]*math.sqrt(2.0*(1.0-cosd(cell[3])))
696        cHx = cell[0]*math.sqrt(3.0*(1.0+2.0*cosd(cell[3])))
697        Hmax[0] = Hmax[1] = int(round(aHx/dmin))
698        Hmax[2] = int(round(cHx/dmin))
699        #print Hmax,aHx,cHx
700    else:                           # all others
701        Hmax = MaxIndex(dmin,A)
702    return Hmax
703   
704def GenHLaue(dmin,SGData,A):
705    """Generate the crystallographically unique powder diffraction reflections
706    for a lattice and Bravais type
707   
708    :param dmin: minimum d-spacing
709    :param SGData: space group dictionary with at least
710   
711        * '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'
712        * 'SGLatt': lattice centering: one of 'P','A','B','C','I','F'
713        * 'SGUniq': code for unique monoclinic axis one of 'a','b','c' (only if 'SGLaue' is '2/m') otherwise an empty string
714       
715    :param A: reciprocal metric tensor elements as [G11,G22,G33,2*G12,2*G13,2*G23]
716    :return: HKL = list of [h,k,l,d] sorted with largest d first and is unique
717            part of reciprocal space ignoring anomalous dispersion
718           
719    """
720    import math
721    SGLaue = SGData['SGLaue']
722    SGLatt = SGData['SGLatt']
723    SGUniq = SGData['SGUniq']
724    #finds maximum allowed hkl for given A within dmin
725    Hmax = getHKLmax(dmin,SGData,A)
726       
727    dminsq = 1./(dmin**2)
728    HKL = []
729    if SGLaue == '-1':                       #triclinic
730        for l in range(-Hmax[2],Hmax[2]+1):
731            for k in range(-Hmax[1],Hmax[1]+1):
732                hmin = 0
733                if (k < 0) or (k ==0 and l < 0): hmin = 1
734                for h in range(hmin,Hmax[0]+1):
735                    H = []
736                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
737                    if H:
738                        rdsq = calc_rDsq(H,A)
739                        if 0 < rdsq <= dminsq:
740                            HKL.append([h,k,l,1/math.sqrt(rdsq)])
741    elif SGLaue == '2/m':                #monoclinic
742        axisnum = 1 + ['a','b','c'].index(SGUniq)
743        Hmax = SwapIndx(axisnum,Hmax)
744        for h in range(Hmax[0]+1):
745            for k in range(-Hmax[1],Hmax[1]+1):
746                lmin = 0
747                if k < 0:lmin = 1
748                for l in range(lmin,Hmax[2]+1):
749                    [h,k,l] = SwapIndx(-axisnum,[h,k,l])
750                    H = []
751                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
752                    if H:
753                        rdsq = calc_rDsq(H,A)
754                        if 0 < rdsq <= dminsq:
755                            HKL.append([h,k,l,1/math.sqrt(rdsq)])
756                    [h,k,l] = SwapIndx(axisnum,[h,k,l])
757    elif SGLaue in ['mmm','4/m','6/m']:            #orthorhombic
758        for l in range(Hmax[2]+1):
759            for h in range(Hmax[0]+1):
760                kmin = 1
761                if SGLaue == 'mmm' or h ==0: kmin = 0
762                for k in range(kmin,Hmax[1]+1):
763                    H = []
764                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
765                    if H:
766                        rdsq = calc_rDsq(H,A)
767                        if 0 < rdsq <= dminsq:
768                            HKL.append([h,k,l,1/math.sqrt(rdsq)])
769    elif SGLaue in ['4/mmm','6/mmm']:                  #tetragonal & hexagonal
770        for l in range(Hmax[2]+1):
771            for h in range(Hmax[0]+1):
772                for k in range(h+1):
773                    H = []
774                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
775                    if H:
776                        rdsq = calc_rDsq(H,A)
777                        if 0 < rdsq <= dminsq:
778                            HKL.append([h,k,l,1/math.sqrt(rdsq)])
779    elif SGLaue in ['3m1','31m','3','3R','3mR']:                  #trigonals
780        for l in range(-Hmax[2],Hmax[2]+1):
781            hmin = 0
782            if l < 0: hmin = 1
783            for h in range(hmin,Hmax[0]+1):
784                if SGLaue in ['3R','3']:
785                    kmax = h
786                    kmin = -int((h-1.)/2.)
787                else:
788                    kmin = 0
789                    kmax = h
790                    if SGLaue in ['3m1','3mR'] and l < 0: kmax = h-1
791                    if SGLaue == '31m' and l < 0: kmin = 1
792                for k in range(kmin,kmax+1):
793                    H = []
794                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
795                    if SGLaue in ['3R','3mR']:
796                        H = Hx2Rh(H)
797                    if H:
798                        rdsq = calc_rDsq(H,A)
799                        if 0 < rdsq <= dminsq:
800                            HKL.append([H[0],H[1],H[2],1/math.sqrt(rdsq)])
801    else:                                   #cubic
802        for h in range(Hmax[0]+1):
803            for k in range(h+1):
804                lmin = 0
805                lmax = k
806                if SGLaue =='m3':
807                    lmax = h-1
808                    if h == k: lmax += 1
809                for l in range(lmin,lmax+1):
810                    H = []
811                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
812                    if H:
813                        rdsq = calc_rDsq(H,A)
814                        if 0 < rdsq <= dminsq:
815                            HKL.append([h,k,l,1/math.sqrt(rdsq)])
816    return sortHKLd(HKL,True,True)
817
818#Spherical harmonics routines
819def OdfChk(SGLaue,L,M):
820    'needs doc string'
821    if not L%2 and abs(M) <= L:
822        if SGLaue == '0':                      #cylindrical symmetry
823            if M == 0: return True
824        elif SGLaue == '-1':
825            return True
826        elif SGLaue == '2/m':
827            if not abs(M)%2: return True
828        elif SGLaue == 'mmm':
829            if not abs(M)%2 and M >= 0: return True
830        elif SGLaue == '4/m':
831            if not abs(M)%4: return True
832        elif SGLaue == '4/mmm':
833            if not abs(M)%4 and M >= 0: return True
834        elif SGLaue in ['3R','3']:
835            if not abs(M)%3: return True
836        elif SGLaue in ['3mR','3m1','31m']:
837            if not abs(M)%3 and M >= 0: return True
838        elif SGLaue == '6/m':
839            if not abs(M)%6: return True
840        elif SGLaue == '6/mmm':
841            if not abs(M)%6 and M >= 0: return True
842        elif SGLaue == 'm3':
843            if M > 0:
844                if L%12 == 2:
845                    if M <= L/12: return True
846                else:
847                    if M <= L/12+1: return True
848        elif SGLaue == 'm3m':
849            if M > 0:
850                if L%12 == 2:
851                    if M <= L/12: return True
852                else:
853                    if M <= L/12+1: return True
854    return False
855       
856def GenSHCoeff(SGLaue,SamSym,L,IfLMN=True):
857    'needs doc string'
858    coeffNames = []
859    for iord in [2*i+2 for i in range(L/2)]:
860        for m in [i-iord for i in range(2*iord+1)]:
861            if OdfChk(SamSym,iord,m):
862                for n in [i-iord for i in range(2*iord+1)]:
863                    if OdfChk(SGLaue,iord,n):
864                        if IfLMN:
865                            coeffNames.append('C(%d,%d,%d)'%(iord,m,n))
866                        else:
867                            coeffNames.append('C(%d,%d)'%(iord,n))
868    return coeffNames
869   
870def CrsAng(H,cell,SGData):
871    'needs doc string'
872    a,b,c,al,be,ga = cell
873    SQ3 = 1.732050807569
874    H1 = np.array([1,0,0])
875    H2 = np.array([0,1,0])
876    H3 = np.array([0,0,1])
877    H4 = np.array([1,1,1])
878    G,g = cell2Gmat(cell)
879    Laue = SGData['SGLaue']
880    Naxis = SGData['SGUniq']
881    DH = np.inner(H,np.inner(G,H))
882    if Laue == '2/m':
883        if Naxis == 'a':
884            DR = np.inner(H1,np.inner(G,H1))
885            DHR = np.inner(H,np.inner(G,H1))
886        elif Naxis == 'b':
887            DR = np.inner(H2,np.inner(G,H2))
888            DHR = np.inner(H,np.inner(G,H2))
889        else:
890            DR = np.inner(H3,np.inner(G,H3))
891            DHR = np.inner(H,np.inner(G,H3))
892    elif Laue in ['R3','R3m']:
893        DR = np.inner(H4,np.inner(G,H4))
894        DHR = np.inner(H,np.inner(G,H4))
895       
896    else:
897        DR = np.inner(H3,np.inner(G,H3))
898        DHR = np.inner(H,np.inner(G,H3))
899    DHR /= np.sqrt(DR*DH)
900    phi = np.where(DHR <= 1.0,acosd(DHR),0.0)
901    if Laue == '-1':
902        BA = H[1]*a/(b-H[0]*cosd(ga))
903        BB = H[0]*sind(ga)**2
904    elif Laue == '2/m':
905        if Naxis == 'a':
906            BA = H[2]*b/(c-H[1]*cosd(al))
907            BB = H[1]*sind(al)**2
908        elif Naxis == 'b':
909            BA = H[0]*c/(a-H[2]*cosd(be))
910            BB = H[2]*sind(be)**2
911        else:
912            BA = H[1]*a/(b-H[0]*cosd(ga))
913            BB = H[0]*sind(ga)**2
914    elif Laue in ['mmm','4/m','4/mmm']:
915        BA = H[1]*a
916        BB = H[0]*b
917   
918    elif Laue in ['3R','3mR']:
919        BA = H[0]+H[1]-2.0*H[2]
920        BB = SQ3*(H[0]-H[1])
921    elif Laue in ['m3','m3m']:
922        BA = H[1]
923        BB = H[0]
924    else:
925        BA = H[0]+2.0*H[1]
926        BB = SQ3*H[0]
927    beta = atan2d(BA,BB)
928    return phi,beta
929   
930def SamAng(Tth,Gangls,Sangl,IFCoup):
931    """Compute sample orientation angles vs laboratory coord. system
932
933    :param Tth:        Signed theta                                   
934    :param Gangls:     Sample goniometer angles phi,chi,omega,azmuth 
935    :param Sangl:      Sample angle zeros om-0, chi-0, phi-0         
936    :param IFCoup:     True if omega & 2-theta coupled in CW scan
937    :returns: 
938        psi,gam:    Sample odf angles                             
939        dPSdA,dGMdA:    Angle zero derivatives
940    """                         
941   
942    rpd = math.pi/180.
943    if IFCoup:
944        GSomeg = sind(Gangls[2]+Tth)
945        GComeg = cosd(Gangls[2]+Tth)
946    else:
947        GSomeg = sind(Gangls[2])
948        GComeg = cosd(Gangls[2])
949    GSTth = sind(Tth)
950    GCTth = cosd(Tth)     
951    GSazm = sind(Gangls[3])
952    GCazm = cosd(Gangls[3])
953    GSchi = sind(Gangls[1])
954    GCchi = cosd(Gangls[1])
955    GSphi = sind(Gangls[0]+Sangl[2])
956    GCphi = cosd(Gangls[0]+Sangl[2])
957    SSomeg = sind(Sangl[0])
958    SComeg = cosd(Sangl[0])
959    SSchi = sind(Sangl[1])
960    SCchi = cosd(Sangl[1])
961    AT = -GSTth*GComeg+GCTth*GCazm*GSomeg
962    BT = GSTth*GSomeg+GCTth*GCazm*GComeg
963    CT = -GCTth*GSazm*GSchi
964    DT = -GCTth*GSazm*GCchi
965   
966    BC1 = -AT*GSphi+(CT+BT*GCchi)*GCphi
967    BC2 = DT-BT*GSchi
968    BC3 = AT*GCphi+(CT+BT*GCchi)*GSphi
969     
970    BC = BC1*SComeg*SCchi+BC2*SComeg*SSchi-BC3*SSomeg     
971    psi = acosd(BC)
972   
973    BD = 1.0-BC**2
974    if BD > 0.:
975        C = rpd/math.sqrt(BD)
976    else:
977        C = 0.
978    dPSdA = [-C*(-BC1*SSomeg*SCchi-BC2*SSomeg*SSchi-BC3*SComeg),
979        -C*(-BC1*SComeg*SSchi+BC2*SComeg*SCchi),
980        -C*(-BC1*SSomeg-BC3*SComeg*SCchi)]
981     
982    BA = -BC1*SSchi+BC2*SCchi
983    BB = BC1*SSomeg*SCchi+BC2*SSomeg*SSchi+BC3*SComeg
984    gam = atan2d(BB,BA)
985
986    BD = (BA**2+BB**2)/rpd
987
988    dBAdO = 0
989    dBAdC = -BC1*SCchi-BC2*SSchi
990    dBAdF = BC3*SSchi
991   
992    dBBdO = BC1*SComeg*SCchi+BC2*SComeg*SSchi-BC3*SSomeg
993    dBBdC = -BC1*SSomeg*SSchi+BC2*SSomeg*SCchi
994    dBBdF = BC1*SComeg-BC3*SSomeg*SCchi
995   
996    if BD > 0.:
997        dGMdA = [(BA*dBBdO-BB*dBAdO)/BD,(BA*dBBdC-BB*dBAdC)/BD,(BA*dBBdF-BB*dBAdF)/BD]
998    else:
999        dGMdA = [0.0,0.0,0.0]
1000
1001       
1002    return psi,gam,dPSdA,dGMdA
1003
1004BOH = {
1005'L=2':[[],[],[]],
1006'L=4':[[0.30469720,0.36418281],[],[]],
1007'L=6':[[-0.14104740,0.52775103],[],[]],
1008'L=8':[[0.28646862,0.21545346,0.32826995],[],[]],
1009'L=10':[[-0.16413497,0.33078546,0.39371345],[],[]],
1010'L=12':[[0.26141975,0.27266871,0.03277460,0.32589402],
1011    [0.09298802,-0.23773812,0.49446631,0.0],[]],
1012'L=14':[[-0.17557309,0.25821932,0.27709173,0.33645360],[],[]],
1013'L=16':[[0.24370673,0.29873515,0.06447688,0.00377,0.32574495],
1014    [0.12039646,-0.25330128,0.23950998,0.40962508,0.0],[]],
1015'L=18':[[-0.16914245,0.17017340,0.34598142,0.07433932,0.32696037],
1016    [-0.06901768,0.16006562,-0.24743528,0.47110273,0.0],[]],
1017'L=20':[[0.23067026,0.31151832,0.09287682,0.01089683,0.00037564,0.32573563],
1018    [0.13615420,-0.25048007,0.12882081,0.28642879,0.34620433,0.0],[]],
1019'L=22':[[-0.16109560,0.10244188,0.36285175,0.13377513,0.01314399,0.32585583],
1020    [-0.09620055,0.20244115,-0.22389483,0.17928946,0.42017231,0.0],[]],
1021'L=24':[[0.22050742,0.31770654,0.11661736,0.02049853,0.00150861,0.00003426,0.32573505],
1022    [0.13651722,-0.21386648,0.00522051,0.33939435,0.10837396,0.32914497,0.0],
1023    [0.05378596,-0.11945819,0.16272298,-0.26449730,0.44923956,0.0,0.0]],
1024'L=26':[[-0.15435003,0.05261630,0.35524646,0.18578869,0.03259103,0.00186197,0.32574594],
1025    [-0.11306511,0.22072681,-0.18706142,0.05439948,0.28122966,0.35634355,0.0],[]],
1026'L=28':[[0.21225019,0.32031716,0.13604702,0.03132468,0.00362703,0.00018294,0.00000294,0.32573501],
1027    [0.13219496,-0.17206256,-0.08742608,0.32671661,0.17973107,0.02567515,0.32619598,0.0],
1028    [0.07989184,-0.16735346,0.18839770,-0.20705337,0.12926808,0.42715602,0.0,0.0]],
1029'L=30':[[-0.14878368,0.01524973,0.33628434,0.22632587,0.05790047,0.00609812,0.00022898,0.32573594],
1030    [-0.11721726,0.20915005,-0.11723436,-0.07815329,0.31318947,0.13655742,0.33241385,0.0],
1031    [-0.04297703,0.09317876,-0.11831248,0.17355132,-0.28164031,0.42719361,0.0,0.0]],
1032'L=32':[[0.20533892,0.32087437,0.15187897,0.04249238,0.00670516,0.00054977,0.00002018,0.00000024,0.32573501],
1033    [0.12775091,-0.13523423,-0.14935701,0.28227378,0.23670434,0.05661270,0.00469819,0.32578978,0.0],
1034    [0.09703829,-0.19373733,0.18610682,-0.14407046,0.00220535,0.26897090,0.36633402,0.0,0.0]],
1035'L=34':[[-0.14409234,-0.01343681,0.31248977,0.25557722,0.08571889,0.01351208,0.00095792,0.00002550,0.32573508],
1036    [-0.11527834,0.18472133,-0.04403280,-0.16908618,0.27227021,0.21086614,0.04041752,0.32688152,0.0],
1037    [-0.06773139,0.14120811,-0.15835721,0.18357456,-0.19364673,0.08377174,0.43116318,0.0,0.0]]
1038}
1039
1040Lnorm = lambda L: 4.*np.pi/(2.0*L+1.)
1041
1042def GetKcl(L,N,SGLaue,phi,beta):
1043    'needs doc string'
1044    import pytexture as ptx
1045    RSQ2PI = 0.3989422804014
1046    SQ2 = 1.414213562373
1047    if SGLaue in ['m3','m3m']:
1048        Kcl = 0.0
1049        for j in range(0,L+1,4):
1050            im = j/4+1
1051            pcrs,dum = ptx.pyplmpsi(L,j,1,phi)
1052            Kcl += BOH['L='+str(L)][N-1][im-1]*pcrs*cosd(j*beta)       
1053    else:
1054        pcrs,dum = ptx.pyplmpsi(L,N,1,phi)
1055        pcrs *= RSQ2PI
1056        if N:
1057            pcrs *= SQ2
1058        if SGLaue in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']:
1059            if SGLaue in ['3mR','3m1','31m']: 
1060                if N%6 == 3:
1061                    Kcl = pcrs*sind(N*beta)
1062                else:
1063                    Kcl = pcrs*cosd(N*beta)
1064            else:
1065                Kcl = pcrs*cosd(N*beta)
1066        else:
1067            Kcl = pcrs*(cosd(N*beta)+sind(N*beta))
1068    return Kcl
1069   
1070def GetKsl(L,M,SamSym,psi,gam):
1071    'needs doc string'
1072    import pytexture as ptx
1073    RSQPI = 0.5641895835478
1074    SQ2 = 1.414213562373
1075    psrs,dpdps = ptx.pyplmpsi(L,M,1,psi)
1076    psrs *= RSQPI
1077    dpdps *= RSQPI
1078    if M == 0:
1079        psrs /= SQ2
1080        dpdps /= SQ2
1081    if SamSym in ['mmm',]:
1082        dum = cosd(M*gam)
1083        Ksl = psrs*dum
1084        dKsdp = dpdps*dum
1085        dKsdg = -psrs*M*sind(M*gam)
1086    else:
1087        dum = cosd(M*gam)+sind(M*gam)
1088        Ksl = psrs*dum
1089        dKsdp = dpdps*dum
1090        dKsdg = psrs*M*(-sind(M*gam)+cosd(M*gam))
1091    return Ksl,dKsdp,dKsdg
1092   
1093def GetKclKsl(L,N,SGLaue,psi,phi,beta):
1094    """
1095    This is used for spherical harmonics description of preferred orientation;
1096        cylindrical symmetry only (M=0) and no sample angle derivatives returned
1097    """
1098    import pytexture as ptx
1099    RSQ2PI = 0.3989422804014
1100    SQ2 = 1.414213562373
1101    Ksl,x = ptx.pyplmpsi(L,0,1,psi)
1102    Ksl *= RSQ2PI
1103    if SGLaue in ['m3','m3m']:
1104        Kcl = 0.0
1105        for j in range(0,L+1,4):
1106            im = j/4+1
1107            pcrs,dum = ptx.pyplmpsi(L,j,1,phi)
1108            Kcl += BOH['L='+str(L)][N-1][im-1]*pcrs*cosd(j*beta)       
1109    else:
1110        pcrs,dum = ptx.pyplmpsi(L,N,1,phi)
1111        pcrs *= RSQ2PI
1112        if N:
1113            pcrs *= SQ2
1114        if SGLaue in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']:
1115            if SGLaue in ['3mR','3m1','31m']: 
1116                if N%6 == 3:
1117                    Kcl = pcrs*sind(N*beta)
1118                else:
1119                    Kcl = pcrs*cosd(N*beta)
1120            else:
1121                Kcl = pcrs*cosd(N*beta)
1122        else:
1123            Kcl = pcrs*(cosd(N*beta)+sind(N*beta))
1124    return Kcl*Ksl,Lnorm(L)
1125   
1126def Glnh(Start,SHCoef,psi,gam,SamSym):
1127    'needs doc string'
1128    import pytexture as ptx
1129    RSQPI = 0.5641895835478
1130    SQ2 = 1.414213562373
1131
1132    if Start:
1133        ptx.pyqlmninit()
1134        Start = False
1135    Fln = np.zeros(len(SHCoef))
1136    for i,term in enumerate(SHCoef):
1137        l,m,n = eval(term.strip('C'))
1138        pcrs,dum = ptx.pyplmpsi(l,m,1,psi)
1139        pcrs *= RSQPI
1140        if m == 0:
1141            pcrs /= SQ2
1142        if SamSym in ['mmm',]:
1143            Ksl = pcrs*cosd(m*gam)
1144        else:
1145            Ksl = pcrs*(cosd(m*gam)+sind(m*gam))
1146        Fln[i] = SHCoef[term]*Ksl*Lnorm(l)
1147    ODFln = dict(zip(SHCoef.keys(),list(zip(SHCoef.values(),Fln))))
1148    return ODFln
1149
1150def Flnh(Start,SHCoef,phi,beta,SGData):
1151    'needs doc string'
1152    import pytexture as ptx
1153   
1154    FORPI = 12.5663706143592
1155    RSQPI = 0.5641895835478
1156    SQ2 = 1.414213562373
1157
1158    if Start:
1159        ptx.pyqlmninit()
1160        Start = False
1161    Fln = np.zeros(len(SHCoef))
1162    for i,term in enumerate(SHCoef):
1163        l,m,n = eval(term.strip('C'))
1164        if SGData['SGLaue'] in ['m3','m3m']:
1165            Kcl = 0.0
1166            for j in range(0,l+1,4):
1167                im = j/4+1
1168                pcrs,dum = ptx.pyplmpsi(l,j,1,phi)
1169                Kcl += BOH['L='+str(l)][n-1][im-1]*pcrs*cosd(j*beta)       
1170        else:                #all but cubic
1171            pcrs,dum = ptx.pyplmpsi(l,n,1,phi)
1172            pcrs *= RSQPI
1173            if n == 0:
1174                pcrs /= SQ2
1175            if SGData['SGLaue'] in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']:
1176               if SGData['SGLaue'] in ['3mR','3m1','31m']: 
1177                   if n%6 == 3:
1178                       Kcl = pcrs*sind(n*beta)
1179                   else:
1180                       Kcl = pcrs*cosd(n*beta)
1181               else:
1182                   Kcl = pcrs*cosd(n*beta)
1183            else:
1184                Kcl = pcrs*(cosd(n*beta)+sind(n*beta))
1185        Fln[i] = SHCoef[term]*Kcl*Lnorm(l)
1186    ODFln = dict(zip(SHCoef.keys(),list(zip(SHCoef.values(),Fln))))
1187    return ODFln
1188   
1189def polfcal(ODFln,SamSym,psi,gam):
1190    'needs doc string'
1191    import pytexture as ptx
1192    RSQPI = 0.5641895835478
1193    SQ2 = 1.414213562373
1194    PolVal = np.ones_like(gam)
1195    for term in ODFln:
1196        if abs(ODFln[term][1]) > 1.e-3:
1197            l,m,n = eval(term.strip('C'))
1198            psrs,dum = ptx.pyplmpsi(l,m,len(psi),psi)
1199            if SamSym in ['-1','2/m']:
1200                if m != 0:
1201                    Ksl = RSQPI*psrs*(cosd(m*gam)+sind(m*gam))
1202                else:
1203                    Ksl = RSQPI*psrs/SQ2
1204            else:
1205                if m != 0:
1206                    Ksl = RSQPI*psrs*cosd(m*gam)
1207                else:
1208                    Ksl = RSQPI*psrs/SQ2
1209            PolVal += ODFln[term][1]*Ksl
1210    return PolVal
1211   
1212def invpolfcal(ODFln,SGData,phi,beta):
1213    'needs doc string'
1214    import pytexture as ptx
1215   
1216    FORPI = 12.5663706143592
1217    RSQPI = 0.5641895835478
1218    SQ2 = 1.414213562373
1219
1220    invPolVal = np.ones_like(beta)
1221    for term in ODFln:
1222        if abs(ODFln[term][1]) > 1.e-3:
1223            l,m,n = eval(term.strip('C'))
1224            if SGData['SGLaue'] in ['m3','m3m']:
1225                Kcl = 0.0
1226                for j in range(0,l+1,4):
1227                    im = j/4+1
1228                    pcrs,dum = ptx.pyplmpsi(l,j,len(beta),phi)
1229                    Kcl += BOH['L='+str(l)][n-1][im-1]*pcrs*cosd(j*beta)       
1230            else:                #all but cubic
1231                pcrs,dum = ptx.pyplmpsi(l,n,len(beta),phi)
1232                pcrs *= RSQPI
1233                if n == 0:
1234                    pcrs /= SQ2
1235                if SGData['SGLaue'] in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']:
1236                   if SGData['SGLaue'] in ['3mR','3m1','31m']: 
1237                       if n%6 == 3:
1238                           Kcl = pcrs*sind(n*beta)
1239                       else:
1240                           Kcl = pcrs*cosd(n*beta)
1241                   else:
1242                       Kcl = pcrs*cosd(n*beta)
1243                else:
1244                    Kcl = pcrs*(cosd(n*beta)+sind(n*beta))
1245            invPolVal += ODFln[term][1]*Kcl
1246    return invPolVal
1247   
1248   
1249def textureIndex(SHCoef):
1250    'needs doc string'
1251    Tindx = 1.0
1252    for term in SHCoef:
1253        l = eval(term.strip('C'))[0]
1254        Tindx += SHCoef[term]**2/(2.0*l+1.)
1255    return Tindx
1256   
1257# self-test materials follow.
1258selftestlist = []
1259'''Defines a list of self-tests'''
1260selftestquiet = True
1261def _ReportTest():
1262    'Report name and doc string of current routine when ``selftestquiet`` is False'
1263    if not selftestquiet:
1264        import inspect
1265        caller = inspect.stack()[1][3]
1266        doc = eval(caller).__doc__
1267        if doc is not None:
1268            print('testing '+__file__+' with '+caller+' ('+doc+')')
1269        else:
1270            print('testing '+__file__()+" with "+caller)
1271NeedTestData = True
1272def TestData():
1273    array = np.array
1274    global NeedTestData
1275    NeedTestData = False
1276    global CellTestData
1277    # output from uctbx computed on platform darwin on 2010-05-28
1278    CellTestData = [
1279# cell, g, G, cell*, V, V*
1280  [(4, 4, 4, 90, 90, 90), 
1281   array([[  1.60000000e+01,   9.79717439e-16,   9.79717439e-16],
1282       [  9.79717439e-16,   1.60000000e+01,   9.79717439e-16],
1283       [  9.79717439e-16,   9.79717439e-16,   1.60000000e+01]]), array([[  6.25000000e-02,   3.82702125e-18,   3.82702125e-18],
1284       [  3.82702125e-18,   6.25000000e-02,   3.82702125e-18],
1285       [  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],
1286# cell, g, G, cell*, V, V*
1287  [(4.0999999999999996, 5.2000000000000002, 6.2999999999999998, 100, 80, 130), 
1288   array([[ 16.81      , -13.70423184,   4.48533243],
1289       [-13.70423184,  27.04      ,  -5.6887143 ],
1290       [  4.48533243,  -5.6887143 ,  39.69      ]]), array([[ 0.10206349,  0.05083339, -0.00424823],
1291       [ 0.05083339,  0.06344997,  0.00334956],
1292       [-0.00424823,  0.00334956,  0.02615544]]), (0.31947376387537696, 0.25189277536327803, 0.16172643497798223, 85.283666420376008, 94.716333579624006, 50.825714168082683), 100.98576357983838, 0.0099023858863968445],
1293# cell, g, G, cell*, V, V*
1294  [(3.5, 3.5, 6, 90, 90, 120), 
1295   array([[  1.22500000e+01,  -6.12500000e+00,   1.28587914e-15],
1296       [ -6.12500000e+00,   1.22500000e+01,   1.28587914e-15],
1297       [  1.28587914e-15,   1.28587914e-15,   3.60000000e+01]]), array([[  1.08843537e-01,   5.44217687e-02,   3.36690552e-18],
1298       [  5.44217687e-02,   1.08843537e-01,   3.36690552e-18],
1299       [  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],
1300  ]
1301    global CoordTestData
1302    CoordTestData = [
1303# cell, ((frac, ortho),...)
1304  ((4,4,4,90,90,90,), [
1305 ((0.10000000000000001, 0.0, 0.0),(0.40000000000000002, 0.0, 0.0)),
1306 ((0.0, 0.10000000000000001, 0.0),(2.4492935982947065e-17, 0.40000000000000002, 0.0)),
1307 ((0.0, 0.0, 0.10000000000000001),(2.4492935982947065e-17, -2.4492935982947065e-17, 0.40000000000000002)),
1308 ((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(0.40000000000000013, 0.79999999999999993, 1.2)),
1309 ((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(0.80000000000000016, 1.2, 0.40000000000000002)),
1310 ((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(1.2, 0.80000000000000004, 0.40000000000000002)),
1311 ((0.5, 0.5, 0.5),(2.0, 1.9999999999999998, 2.0)),
1312]),
1313# cell, ((frac, ortho),...)
1314  ((4.1,5.2,6.3,100,80,130,), [
1315 ((0.10000000000000001, 0.0, 0.0),(0.40999999999999998, 0.0, 0.0)),
1316 ((0.0, 0.10000000000000001, 0.0),(-0.33424955703700043, 0.39834311042186865, 0.0)),
1317 ((0.0, 0.0, 0.10000000000000001),(0.10939835193016617, -0.051013289294572106, 0.6183281045774256)),
1318 ((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(0.069695941716497567, 0.64364635296002093, 1.8549843137322766)),
1319 ((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(-0.073350319180835066, 1.1440160419710339, 0.6183281045774256)),
1320 ((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(0.67089923785616512, 0.74567293154916525, 0.6183281045774256)),
1321 ((0.5, 0.5, 0.5),(0.92574397446582857, 1.7366491056364828, 3.0916405228871278)),
1322]),
1323# cell, ((frac, ortho),...)
1324  ((3.5,3.5,6,90,90,120,), [
1325 ((0.10000000000000001, 0.0, 0.0),(0.35000000000000003, 0.0, 0.0)),
1326 ((0.0, 0.10000000000000001, 0.0),(-0.17499999999999993, 0.3031088913245536, 0.0)),
1327 ((0.0, 0.0, 0.10000000000000001),(3.6739403974420595e-17, -3.6739403974420595e-17, 0.60000000000000009)),
1328 ((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(2.7675166561703527e-16, 0.60621778264910708, 1.7999999999999998)),
1329 ((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(0.17500000000000041, 0.90932667397366063, 0.60000000000000009)),
1330 ((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(0.70000000000000018, 0.6062177826491072, 0.60000000000000009)),
1331 ((0.5, 0.5, 0.5),(0.87500000000000067, 1.5155444566227676, 3.0)),
1332]),
1333]
1334    global LaueTestData             #generated by GSAS
1335    LaueTestData = {
1336    '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),
1337        (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),
1338        (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))],
1339    '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),
1340        (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),
1341        (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),
1342        (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))],
1343    '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),
1344        (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),
1345        (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),
1346        (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),
1347        (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),
1348        (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),
1349        (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),
1350        (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),
1351        (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),
1352        (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),
1353        (2,1,5,6),(2,1,-5,6),(3,-1,-5,6))],
1354    '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),
1355        (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),
1356        (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),
1357        (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),
1358        (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),
1359        (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),
1360        (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),
1361        (3,0,4,6),(3,1,-2,12),(3,0,-4,6),(1,1,6,12),(2,2,3,12))],
1362    '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),
1363        (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),
1364        (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),
1365        (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),
1366        (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),
1367        (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),
1368        (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))],
1369    }
1370   
1371    global FLnhTestData
1372    FLnhTestData = [{
1373    'C(4,0,0)': (0.965, 0.42760447),
1374    'C(2,0,0)': (1.0122, -0.80233610),
1375    'C(2,0,2)': (0.0061, 8.37491546E-03),
1376    'C(6,0,4)': (-0.0898, 4.37985696E-02),
1377    'C(6,0,6)': (-0.1369, -9.04081762E-02),
1378    'C(6,0,0)': (0.5935, -0.18234928),
1379    'C(4,0,4)': (0.1872, 0.16358127),
1380    'C(6,0,2)': (0.6193, 0.27573633),
1381    'C(4,0,2)': (-0.1897, 0.12530720)},[1,0,0]]
1382def test0():
1383    if NeedTestData: TestData()
1384    msg = 'test cell2Gmat, fillgmat, Gmat2cell'
1385    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
1386        G, g = cell2Gmat(cell)
1387        assert np.allclose(G,tG),msg
1388        assert np.allclose(g,tg),msg
1389        tcell = Gmat2cell(g)
1390        assert np.allclose(cell,tcell),msg
1391        tcell = Gmat2cell(G)
1392        assert np.allclose(tcell,trcell),msg
1393selftestlist.append(test0)
1394
1395def test1():
1396    'test cell2A and A2Gmat'
1397    _ReportTest()
1398    if NeedTestData: TestData()
1399    msg = 'test cell2A and A2Gmat'
1400    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
1401        G, g = A2Gmat(cell2A(cell))
1402        assert np.allclose(G,tG),msg
1403        assert np.allclose(g,tg),msg
1404selftestlist.append(test1)
1405
1406def test2():
1407    'test Gmat2A, A2cell, A2Gmat, Gmat2cell'
1408    _ReportTest()
1409    if NeedTestData: TestData()
1410    msg = 'test Gmat2A, A2cell, A2Gmat, Gmat2cell'
1411    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
1412        G, g = cell2Gmat(cell)
1413        tcell = A2cell(Gmat2A(G))
1414        assert np.allclose(cell,tcell),msg
1415selftestlist.append(test2)
1416
1417def test3():
1418    'test invcell2Gmat'
1419    _ReportTest()
1420    if NeedTestData: TestData()
1421    msg = 'test invcell2Gmat'
1422    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
1423        G, g = invcell2Gmat(trcell)
1424        assert np.allclose(G,tG),msg
1425        assert np.allclose(g,tg),msg
1426selftestlist.append(test3)
1427
1428def test4():
1429    'test calc_rVsq, calc_rV, calc_V'
1430    _ReportTest()
1431    if NeedTestData: TestData()
1432    msg = 'test calc_rVsq, calc_rV, calc_V'
1433    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
1434        assert np.allclose(calc_rV(cell2A(cell)),trV), msg
1435        assert np.allclose(calc_V(cell2A(cell)),tV), msg
1436selftestlist.append(test4)
1437
1438def test5():
1439    'test A2invcell'
1440    _ReportTest()
1441    if NeedTestData: TestData()
1442    msg = 'test A2invcell'
1443    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
1444        rcell = A2invcell(cell2A(cell))
1445        assert np.allclose(rcell,trcell),msg
1446selftestlist.append(test5)
1447
1448def test6():
1449    'test cell2AB'
1450    _ReportTest()
1451    if NeedTestData: TestData()
1452    msg = 'test cell2AB'
1453    for (cell,coordlist) in CoordTestData:
1454        A,B = cell2AB(cell)
1455        for (frac,ortho) in coordlist:
1456            to = np.inner(A,frac)
1457            tf = np.inner(B,to)
1458            assert np.allclose(ortho,to), msg
1459            assert np.allclose(frac,tf), msg
1460            to = np.sum(A*frac,axis=1)
1461            tf = np.sum(B*to,axis=1)
1462            assert np.allclose(ortho,to), msg
1463            assert np.allclose(frac,tf), msg
1464selftestlist.append(test6)
1465
1466def test7():
1467    'test GetBraviasNum(...) and GenHBravais(...)'
1468    _ReportTest()
1469    import os.path
1470    import sys
1471    import GSASIIspc as spc
1472    testdir = os.path.join(os.path.split(os.path.abspath( __file__ ))[0],'testinp')
1473    if os.path.exists(testdir):
1474        if testdir not in sys.path: sys.path.insert(0,testdir)
1475    import sgtbxlattinp
1476    derror = 1e-4
1477    def indexmatch(hklin, hkllist, system):
1478        for hklref in hkllist:
1479            hklref = list(hklref)
1480            # these permutations are far from complete, but are sufficient to
1481            # allow the test to complete
1482            if system == 'cubic':
1483                permlist = [(1,2,3),(1,3,2),(2,1,3),(2,3,1),(3,1,2),(3,2,1),]
1484            elif system == 'monoclinic':
1485                permlist = [(1,2,3),(-1,2,-3)]
1486            else:
1487                permlist = [(1,2,3)]
1488
1489            for perm in permlist:
1490                hkl = [abs(i) * hklin[abs(i)-1] / i for i in perm]
1491                if hkl == hklref: return True
1492                if [-i for i in hkl] == hklref: return True
1493        else:
1494            return False
1495
1496    for key in sgtbxlattinp.sgtbx7:
1497        spdict = spc.SpcGroup(key)
1498        cell = sgtbxlattinp.sgtbx7[key][0]
1499        system = spdict[1]['SGSys']
1500        center = spdict[1]['SGLatt']
1501
1502        bravcode = GetBraviasNum(center, system)
1503
1504        g2list = GenHBravais(sgtbxlattinp.dmin, bravcode, cell2A(cell))
1505
1506        assert len(sgtbxlattinp.sgtbx7[key][1]) == len(g2list), 'Reflection lists differ for %s' % key
1507        for h,k,l,d,num in g2list:
1508            for hkllist,dref in sgtbxlattinp.sgtbx7[key][1]: 
1509                if abs(d-dref) < derror:
1510                    if indexmatch((h,k,l,), hkllist, system):
1511                        break
1512            else:
1513                assert 0,'No match for %s at %s (%s)' % ((h,k,l),d,key)
1514selftestlist.append(test7)
1515
1516def test8():
1517    'test GenHLaue'
1518    _ReportTest()
1519    import GSASIIspc as spc
1520    import sgtbxlattinp
1521    derror = 1e-4
1522    dmin = sgtbxlattinp.dmin
1523
1524    def indexmatch(hklin, hklref, system, axis):
1525        # these permutations are far from complete, but are sufficient to
1526        # allow the test to complete
1527        if system == 'cubic':
1528            permlist = [(1,2,3),(1,3,2),(2,1,3),(2,3,1),(3,1,2),(3,2,1),]
1529        elif system == 'monoclinic' and axis=='b':
1530            permlist = [(1,2,3),(-1,2,-3)]
1531        elif system == 'monoclinic' and axis=='a':
1532            permlist = [(1,2,3),(1,-2,-3)]
1533        elif system == 'monoclinic' and axis=='c':
1534            permlist = [(1,2,3),(-1,-2,3)]
1535        elif system == 'trigonal':
1536            permlist = [(1,2,3),(2,1,3),(-1,-2,3),(-2,-1,3)]
1537        elif system == 'rhombohedral':
1538            permlist = [(1,2,3),(2,3,1),(3,1,2)]
1539        else:
1540            permlist = [(1,2,3)]
1541
1542        hklref = list(hklref)
1543        for perm in permlist:
1544            hkl = [abs(i) * hklin[abs(i)-1] / i for i in perm]
1545            if hkl == hklref: return True
1546            if [-i for i in hkl] == hklref: return True
1547        return False
1548
1549    for key in sgtbxlattinp.sgtbx8:
1550        spdict = spc.SpcGroup(key)[1]
1551        cell = sgtbxlattinp.sgtbx8[key][0]
1552        center = spdict['SGLatt']
1553        Laue = spdict['SGLaue']
1554        Axis = spdict['SGUniq']
1555        system = spdict['SGSys']
1556
1557        g2list = GenHLaue(dmin,spdict,cell2A(cell))
1558        #if len(g2list) != len(sgtbxlattinp.sgtbx8[key][1]):
1559        #    print 'failed',key,':' ,len(g2list),'vs',len(sgtbxlattinp.sgtbx8[key][1])
1560        #    print 'GSAS-II:'
1561        #    for h,k,l,d in g2list: print '  ',(h,k,l),d
1562        #    print 'SGTBX:'
1563        #    for hkllist,dref in sgtbxlattinp.sgtbx8[key][1]: print '  ',hkllist,dref
1564        assert len(g2list) == len(sgtbxlattinp.sgtbx8[key][1]), (
1565            'Reflection lists differ for %s' % key
1566            )
1567        #match = True
1568        for h,k,l,d in g2list:
1569            for hkllist,dref in sgtbxlattinp.sgtbx8[key][1]: 
1570                if abs(d-dref) < derror:
1571                    if indexmatch((h,k,l,), hkllist, system, Axis): break
1572            else:
1573                assert 0,'No match for %s at %s (%s)' % ((h,k,l),d,key)
1574                #match = False
1575        #if not match:
1576            #for hkllist,dref in sgtbxlattinp.sgtbx8[key][1]: print '  ',hkllist,dref
1577            #print center, Laue, Axis, system
1578selftestlist.append(test8)
1579           
1580def test9():
1581    'test GenHLaue'
1582    _ReportTest()
1583    import GSASIIspc as G2spc
1584    if NeedTestData: TestData()
1585    for spc in LaueTestData:
1586        data = LaueTestData[spc]
1587        cell = data[0]
1588        hklm = np.array(data[1])
1589        H = hklm[-1][:3]
1590        hklO = hklm.T[:3].T
1591        A = cell2A(cell)
1592        dmin = 1./np.sqrt(calc_rDsq(H,A))
1593        SGData = G2spc.SpcGroup(spc)[1]
1594        hkls = np.array(GenHLaue(dmin,SGData,A))
1595        hklN = hkls.T[:3].T
1596        #print spc,hklO.shape,hklN.shape
1597        err = True
1598        for H in hklO:
1599            if H not in hklN:
1600                print H,' missing from hkl from GSASII'
1601                err = False
1602        assert(err)
1603selftestlist.append(test9)
1604       
1605       
1606   
1607
1608if __name__ == '__main__':
1609    # run self-tests
1610    selftestquiet = False
1611    for test in selftestlist:
1612        test()
1613    print "OK"
Note: See TracBrowser for help on using the repository browser.