source: trunk/GSASIIlattice.py @ 1597

Last change on this file since 1597 was 1597, checked in by vondreele, 7 years ago

Incommensurate Pawley refinement works, including refinement of modulation vector.

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