source: trunk/GSASIIlattice.py @ 2802

Last change on this file since 2802 was 2802, checked in by toby, 6 years ago

updates for doc build

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