# source:trunk/GSASIIlattice.py@1587

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

use successive approx. to get d from TOF
allow sorting on a,b,c,alp,bet,gam, vol & m20 in cell indexed table
put in derivatives for fitCell inside indexing routine - better accuracy & speed
add a couple more orthos to SS table

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