source: trunk/GSASIIlattice.py @ 2481

Last change on this file since 2481 was 2481, checked in by vondreele, 5 years ago

force Transform to delete nonmagnetic atoms if phase made magnetic & add 'mag' to new phase name
fix TOF cosine background function
magnetic structure refinement works (in numeric mode only)
magnetic structure factor calculation correct

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