source: trunk/GSASIIlattice.py @ 2775

Last change on this file since 2775 was 2775, checked in by toby, 5 years ago

make autoint polltime adjustable

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 95.3 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-10 16:33:13 +0000 (Mon, 10 Apr 2017) $
27# $Author: toby $
28# $Revision: 2775 $
29# $URL: trunk/GSASIIlattice.py $
30# $Id: GSASIIlattice.py 2775 2017-04-10 16:33:13Z 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: 2775 $")
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     ['-1','2/m','112/m','2/m11','mmm','-42m','-4m2','4/mmm','-3','-31m','-3m1','6/m','6/mmm','m3','m3m']
1213     
1214      noncentrosymmetric Laue groups::
1215     
1216     ['1','2','211','112','m','m11','11m','222','mm2','m2m','2mm','4','-4','422','4mm','3','312','321','31m','3m1','6','-6','622','6mm','-62m','-6m2','23','432','-43m']
1217     
1218    :param HKLF: np.array([[h,k,l,...]]) reflection set to be converted
1219   
1220    :return: HKLF new reflection array with imposed Laue symmetry
1221    '''
1222   
1223    HKLFT = HKLF.T
1224    mat41 = np.array([[0,1,0],[-1,0,0],[0,0,1]])    #hkl -> k,-h,l
1225    mat43 = np.array([[0,-1,0],[1,0,0],[0,0,1]])    #hkl -> -k,h,l
1226    mat4bar = np.array([[0,-1,0],[1,0,0],[0,0,-1]]) #hkl -> k,-h,-l
1227    mat31 = np.array([[-1,-1,0],[1,0,0],[0,0,1]])   #hkl -> ihl = -h-k,h,l
1228    mat32 = np.array([[0,1,0],[-1,-1,0],[0,0,1]])   #hkl -> kil = k,-h-k,l
1229    matd3 = np.array([[0,1,0],[0,0,1],[1,0,0]])     #hkl -> k,l,h
1230    matd3q = np.array([[0,0,-1],[-1,0,0],[0,1,0]])  #hkl -> -l,-h,k
1231    matd3t = np.array([[0,0,-1],[1,0,0],[0,-1,0]])  #hkl -> -l,h,-k
1232    mat6 = np.array([[1,1,0],[-1,0,0],[0,0,1]])     #hkl -> h+k,-h,l really 65
1233    matdm = np.array([[0,1,0],[1,0,0],[0,0,1]])     #hkl -> k,h,l
1234    matdmp = np.array([[-1,-1,0],[0,1,0],[0,0,1]])  #hkl -> -h-k,k,l
1235    matkm = np.array([[-1,0,0],[1,1,0],[0,0,1]])    #hkl -> -h,h+k,l
1236    matd2 = np.array([[0,1,0],[1,0,0],[0,0,-1]])    #hkl -> k,h,-l
1237    matdm3 = np.array([[1,0,0],[0,0,1],[0,1,0]])    #hkl -> h,l,k
1238    mat2d43 = np.array([[0,1,0],[1,0,0],[0,0,1]])   #hkl -> k,-h,l
1239    matk2 = np.array([[-1,0,0],[1,1,0],[0,0,-1]])   #hkl -> -h,-i,-l
1240    #triclinic
1241    if Laue == '1': #ok
1242        pass
1243    elif Laue == '-1':  #ok
1244        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
1245        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]<0),HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
1246        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
1247    #monoclinic
1248    #noncentrosymmetric - all ok
1249    elif Laue == '2': 
1250        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3])
1251        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3])
1252    elif Laue == '1 1 2':
1253        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1254        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]<0),HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1255    elif Laue == '2 1 1':   
1256        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1257        HKLFT[:3] = np.where((HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1258    elif Laue == 'm':
1259        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1260    elif Laue == 'm 1 1':
1261        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1262    elif Laue == '1 1 m':
1263        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1264    #centrosymmetric - all ok
1265    elif Laue == '2/m 1 1':       
1266        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1267        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1268        HKLFT[:3] = np.where((HKLFT[2]*HKLFT[0]==0)&(HKLFT[1]<0),HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1269    elif Laue == '2/m':
1270        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3])
1271        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1272        HKLFT[:3] = np.where((HKLFT[0]*HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1273    elif Laue == '1 1 2/m':
1274        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1275        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1276        HKLFT[:3] = np.where((HKLFT[1]*HKLFT[2]==0)&(HKLFT[0]<0),HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1277    #orthorhombic
1278    #noncentrosymmetric - all OK
1279    elif Laue == '2 2 2':
1280        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1281        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1282        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1283        HKLFT[:3] = np.where((HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1284    elif Laue == 'm m 2':
1285        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1286        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1287    elif Laue == '2 m m': 
1288        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1289        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1290    elif Laue == 'm 2 m':
1291        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1292        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1293    #centrosymmetric - all ok
1294    elif Laue == 'm m m':
1295        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1296        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1297        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1298    #tetragonal
1299    #noncentrosymmetric - all ok
1300    elif Laue == '4':
1301        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1302        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1303        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]>0),np.squeeze(np.inner(HKLF[:,:3],mat41[nxs,:,:])).T,HKLFT[:3])
1304    elif Laue == '-4': 
1305        HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])     
1306        HKLFT[:3] = np.where(HKLFT[0]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1307        HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])     
1308        HKLFT[:3] = np.where(HKLFT[1]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1309        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1310    elif Laue == '4 2 2':
1311        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1312        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1313        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1314        HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]<HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1315        HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])   #in lieu od 2-fold
1316    elif Laue == '4 m m':
1317        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1318        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1319        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1320        HKLFT[:3] = np.where(HKLFT[0]<HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1321    elif Laue == '-4 2 m':
1322        HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])     
1323        HKLFT[:3] = np.where(HKLFT[0]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1324        HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])     
1325        HKLFT[:3] = np.where(HKLFT[1]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1326        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1327        HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1328        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1329    elif Laue == '-4 m 2':
1330        HKLFT[:3] = np.where(HKLFT[2]<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[2]==0)&(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[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])     
1334        HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]==0),np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1335        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3]) 
1336        HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[0]>HKLFT[1]),np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1337    #centrosymmetric - all ok
1338    elif Laue == '4/m':
1339        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1340        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1341        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1342        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]>0),np.squeeze(np.inner(HKLF[:,:3],mat41[nxs,:,:])).T,HKLFT[:3])
1343    elif Laue == '4/m m m':
1344        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1345        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1346        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])       
1347        HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],mat41[nxs,:,:])).T,HKLFT[:3])
1348        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1349    #trigonal - all hex cell
1350    #noncentrosymmetric - all ok
1351    elif Laue == '3':
1352        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1353        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1354        HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1355    elif Laue == '3 1 2':
1356        HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matk2[nxs,:,:])).T,HKLFT[:3])
1357        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1358        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1359        HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1360        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],matk2[nxs,:,:])).T,HKLFT[:3])
1361    elif Laue == '3 2 1':
1362        HKLFT[:3] = np.where(HKLFT[0]<=-2*HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1363        HKLFT[:3] = np.where(HKLFT[1]<-2*HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1364        HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1365        HKLFT[:3] = np.where((HKLFT[2]>0)&(HKLFT[1]==HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1366        HKLFT[:3] = np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T
1367        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])
1368    elif Laue == '3 1 m':
1369        HKLFT[:3] = np.where(HKLFT[0]>=HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1370        HKLFT[:3] = np.where(2*HKLFT[1]<-HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1371        HKLFT[:3] = np.where(HKLFT[1]>-2*HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matdmp[nxs,:,:])).T,HKLFT[:3])
1372        HKLFT[:3] = np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T
1373    elif Laue == '3 m 1':
1374        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1375        HKLFT[:3] = np.where((HKLFT[1]+HKLFT[0])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1376        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],matkm[nxs,:,:])).T,HKLFT[:3])
1377    #centrosymmetric
1378    elif Laue == '-3':  #ok
1379        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
1380        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1381        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1382        HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1383        HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[0]<0),-np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1384        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],-mat31[nxs,:,:])).T,HKLFT[:3])   
1385    elif Laue == '-3 m 1':  #ok
1386        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1387        HKLFT[:3] = np.where((HKLFT[1]+HKLFT[0])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1388        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],matkm[nxs,:,:])).T,HKLFT[:3])
1389        HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1390        HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]<HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1391    elif Laue == '-3 1 m':  #ok
1392        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
1393        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1394        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1395        HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1396        HKLFT[:3] = np.where(HKLFT[0]<=0,np.squeeze(np.inner(HKLF[:,:3],-mat31[nxs,:,:])).T,HKLFT[:3])   
1397        HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1398    #hexagonal
1399    #noncentrosymmetric
1400    elif Laue == '6':   #ok
1401        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1402        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1403        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1404        HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1405    elif Laue == '-6':  #ok
1406        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1407        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1408        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1409        HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1410    elif Laue == '6 2 2':   #ok
1411        HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1412        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1413        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1414        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1415        HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1416        HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[0]>HKLFT[1]),np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1417    elif Laue == '6 m m':   #ok
1418        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1419        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1420        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[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[0]>HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1423    elif Laue == '-6 m 2':  #ok
1424        HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matk2[nxs,:,:])).T,HKLFT[:3])
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[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1428        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],matk2[nxs,:,:])).T,HKLFT[:3])
1429        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1430    elif Laue == '-6 2 m':  #ok
1431        HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1432        HKLFT[:3] = np.where(HKLFT[0]<=-2*HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1433        HKLFT[:3] = np.where(HKLFT[1]<-2*HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1434        HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1435        HKLFT[:3] = np.where((HKLFT[2]>0)&(HKLFT[1]==HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1436        HKLFT[:3] = np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T
1437        HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1438        HKLFT[:3] = np.where(HKLFT[0]>HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1439    #centrosymmetric
1440    elif Laue == '6/m': #ok
1441        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1442        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1443        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1444        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1445        HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1446    elif Laue == '6/m m m': #ok
1447        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1448        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1449        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1450        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1451        HKLFT[:3] = np.where(HKLFT[0]>HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm.T[nxs,:,:])).T,HKLFT[:3])
1452    #cubic - all ok
1453    #noncentrosymmetric -
1454    elif Laue == '2 3': 
1455        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1456        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1457        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1458        HKLFT[:3] = np.where((HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1459        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])
1460        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])
1461        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])
1462        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])
1463        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3])       
1464    elif Laue == '4 3 2':   
1465        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1466        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1467        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1468        HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]<HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1469        HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])   #in lieu od 2-fold
1470        HKLFT[:3] = np.where((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2]),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
1471        HKLFT[:3] = np.where((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2]),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
1472        HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat2d43[nxs,:,:])).T,HKLFT[:3])
1473    elif Laue == '-4 3 m': 
1474        HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])     
1475        HKLFT[:3] = np.where(HKLFT[0]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1476        HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])     
1477        HKLFT[:3] = np.where(HKLFT[1]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1478        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1479        HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1480        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1481        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])
1482        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])
1483        HKLFT[:3] = np.where((HKLFT[2]>=0)&(HKLFT[1]<HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1484        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3]) 
1485        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])
1486        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])
1487    #centrosymmetric
1488    elif Laue == 'm 3':
1489        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1490        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],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[2]>=0)&((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
1493        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])
1494    elif Laue == 'm 3 m':
1495        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1496        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1497        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])           
1498        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])
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[0]>HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1501    return HKLFT.T
1502       
1503
1504#Spherical harmonics routines
1505def OdfChk(SGLaue,L,M):
1506    'needs doc string'
1507    if not L%2 and abs(M) <= L:
1508        if SGLaue == '0':                      #cylindrical symmetry
1509            if M == 0: return True
1510        elif SGLaue == '-1':
1511            return True
1512        elif SGLaue == '2/m':
1513            if not abs(M)%2: return True
1514        elif SGLaue == 'mmm':
1515            if not abs(M)%2 and M >= 0: return True
1516        elif SGLaue == '4/m':
1517            if not abs(M)%4: return True
1518        elif SGLaue == '4/mmm':
1519            if not abs(M)%4 and M >= 0: return True
1520        elif SGLaue in ['3R','3']:
1521            if not abs(M)%3: return True
1522        elif SGLaue in ['3mR','3m1','31m']:
1523            if not abs(M)%3 and M >= 0: return True
1524        elif SGLaue == '6/m':
1525            if not abs(M)%6: return True
1526        elif SGLaue == '6/mmm':
1527            if not abs(M)%6 and M >= 0: return True
1528        elif SGLaue == 'm3':
1529            if M > 0:
1530                if L%12 == 2:
1531                    if M <= L/12: return True
1532                else:
1533                    if M <= L/12+1: return True
1534        elif SGLaue == 'm3m':
1535            if M > 0:
1536                if L%12 == 2:
1537                    if M <= L/12: return True
1538                else:
1539                    if M <= L/12+1: return True
1540    return False
1541       
1542def GenSHCoeff(SGLaue,SamSym,L,IfLMN=True):
1543    'needs doc string'
1544    coeffNames = []
1545    for iord in [2*i+2 for i in range(L/2)]:
1546        for m in [i-iord for i in range(2*iord+1)]:
1547            if OdfChk(SamSym,iord,m):
1548                for n in [i-iord for i in range(2*iord+1)]:
1549                    if OdfChk(SGLaue,iord,n):
1550                        if IfLMN:
1551                            coeffNames.append('C(%d,%d,%d)'%(iord,m,n))
1552                        else:
1553                            coeffNames.append('C(%d,%d)'%(iord,n))
1554    return coeffNames
1555   
1556def CrsAng(H,cell,SGData):
1557    'needs doc string'
1558    a,b,c,al,be,ga = cell
1559    SQ3 = 1.732050807569
1560    H1 = np.array([1,0,0])
1561    H2 = np.array([0,1,0])
1562    H3 = np.array([0,0,1])
1563    H4 = np.array([1,1,1])
1564    G,g = cell2Gmat(cell)
1565    Laue = SGData['SGLaue']
1566    Naxis = SGData['SGUniq']
1567    if len(H.shape) == 1:
1568        DH = np.inner(H,np.inner(G,H))
1569    else:
1570        DH = np.array([np.inner(h,np.inner(G,h)) for h in H])
1571    if Laue == '2/m':
1572        if Naxis == 'a':
1573            DR = np.inner(H1,np.inner(G,H1))
1574            DHR = np.inner(H,np.inner(G,H1))
1575        elif Naxis == 'b':
1576            DR = np.inner(H2,np.inner(G,H2))
1577            DHR = np.inner(H,np.inner(G,H2))
1578        else:
1579            DR = np.inner(H3,np.inner(G,H3))
1580            DHR = np.inner(H,np.inner(G,H3))
1581    elif Laue in ['R3','R3m']:
1582        DR = np.inner(H4,np.inner(G,H4))
1583        DHR = np.inner(H,np.inner(G,H4))
1584    else:
1585        DR = np.inner(H3,np.inner(G,H3))
1586        DHR = np.inner(H,np.inner(G,H3))
1587    DHR /= np.sqrt(DR*DH)
1588    phi = np.where(DHR <= 1.0,acosd(DHR),0.0)
1589    if Laue == '-1':
1590        BA = H.T[1]*a/(b-H.T[0]*cosd(ga))
1591        BB = H.T[0]*sind(ga)**2
1592    elif Laue == '2/m':
1593        if Naxis == 'a':
1594            BA = H.T[2]*b/(c-H.T[1]*cosd(al))
1595            BB = H.T[1]*sind(al)**2
1596        elif Naxis == 'b':
1597            BA = H.T[0]*c/(a-H.T[2]*cosd(be))
1598            BB = H.T[2]*sind(be)**2
1599        else:
1600            BA = H.T[1]*a/(b-H.T[0]*cosd(ga))
1601            BB = H.T[0]*sind(ga)**2
1602    elif Laue in ['mmm','4/m','4/mmm']:
1603        BA = H.T[1]*a
1604        BB = H.T[0]*b
1605    elif Laue in ['3R','3mR']:
1606        BA = H.T[0]+H.T[1]-2.0*H.T[2]
1607        BB = SQ3*(H.T[0]-H.T[1])
1608    elif Laue in ['m3','m3m']:
1609        BA = H.T[1]
1610        BB = H.T[0]
1611    else:
1612        BA = H.T[0]+2.0*H.T[1]
1613        BB = SQ3*H.T[0]
1614    beta = atan2d(BA,BB)
1615    return phi,beta
1616   
1617def SamAng(Tth,Gangls,Sangl,IFCoup):
1618    """Compute sample orientation angles vs laboratory coord. system
1619
1620    :param Tth:        Signed theta                                   
1621    :param Gangls:     Sample goniometer angles phi,chi,omega,azmuth 
1622    :param Sangl:      Sample angle zeros om-0, chi-0, phi-0         
1623    :param IFCoup:     True if omega & 2-theta coupled in CW scan
1624    :returns: 
1625        psi,gam:    Sample odf angles                             
1626        dPSdA,dGMdA:    Angle zero derivatives
1627    """                         
1628   
1629    if IFCoup:
1630        GSomeg = sind(Gangls[2]+Tth)
1631        GComeg = cosd(Gangls[2]+Tth)
1632    else:
1633        GSomeg = sind(Gangls[2])
1634        GComeg = cosd(Gangls[2])
1635    GSTth = sind(Tth)
1636    GCTth = cosd(Tth)     
1637    GSazm = sind(Gangls[3])
1638    GCazm = cosd(Gangls[3])
1639    GSchi = sind(Gangls[1])
1640    GCchi = cosd(Gangls[1])
1641    GSphi = sind(Gangls[0]+Sangl[2])
1642    GCphi = cosd(Gangls[0]+Sangl[2])
1643    SSomeg = sind(Sangl[0])
1644    SComeg = cosd(Sangl[0])
1645    SSchi = sind(Sangl[1])
1646    SCchi = cosd(Sangl[1])
1647    AT = -GSTth*GComeg+GCTth*GCazm*GSomeg
1648    BT = GSTth*GSomeg+GCTth*GCazm*GComeg
1649    CT = -GCTth*GSazm*GSchi
1650    DT = -GCTth*GSazm*GCchi
1651   
1652    BC1 = -AT*GSphi+(CT+BT*GCchi)*GCphi
1653    BC2 = DT-BT*GSchi
1654    BC3 = AT*GCphi+(CT+BT*GCchi)*GSphi
1655     
1656    BC = BC1*SComeg*SCchi+BC2*SComeg*SSchi-BC3*SSomeg     
1657    psi = acosd(BC)
1658   
1659    BD = 1.0-BC**2
1660    C = np.where(BD>1.e-6,rpd/np.sqrt(BD),0.)
1661    dPSdA = [-C*(-BC1*SSomeg*SCchi-BC2*SSomeg*SSchi-BC3*SComeg),
1662        -C*(-BC1*SComeg*SSchi+BC2*SComeg*SCchi),
1663        -C*(-BC1*SSomeg-BC3*SComeg*SCchi)]
1664     
1665    BA = -BC1*SSchi+BC2*SCchi
1666    BB = BC1*SSomeg*SCchi+BC2*SSomeg*SSchi+BC3*SComeg
1667    gam = atan2d(BB,BA)
1668
1669    BD = (BA**2+BB**2)/rpd
1670
1671    dBAdO = 0
1672    dBAdC = -BC1*SCchi-BC2*SSchi
1673    dBAdF = BC3*SSchi
1674   
1675    dBBdO = BC1*SComeg*SCchi+BC2*SComeg*SSchi-BC3*SSomeg
1676    dBBdC = -BC1*SSomeg*SSchi+BC2*SSomeg*SCchi
1677    dBBdF = BC1*SComeg-BC3*SSomeg*SCchi
1678   
1679    dGMdA = np.where(BD > 1.e-6,[(BA*dBBdO-BB*dBAdO)/BD,(BA*dBBdC-BB*dBAdC)/BD, \
1680        (BA*dBBdF-BB*dBAdF)/BD],[np.zeros_like(BD),np.zeros_like(BD),np.zeros_like(BD)])
1681       
1682    return psi,gam,dPSdA,dGMdA
1683
1684BOH = {
1685'L=2':[[],[],[]],
1686'L=4':[[0.30469720,0.36418281],[],[]],
1687'L=6':[[-0.14104740,0.52775103],[],[]],
1688'L=8':[[0.28646862,0.21545346,0.32826995],[],[]],
1689'L=10':[[-0.16413497,0.33078546,0.39371345],[],[]],
1690'L=12':[[0.26141975,0.27266871,0.03277460,0.32589402],
1691    [0.09298802,-0.23773812,0.49446631,0.0],[]],
1692'L=14':[[-0.17557309,0.25821932,0.27709173,0.33645360],[],[]],
1693'L=16':[[0.24370673,0.29873515,0.06447688,0.00377,0.32574495],
1694    [0.12039646,-0.25330128,0.23950998,0.40962508,0.0],[]],
1695'L=18':[[-0.16914245,0.17017340,0.34598142,0.07433932,0.32696037],
1696    [-0.06901768,0.16006562,-0.24743528,0.47110273,0.0],[]],
1697'L=20':[[0.23067026,0.31151832,0.09287682,0.01089683,0.00037564,0.32573563],
1698    [0.13615420,-0.25048007,0.12882081,0.28642879,0.34620433,0.0],[]],
1699'L=22':[[-0.16109560,0.10244188,0.36285175,0.13377513,0.01314399,0.32585583],
1700    [-0.09620055,0.20244115,-0.22389483,0.17928946,0.42017231,0.0],[]],
1701'L=24':[[0.22050742,0.31770654,0.11661736,0.02049853,0.00150861,0.00003426,0.32573505],
1702    [0.13651722,-0.21386648,0.00522051,0.33939435,0.10837396,0.32914497,0.0],
1703    [0.05378596,-0.11945819,0.16272298,-0.26449730,0.44923956,0.0,0.0]],
1704'L=26':[[-0.15435003,0.05261630,0.35524646,0.18578869,0.03259103,0.00186197,0.32574594],
1705    [-0.11306511,0.22072681,-0.18706142,0.05439948,0.28122966,0.35634355,0.0],[]],
1706'L=28':[[0.21225019,0.32031716,0.13604702,0.03132468,0.00362703,0.00018294,0.00000294,0.32573501],
1707    [0.13219496,-0.17206256,-0.08742608,0.32671661,0.17973107,0.02567515,0.32619598,0.0],
1708    [0.07989184,-0.16735346,0.18839770,-0.20705337,0.12926808,0.42715602,0.0,0.0]],
1709'L=30':[[-0.14878368,0.01524973,0.33628434,0.22632587,0.05790047,0.00609812,0.00022898,0.32573594],
1710    [-0.11721726,0.20915005,-0.11723436,-0.07815329,0.31318947,0.13655742,0.33241385,0.0],
1711    [-0.04297703,0.09317876,-0.11831248,0.17355132,-0.28164031,0.42719361,0.0,0.0]],
1712'L=32':[[0.20533892,0.32087437,0.15187897,0.04249238,0.00670516,0.00054977,0.00002018,0.00000024,0.32573501],
1713    [0.12775091,-0.13523423,-0.14935701,0.28227378,0.23670434,0.05661270,0.00469819,0.32578978,0.0],
1714    [0.09703829,-0.19373733,0.18610682,-0.14407046,0.00220535,0.26897090,0.36633402,0.0,0.0]],
1715'L=34':[[-0.14409234,-0.01343681,0.31248977,0.25557722,0.08571889,0.01351208,0.00095792,0.00002550,0.32573508],
1716    [-0.11527834,0.18472133,-0.04403280,-0.16908618,0.27227021,0.21086614,0.04041752,0.32688152,0.0],
1717    [-0.06773139,0.14120811,-0.15835721,0.18357456,-0.19364673,0.08377174,0.43116318,0.0,0.0]]
1718}
1719
1720Lnorm = lambda L: 4.*np.pi/(2.0*L+1.)
1721
1722def GetKcl(L,N,SGLaue,phi,beta):
1723    'needs doc string'
1724    import pytexture as ptx
1725    if SGLaue in ['m3','m3m']:
1726        if 'array' in str(type(phi)) and np.any(phi.shape):
1727            Kcl = np.zeros_like(phi)
1728        else:
1729            Kcl = 0.
1730        for j in range(0,L+1,4):
1731            im = j/4
1732            if 'array' in str(type(phi)) and np.any(phi.shape):
1733                pcrs = ptx.pyplmpsi(L,j,len(phi),phi)[0]
1734            else:
1735                pcrs = ptx.pyplmpsi(L,j,1,phi)[0]
1736            Kcl += BOH['L=%d'%(L)][N-1][im]*pcrs*cosd(j*beta)       
1737    else:
1738        if 'array' in str(type(phi)) and np.any(phi.shape):
1739            pcrs = ptx.pyplmpsi(L,N,len(phi),phi)[0]
1740        else:
1741            pcrs = ptx.pyplmpsi(L,N,1,phi)[0]
1742        pcrs *= RSQ2PI
1743        if N:
1744            pcrs *= SQ2
1745        if SGLaue in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']:
1746            if SGLaue in ['3mR','3m1','31m']: 
1747                if N%6 == 3:
1748                    Kcl = pcrs*sind(N*beta)
1749                else:
1750                    Kcl = pcrs*cosd(N*beta)
1751            else:
1752                Kcl = pcrs*cosd(N*beta)
1753        else:
1754            Kcl = pcrs*(cosd(N*beta)+sind(N*beta))
1755    return Kcl
1756   
1757def GetKsl(L,M,SamSym,psi,gam):
1758    'needs doc string'
1759    import pytexture as ptx
1760    if 'array' in str(type(psi)) and np.any(psi.shape):
1761        psrs,dpdps = ptx.pyplmpsi(L,M,len(psi),psi)
1762    else:
1763        psrs,dpdps = ptx.pyplmpsi(L,M,1,psi)
1764    psrs *= RSQ2PI
1765    dpdps *= RSQ2PI
1766    if M:
1767        psrs *= SQ2
1768        dpdps *= SQ2
1769    if SamSym in ['mmm',]:
1770        dum = cosd(M*gam)
1771        Ksl = psrs*dum
1772        dKsdp = dpdps*dum
1773        dKsdg = -psrs*M*sind(M*gam)
1774    else:
1775        dum = cosd(M*gam)+sind(M*gam)
1776        Ksl = psrs*dum
1777        dKsdp = dpdps*dum
1778        dKsdg = psrs*M*(-sind(M*gam)+cosd(M*gam))
1779    return Ksl,dKsdp,dKsdg
1780   
1781def GetKclKsl(L,N,SGLaue,psi,phi,beta):
1782    """
1783    This is used for spherical harmonics description of preferred orientation;
1784        cylindrical symmetry only (M=0) and no sample angle derivatives returned
1785    """
1786    import pytexture as ptx
1787    Ksl,x = ptx.pyplmpsi(L,0,1,psi)
1788    Ksl *= RSQ2PI
1789    if SGLaue in ['m3','m3m']:
1790        Kcl = 0.0
1791        for j in range(0,L+1,4):
1792            im = j/4
1793            pcrs,dum = ptx.pyplmpsi(L,j,1,phi)
1794            Kcl += BOH['L=%d'%(L)][N-1][im]*pcrs*cosd(j*beta)       
1795    else:
1796        pcrs,dum = ptx.pyplmpsi(L,N,1,phi)
1797        pcrs *= RSQ2PI
1798        if N:
1799            pcrs *= SQ2
1800        if SGLaue in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']:
1801            if SGLaue in ['3mR','3m1','31m']: 
1802                if N%6 == 3:
1803                    Kcl = pcrs*sind(N*beta)
1804                else:
1805                    Kcl = pcrs*cosd(N*beta)
1806            else:
1807                Kcl = pcrs*cosd(N*beta)
1808        else:
1809            Kcl = pcrs*(cosd(N*beta)+sind(N*beta))
1810    return Kcl*Ksl,Lnorm(L)
1811   
1812def Glnh(Start,SHCoef,psi,gam,SamSym):
1813    'needs doc string'
1814    import pytexture as ptx
1815
1816    if Start:
1817        ptx.pyqlmninit()
1818        Start = False
1819    Fln = np.zeros(len(SHCoef))
1820    for i,term in enumerate(SHCoef):
1821        l,m,n = eval(term.strip('C'))
1822        pcrs,dum = ptx.pyplmpsi(l,m,1,psi)
1823        pcrs *= RSQPI
1824        if m == 0:
1825            pcrs /= SQ2
1826        if SamSym in ['mmm',]:
1827            Ksl = pcrs*cosd(m*gam)
1828        else:
1829            Ksl = pcrs*(cosd(m*gam)+sind(m*gam))
1830        Fln[i] = SHCoef[term]*Ksl*Lnorm(l)
1831    ODFln = dict(zip(SHCoef.keys(),list(zip(SHCoef.values(),Fln))))
1832    return ODFln
1833
1834def Flnh(Start,SHCoef,phi,beta,SGData):
1835    'needs doc string'
1836    import pytexture as ptx
1837   
1838    if Start:
1839        ptx.pyqlmninit()
1840        Start = False
1841    Fln = np.zeros(len(SHCoef))
1842    for i,term in enumerate(SHCoef):
1843        l,m,n = eval(term.strip('C'))
1844        if SGData['SGLaue'] in ['m3','m3m']:
1845            Kcl = 0.0
1846            for j in range(0,l+1,4):
1847                im = j/4
1848                pcrs,dum = ptx.pyplmpsi(l,j,1,phi)
1849                Kcl += BOH['L='+str(l)][n-1][im]*pcrs*cosd(j*beta)       
1850        else:                #all but cubic
1851            pcrs,dum = ptx.pyplmpsi(l,n,1,phi)
1852            pcrs *= RSQPI
1853            if n == 0:
1854                pcrs /= SQ2
1855            if SGData['SGLaue'] in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']:
1856               if SGData['SGLaue'] in ['3mR','3m1','31m']: 
1857                   if n%6 == 3:
1858                       Kcl = pcrs*sind(n*beta)
1859                   else:
1860                       Kcl = pcrs*cosd(n*beta)
1861               else:
1862                   Kcl = pcrs*cosd(n*beta)
1863            else:
1864                Kcl = pcrs*(cosd(n*beta)+sind(n*beta))
1865        Fln[i] = SHCoef[term]*Kcl*Lnorm(l)
1866    ODFln = dict(zip(SHCoef.keys(),list(zip(SHCoef.values(),Fln))))
1867    return ODFln
1868   
1869def polfcal(ODFln,SamSym,psi,gam):
1870    '''Perform a pole figure computation.
1871    Note that the the number of gam values must either be 1 or must
1872    match psi. Updated for numpy 1.8.0
1873    '''
1874    import pytexture as ptx
1875    PolVal = np.ones_like(psi)
1876    for term in ODFln:
1877        if abs(ODFln[term][1]) > 1.e-3:
1878            l,m,n = eval(term.strip('C'))
1879            psrs,dum = ptx.pyplmpsi(l,m,len(psi),psi)
1880            if SamSym in ['-1','2/m']:
1881                if m:
1882                    Ksl = RSQPI*psrs*(cosd(m*gam)+sind(m*gam))
1883                else:
1884                    Ksl = RSQPI*psrs/SQ2
1885            else:
1886                if m:
1887                    Ksl = RSQPI*psrs*cosd(m*gam)
1888                else:
1889                    Ksl = RSQPI*psrs/SQ2
1890            PolVal += ODFln[term][1]*Ksl
1891    return PolVal
1892   
1893def invpolfcal(ODFln,SGData,phi,beta):
1894    'needs doc string'
1895    import pytexture as ptx
1896   
1897    invPolVal = np.ones_like(beta)
1898    for term in ODFln:
1899        if abs(ODFln[term][1]) > 1.e-3:
1900            l,m,n = eval(term.strip('C'))
1901            if SGData['SGLaue'] in ['m3','m3m']:
1902                Kcl = 0.0
1903                for j in range(0,l+1,4):
1904                    im = j/4
1905                    pcrs,dum = ptx.pyplmpsi(l,j,len(beta),phi)
1906                    Kcl += BOH['L=%d'%(l)][n-1][im]*pcrs*cosd(j*beta)       
1907            else:                #all but cubic
1908                pcrs,dum = ptx.pyplmpsi(l,n,len(beta),phi)
1909                pcrs *= RSQPI
1910                if n == 0:
1911                    pcrs /= SQ2
1912                if SGData['SGLaue'] in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']:
1913                   if SGData['SGLaue'] in ['3mR','3m1','31m']: 
1914                       if n%6 == 3:
1915                           Kcl = pcrs*sind(n*beta)
1916                       else:
1917                           Kcl = pcrs*cosd(n*beta)
1918                   else:
1919                       Kcl = pcrs*cosd(n*beta)
1920                else:
1921                    Kcl = pcrs*(cosd(n*beta)+sind(n*beta))
1922            invPolVal += ODFln[term][1]*Kcl
1923    return invPolVal
1924   
1925   
1926def textureIndex(SHCoef):
1927    'needs doc string'
1928    Tindx = 1.0
1929    for term in SHCoef:
1930        l = eval(term.strip('C'))[0]
1931        Tindx += SHCoef[term]**2/(2.0*l+1.)
1932    return Tindx
1933   
1934# self-test materials follow.
1935selftestlist = []
1936'''Defines a list of self-tests'''
1937selftestquiet = True
1938def _ReportTest():
1939    'Report name and doc string of current routine when ``selftestquiet`` is False'
1940    if not selftestquiet:
1941        import inspect
1942        caller = inspect.stack()[1][3]
1943        doc = eval(caller).__doc__
1944        if doc is not None:
1945            print('testing '+__file__+' with '+caller+' ('+doc+')')
1946        else:
1947            print('testing '+__file__()+" with "+caller)
1948NeedTestData = True
1949def TestData():
1950    array = np.array
1951    global NeedTestData
1952    NeedTestData = False
1953    global CellTestData
1954    # output from uctbx computed on platform darwin on 2010-05-28
1955    CellTestData = [
1956# cell, g, G, cell*, V, V*
1957  [(4, 4, 4, 90, 90, 90), 
1958   array([[  1.60000000e+01,   9.79717439e-16,   9.79717439e-16],
1959       [  9.79717439e-16,   1.60000000e+01,   9.79717439e-16],
1960       [  9.79717439e-16,   9.79717439e-16,   1.60000000e+01]]), array([[  6.25000000e-02,   3.82702125e-18,   3.82702125e-18],
1961       [  3.82702125e-18,   6.25000000e-02,   3.82702125e-18],
1962       [  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],
1963# cell, g, G, cell*, V, V*
1964  [(4.0999999999999996, 5.2000000000000002, 6.2999999999999998, 100, 80, 130), 
1965   array([[ 16.81      , -13.70423184,   4.48533243],
1966       [-13.70423184,  27.04      ,  -5.6887143 ],
1967       [  4.48533243,  -5.6887143 ,  39.69      ]]), array([[ 0.10206349,  0.05083339, -0.00424823],
1968       [ 0.05083339,  0.06344997,  0.00334956],
1969       [-0.00424823,  0.00334956,  0.02615544]]), (0.31947376387537696, 0.25189277536327803, 0.16172643497798223, 85.283666420376008, 94.716333579624006, 50.825714168082683), 100.98576357983838, 0.0099023858863968445],
1970# cell, g, G, cell*, V, V*
1971  [(3.5, 3.5, 6, 90, 90, 120), 
1972   array([[  1.22500000e+01,  -6.12500000e+00,   1.28587914e-15],
1973       [ -6.12500000e+00,   1.22500000e+01,   1.28587914e-15],
1974       [  1.28587914e-15,   1.28587914e-15,   3.60000000e+01]]), array([[  1.08843537e-01,   5.44217687e-02,   3.36690552e-18],
1975       [  5.44217687e-02,   1.08843537e-01,   3.36690552e-18],
1976       [  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],
1977  ]
1978    global CoordTestData
1979    CoordTestData = [
1980# cell, ((frac, ortho),...)
1981  ((4,4,4,90,90,90,), [
1982 ((0.10000000000000001, 0.0, 0.0),(0.40000000000000002, 0.0, 0.0)),
1983 ((0.0, 0.10000000000000001, 0.0),(2.4492935982947065e-17, 0.40000000000000002, 0.0)),
1984 ((0.0, 0.0, 0.10000000000000001),(2.4492935982947065e-17, -2.4492935982947065e-17, 0.40000000000000002)),
1985 ((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(0.40000000000000013, 0.79999999999999993, 1.2)),
1986 ((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(0.80000000000000016, 1.2, 0.40000000000000002)),
1987 ((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(1.2, 0.80000000000000004, 0.40000000000000002)),
1988 ((0.5, 0.5, 0.5),(2.0, 1.9999999999999998, 2.0)),
1989]),
1990# cell, ((frac, ortho),...)
1991  ((4.1,5.2,6.3,100,80,130,), [
1992 ((0.10000000000000001, 0.0, 0.0),(0.40999999999999998, 0.0, 0.0)),
1993 ((0.0, 0.10000000000000001, 0.0),(-0.33424955703700043, 0.39834311042186865, 0.0)),
1994 ((0.0, 0.0, 0.10000000000000001),(0.10939835193016617, -0.051013289294572106, 0.6183281045774256)),
1995 ((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(0.069695941716497567, 0.64364635296002093, 1.8549843137322766)),
1996 ((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(-0.073350319180835066, 1.1440160419710339, 0.6183281045774256)),
1997 ((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(0.67089923785616512, 0.74567293154916525, 0.6183281045774256)),
1998 ((0.5, 0.5, 0.5),(0.92574397446582857, 1.7366491056364828, 3.0916405228871278)),
1999]),
2000# cell, ((frac, ortho),...)
2001  ((3.5,3.5,6,90,90,120,), [
2002 ((0.10000000000000001, 0.0, 0.0),(0.35000000000000003, 0.0, 0.0)),
2003 ((0.0, 0.10000000000000001, 0.0),(-0.17499999999999993, 0.3031088913245536, 0.0)),
2004 ((0.0, 0.0, 0.10000000000000001),(3.6739403974420595e-17, -3.6739403974420595e-17, 0.60000000000000009)),
2005 ((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(2.7675166561703527e-16, 0.60621778264910708, 1.7999999999999998)),
2006 ((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(0.17500000000000041, 0.90932667397366063, 0.60000000000000009)),
2007 ((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(0.70000000000000018, 0.6062177826491072, 0.60000000000000009)),
2008 ((0.5, 0.5, 0.5),(0.87500000000000067, 1.5155444566227676, 3.0)),
2009]),
2010]
2011    global LaueTestData             #generated by GSAS
2012    LaueTestData = {
2013    '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),
2014        (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),
2015        (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))],
2016    '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),
2017        (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),
2018        (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),
2019        (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))],
2020    '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),
2021        (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),
2022        (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),
2023        (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),
2024        (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),
2025        (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),
2026        (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),
2027        (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),
2028        (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),
2029        (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),
2030        (2,1,5,6),(2,1,-5,6),(3,-1,-5,6))],
2031    '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),
2032        (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),
2033        (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),
2034        (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),
2035        (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),
2036        (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),
2037        (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),
2038        (3,0,4,6),(3,1,-2,12),(3,0,-4,6),(1,1,6,12),(2,2,3,12))],
2039    '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),
2040        (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),
2041        (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),
2042        (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),
2043        (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),
2044        (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),
2045        (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))],
2046    }
2047   
2048    global FLnhTestData
2049    FLnhTestData = [{
2050    'C(4,0,0)': (0.965, 0.42760447),
2051    'C(2,0,0)': (1.0122, -0.80233610),
2052    'C(2,0,2)': (0.0061, 8.37491546E-03),
2053    'C(6,0,4)': (-0.0898, 4.37985696E-02),
2054    'C(6,0,6)': (-0.1369, -9.04081762E-02),
2055    'C(6,0,0)': (0.5935, -0.18234928),
2056    'C(4,0,4)': (0.1872, 0.16358127),
2057    'C(6,0,2)': (0.6193, 0.27573633),
2058    'C(4,0,2)': (-0.1897, 0.12530720)},[1,0,0]]
2059def test0():
2060    if NeedTestData: TestData()
2061    msg = 'test cell2Gmat, fillgmat, Gmat2cell'
2062    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2063        G, g = cell2Gmat(cell)
2064        assert np.allclose(G,tG),msg
2065        assert np.allclose(g,tg),msg
2066        tcell = Gmat2cell(g)
2067        assert np.allclose(cell,tcell),msg
2068        tcell = Gmat2cell(G)
2069        assert np.allclose(tcell,trcell),msg
2070selftestlist.append(test0)
2071
2072def test1():
2073    'test cell2A and A2Gmat'
2074    _ReportTest()
2075    if NeedTestData: TestData()
2076    msg = 'test cell2A and A2Gmat'
2077    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2078        G, g = A2Gmat(cell2A(cell))
2079        assert np.allclose(G,tG),msg
2080        assert np.allclose(g,tg),msg
2081selftestlist.append(test1)
2082
2083def test2():
2084    'test Gmat2A, A2cell, A2Gmat, Gmat2cell'
2085    _ReportTest()
2086    if NeedTestData: TestData()
2087    msg = 'test Gmat2A, A2cell, A2Gmat, Gmat2cell'
2088    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2089        G, g = cell2Gmat(cell)
2090        tcell = A2cell(Gmat2A(G))
2091        assert np.allclose(cell,tcell),msg
2092selftestlist.append(test2)
2093
2094def test3():
2095    'test invcell2Gmat'
2096    _ReportTest()
2097    if NeedTestData: TestData()
2098    msg = 'test invcell2Gmat'
2099    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2100        G, g = invcell2Gmat(trcell)
2101        assert np.allclose(G,tG),msg
2102        assert np.allclose(g,tg),msg
2103selftestlist.append(test3)
2104
2105def test4():
2106    'test calc_rVsq, calc_rV, calc_V'
2107    _ReportTest()
2108    if NeedTestData: TestData()
2109    msg = 'test calc_rVsq, calc_rV, calc_V'
2110    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2111        assert np.allclose(calc_rV(cell2A(cell)),trV), msg
2112        assert np.allclose(calc_V(cell2A(cell)),tV), msg
2113selftestlist.append(test4)
2114
2115def test5():
2116    'test A2invcell'
2117    _ReportTest()
2118    if NeedTestData: TestData()
2119    msg = 'test A2invcell'
2120    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2121        rcell = A2invcell(cell2A(cell))
2122        assert np.allclose(rcell,trcell),msg
2123selftestlist.append(test5)
2124
2125def test6():
2126    'test cell2AB'
2127    _ReportTest()
2128    if NeedTestData: TestData()
2129    msg = 'test cell2AB'
2130    for (cell,coordlist) in CoordTestData:
2131        A,B = cell2AB(cell)
2132        for (frac,ortho) in coordlist:
2133            to = np.inner(A,frac)
2134            tf = np.inner(B,to)
2135            assert np.allclose(ortho,to), msg
2136            assert np.allclose(frac,tf), msg
2137            to = np.sum(A*frac,axis=1)
2138            tf = np.sum(B*to,axis=1)
2139            assert np.allclose(ortho,to), msg
2140            assert np.allclose(frac,tf), msg
2141selftestlist.append(test6)
2142
2143def test7():
2144    'test GetBraviasNum(...) and GenHBravais(...)'
2145    _ReportTest()
2146    import os.path
2147    import sys
2148    import GSASIIspc as spc
2149    testdir = os.path.join(os.path.split(os.path.abspath( __file__ ))[0],'testinp')
2150    if os.path.exists(testdir):
2151        if testdir not in sys.path: sys.path.insert(0,testdir)
2152    import sgtbxlattinp
2153    derror = 1e-4
2154    def indexmatch(hklin, hkllist, system):
2155        for hklref in hkllist:
2156            hklref = list(hklref)
2157            # these permutations are far from complete, but are sufficient to
2158            # allow the test to complete
2159            if system == 'cubic':
2160                permlist = [(1,2,3),(1,3,2),(2,1,3),(2,3,1),(3,1,2),(3,2,1),]
2161            elif system == 'monoclinic':
2162                permlist = [(1,2,3),(-1,2,-3)]
2163            else:
2164                permlist = [(1,2,3)]
2165
2166            for perm in permlist:
2167                hkl = [abs(i) * hklin[abs(i)-1] / i for i in perm]
2168                if hkl == hklref: return True
2169                if [-i for i in hkl] == hklref: return True
2170        else:
2171            return False
2172
2173    for key in sgtbxlattinp.sgtbx7:
2174        spdict = spc.SpcGroup(key)
2175        cell = sgtbxlattinp.sgtbx7[key][0]
2176        system = spdict[1]['SGSys']
2177        center = spdict[1]['SGLatt']
2178
2179        bravcode = GetBraviasNum(center, system)
2180
2181        g2list = GenHBravais(sgtbxlattinp.dmin, bravcode, cell2A(cell))
2182
2183        assert len(sgtbxlattinp.sgtbx7[key][1]) == len(g2list), 'Reflection lists differ for %s' % key
2184        for h,k,l,d,num in g2list:
2185            for hkllist,dref in sgtbxlattinp.sgtbx7[key][1]: 
2186                if abs(d-dref) < derror:
2187                    if indexmatch((h,k,l,), hkllist, system):
2188                        break
2189            else:
2190                assert 0,'No match for %s at %s (%s)' % ((h,k,l),d,key)
2191selftestlist.append(test7)
2192
2193def test8():
2194    'test GenHLaue'
2195    _ReportTest()
2196    import GSASIIspc as spc
2197    import sgtbxlattinp
2198    derror = 1e-4
2199    dmin = sgtbxlattinp.dmin
2200
2201    def indexmatch(hklin, hklref, system, axis):
2202        # these permutations are far from complete, but are sufficient to
2203        # allow the test to complete
2204        if system == 'cubic':
2205            permlist = [(1,2,3),(1,3,2),(2,1,3),(2,3,1),(3,1,2),(3,2,1),]
2206        elif system == 'monoclinic' and axis=='b':
2207            permlist = [(1,2,3),(-1,2,-3)]
2208        elif system == 'monoclinic' and axis=='a':
2209            permlist = [(1,2,3),(1,-2,-3)]
2210        elif system == 'monoclinic' and axis=='c':
2211            permlist = [(1,2,3),(-1,-2,3)]
2212        elif system == 'trigonal':
2213            permlist = [(1,2,3),(2,1,3),(-1,-2,3),(-2,-1,3)]
2214        elif system == 'rhombohedral':
2215            permlist = [(1,2,3),(2,3,1),(3,1,2)]
2216        else:
2217            permlist = [(1,2,3)]
2218
2219        hklref = list(hklref)
2220        for perm in permlist:
2221            hkl = [abs(i) * hklin[abs(i)-1] / i for i in perm]
2222            if hkl == hklref: return True
2223            if [-i for i in hkl] == hklref: return True
2224        return False
2225
2226    for key in sgtbxlattinp.sgtbx8:
2227        spdict = spc.SpcGroup(key)[1]
2228        cell = sgtbxlattinp.sgtbx8[key][0]
2229        Axis = spdict['SGUniq']
2230        system = spdict['SGSys']
2231
2232        g2list = GenHLaue(dmin,spdict,cell2A(cell))
2233        #if len(g2list) != len(sgtbxlattinp.sgtbx8[key][1]):
2234        #    print 'failed',key,':' ,len(g2list),'vs',len(sgtbxlattinp.sgtbx8[key][1])
2235        #    print 'GSAS-II:'
2236        #    for h,k,l,d in g2list: print '  ',(h,k,l),d
2237        #    print 'SGTBX:'
2238        #    for hkllist,dref in sgtbxlattinp.sgtbx8[key][1]: print '  ',hkllist,dref
2239        assert len(g2list) == len(sgtbxlattinp.sgtbx8[key][1]), (
2240            'Reflection lists differ for %s' % key
2241            )
2242        #match = True
2243        for h,k,l,d in g2list:
2244            for hkllist,dref in sgtbxlattinp.sgtbx8[key][1]: 
2245                if abs(d-dref) < derror:
2246                    if indexmatch((h,k,l,), hkllist, system, Axis): break
2247            else:
2248                assert 0,'No match for %s at %s (%s)' % ((h,k,l),d,key)
2249                #match = False
2250        #if not match:
2251            #for hkllist,dref in sgtbxlattinp.sgtbx8[key][1]: print '  ',hkllist,dref
2252            #print center, Laue, Axis, system
2253selftestlist.append(test8)
2254           
2255def test9():
2256    'test GenHLaue'
2257    _ReportTest()
2258    import GSASIIspc as G2spc
2259    if NeedTestData: TestData()
2260    for spc in LaueTestData:
2261        data = LaueTestData[spc]
2262        cell = data[0]
2263        hklm = np.array(data[1])
2264        H = hklm[-1][:3]
2265        hklO = hklm.T[:3].T
2266        A = cell2A(cell)
2267        dmin = 1./np.sqrt(calc_rDsq(H,A))
2268        SGData = G2spc.SpcGroup(spc)[1]
2269        hkls = np.array(GenHLaue(dmin,SGData,A))
2270        hklN = hkls.T[:3].T
2271        #print spc,hklO.shape,hklN.shape
2272        err = True
2273        for H in hklO:
2274            if H not in hklN:
2275                print H,' missing from hkl from GSASII'
2276                err = False
2277        assert(err)
2278selftestlist.append(test9)
2279       
2280       
2281   
2282
2283if __name__ == '__main__':
2284    # run self-tests
2285    selftestquiet = False
2286    for test in selftestlist:
2287        test()
2288    print "OK"
Note: See TracBrowser for help on using the repository browser.