source: trunk/GSASIIlattice.py @ 1884

Last change on this file since 1884 was 1820, checked in by vondreele, 10 years ago

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