source: trunk/GSASIIlattice.py @ 885

Last change on this file since 885 was 885, checked in by vondreele, 9 years ago

add 'Global' to constraints
remove a stray debug print
fix U6toUij so it returns a numpy array
fixes to Q2AV & Q2AVdeg
RB changes - get derivatives in

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