source: trunk/GSASIIlattice.py @ 2486

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

fix mag form factor lookup - some aren't chemical valences (e.g. Fe+4)
fix mag form factor periodic table - nonmagnetic atoms now not shown
fix powder reflection mark & diff curve placement issues
fix issues with mag structure drawings - fill unit cell, etc.
fix structure transform to make magnetic cells (still needs making of constraints)
fix problem of mustrain & hstrain coeff changes with change in space group
add print to file of covariance matrix - user request
fix magnetic moment site symmetry restrictions
work on mag structure factor calcs.
change EXP file importer to ignore magnetic symmetry from GSAS

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