source: trunk/GSASIIlattice.py @ 326

Last change on this file since 326 was 326, checked in by vondreele, 12 years ago

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