source: trunk/GSASIIlattice.py @ 1489

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

correct Uij to U6 conversions Uij = Uk/2 for i != j; k > 2 in G2lattice & G2spc
This fixes derivative issues for aniso tropic atoms in Al2O3
change names of sytsyms e.g. 2(100) to 2(x) for non-hex/trig space groups. Removes possible ambiguity with ortho, etc. cases.
Still checking extinction issues

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