source: trunk/GSASIIlattice.py @ 2560

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

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