source: trunk/GSASIIlattice.py @ 2484

Last change on this file since 2484 was 2484, checked in by vondreele, 5 years ago

fix transform to properly move magnetic moments (if space group is not changed)
some work on magnetic derivatives

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