source: trunk/GSASIIlattice.py @ 2473

Last change on this file since 2473 was 2473, checked in by vondreele, 9 years ago

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