source: trunk/GSASIIlattice.py @ 2107

Last change on this file since 2107 was 2038, checked in by vondreele, 10 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: toby $
28# $Revision: 2038 $
29# $URL: trunk/GSASIIlattice.py $
30# $Id: GSASIIlattice.py 2038 2015-10-30 21:02:02Z toby $
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.