source: trunk/GSASIIlattice.py @ 1598

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

insure super lattice reflection set is unique reflections only
fix PawleyUpDate? for SS reflections
add plot of recip. unit cell box in 3D reflection drawing
fix SS plot math in 3Dplots

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