source: trunk/GSASIIlattice.py @ 432

Last change on this file since 432 was 432, checked in by vondreele, 10 years ago

remove test prints from GSASIIIO
new Gmat2AB routine
work on the DData GUI bug - doesn't crash but bad GUI shown
work on ellipsoidal crystallites - all now work.

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