source: trunk/GSASIIlattice.py @ 2135

Last change on this file since 2135 was 2135, checked in by vondreele, 6 years ago

more fixes to peak exports to get correct FWHM
refactor Merge GUI

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