source: trunk/GSASIIlattice.py @ 3092

Last change on this file since 3092 was 3092, checked in by vondreele, 4 years ago

change # decimal places on pseudovariables from 3 to 5
fix bug on phase Transform command
clean up formatting of General tab

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