source: trunk/GSASIIlattice.py @ 2038

Last change on this file since 2038 was 2038, checked in by vondreele, 8 years ago

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