# source:trunk/GSASIIlattice.py@4880

Last change on this file since 4880 was 4880, checked in by toby, 20 months ago

implement constraints in scriptable, part 3; misc docs fixes

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