source: trunk/GSASIIlattice.py @ 2512

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

remove 'constraint' and 'restraint' from Consraints & Restraints tabs
work on automatic nucl-mag constraints
allow selection of possible mag atoms to use in new mag phase
fix issues with chemical composition restraints

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