source: trunk/GSASIIlattice.py @ 378

Last change on this file since 378 was 378, checked in by vondreele, 11 years ago

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