# source:trunk/GSASIIlattice.py@4792

Last change on this file since 4792 was 4792, checked in by toby, 2 years ago

check in changes to allow cctbx to use DoIndexPeaks and !refinePeaks

• Property svn:eol-style set to native
• Property svn:keywords set to Date Author Revision URL Id
File size: 106.9 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-02-02 17:01:51 +0000 (Tue, 02 Feb 2021)$
36# $Author: toby$
37# $Revision: 4792$
38# $URL: trunk/GSASIIlattice.py$
39# $Id: GSASIIlattice.py 4792 2021-02-02 17:01:51Z 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: 4792$")
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
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] or 'PKS' in Inst['Type'][0]:
845        return TOF2dsp(Inst,pos)
846    else:   #'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] or 'PKS' 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(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)
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 GetBraviasNum(center,system):
1073    """Determine the Bravais lattice number, as used in GenHBravais
1074
1075    :param center: one of: 'P', 'C', 'I', 'F', 'R' (see SGLatt from GSASIIspc.SpcGroup)
1076    :param system: one of 'cubic', 'hexagonal', 'tetragonal', 'orthorhombic', 'trigonal' (for R)
1077      'monoclinic', 'triclinic' (see SGSys from GSASIIspc.SpcGroup)
1078    :return: a number between 0 and 13
1079      or throws a ValueError exception if the combination of center, system is not found (i.e. non-standard)
1080
1081    """
1082    if center.upper() == 'F' and system.lower() == 'cubic':
1083        return 0
1084    elif center.upper() == 'I' and system.lower() == 'cubic':
1085        return 1
1086    elif center.upper() == 'P' and system.lower() == 'cubic':
1087        return 2
1088    elif center.upper() == 'R' and system.lower() == 'trigonal':
1089        return 3
1090    elif center.upper() == 'P' and system.lower() == 'hexagonal':
1091        return 4
1092    elif center.upper() == 'I' and system.lower() == 'tetragonal':
1093        return 5
1094    elif center.upper() == 'P' and system.lower() == 'tetragonal':
1095        return 6
1096    elif center.upper() == 'F' and system.lower() == 'orthorhombic':
1097        return 7
1098    elif center.upper() == 'I' and system.lower() == 'orthorhombic':
1099        return 8
1100    elif center.upper() == 'A' and system.lower() == 'orthorhombic':
1101        return 9
1102    elif center.upper() == 'B' and system.lower() == 'orthorhombic':
1103        return 10
1104    elif center.upper() == 'C' and system.lower() == 'orthorhombic':
1105        return 11
1106    elif center.upper() == 'P' and system.lower() == 'orthorhombic':
1107        return 12
1108    elif center.upper() == 'C' and system.lower() == 'monoclinic':
1109        return 13
1110    elif center.upper() == 'P' and system.lower() == 'monoclinic':
1111        return 14
1112    elif center.upper() == 'P' and system.lower() == 'triclinic':
1113        return 15
1114    raise ValueError('non-standard Bravais lattice center=%s, cell=%s' % (center,system))
1115
1116def _GenHBravais_cctbx(dmin, Bravais, A, sg_type, uctbx_unit_cell, miller_index_generator):
1117    '''Alternate form of :func:GenHBravais that uses CCTBX internals
1118    '''
1119    g_inv = np.array([[A[0],   A[3]/2, A[4]/2],
1120                      [A[3]/2, A[1],   A[5]/2],
1121                      [A[4]/2, A[5]/2, A[2]]])
1122    g = np.linalg.inv(g_inv)
1123    g_elems = (g[0][0], g[1][1], g[2][2], g[0][1], g[0][2], g[1][2])
1124    try:
1125        uc = uctbx_unit_cell(metrical_matrix=g_elems)
1126    except ValueError: # this function sometimes receives an A matrix that gives
1127                           # numbers <0 in the diagonal elems of g. Not sure why.
1128        return []
1129    #if sg_type is None:
1130    #    sg_type = make_sgtype(Bravais)
1131    mig = miller_index_generator(uc, sg_type, 0, dmin)
1132    result = []
1133    for h,k,l in mig:
1134        d = uc.d((h,k,l))
1135        result.append([h, k, l, d, -1])
1136    result.sort(key=lambda l: l[3], reverse=True)
1137    return result
1138
1139def GenHBravais(dmin, Bravais, A, cctbx_args=None):
1140    """Generate the positionally unique powder diffraction reflections
1141
1142    :param dmin: minimum d-spacing in A
1143    :param Bravais: lattice type (see GetBraviasNum). Bravais is one of:
1144
1145            * 0 F cubic
1146            * 1 I cubic
1147            * 2 P cubic
1148            * 3 R hexagonal (trigonal not rhombohedral)
1149            * 4 P hexagonal
1150            * 5 I tetragonal
1151            * 6 P tetragonal
1152            * 7 F orthorhombic
1153            * 8 I orthorhombic
1154            * 9 A orthorhombic
1155            * 10 B orthorhombic
1156            * 11 C orthorhombic
1157            * 12 P orthorhombic
1158            * 13 I monoclinic
1159            * 14 C monoclinic
1160            * 15 P monoclinic
1161            * 16 P triclinic
1162
1163    :param A: reciprocal metric tensor elements as [G11,G22,G33,2*G12,2*G13,2*G23]
1164    :param dict cctbx_args: items defined in CCTBX:
1165
1166         * 'sg_type': value from cctbx.sgtbx.space_group_type(symmorphic_sgs[ibrav])
1167         * 'uctbx_unit_cell': pointer to :meth:cctbx.uctbx.unit_cell
1168         * 'miller_index_generator':  pointer to :meth:cctbx.miller.index_generator
1169
1170    :returns: HKL unique d list of [h,k,l,d,-1] sorted with largest d first
1171
1172    """
1173    if cctbx_args:
1174        return _GenHBravais_cctbx(dmin, Bravais, A,
1175                    cctbx_args['sg_type'], cctbx_args['uctbx_unit_cell'], cctbx_args['miller_index_generator'])
1176
1177    if Bravais in [9,]:
1178        Cent = 'A'
1179    elif Bravais in [10,]:
1180        Cent = 'B'
1181    elif Bravais in [11,14]:
1182        Cent = 'C'
1183    elif Bravais in [1,5,8,13]:
1184        Cent = 'I'
1185    elif Bravais in [0,7]:
1186        Cent = 'F'
1187    elif Bravais in [3]:
1188        Cent = 'R'
1189    else:
1190        Cent = 'P'
1191    Hmax = MaxIndex(dmin,A)
1192    dminsq = 1./(dmin**2)
1193    HKL = []
1194    if Bravais == 16:                       #triclinic
1195        for l in range(-Hmax[2],Hmax[2]+1):
1196            for k in range(-Hmax[1],Hmax[1]+1):
1197                hmin = 0
1198                if (k < 0): hmin = 1
1199                if (k ==0 and l < 0): hmin = 1
1200                for h in range(hmin,Hmax[0]+1):
1201                    H=[h,k,l]
1202                    rdsq = calc_rDsq(H,A)
1203                    if 0 < rdsq <= dminsq:
1204                        HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
1205    elif Bravais in [13,14,15]:                #monoclinic - b unique
1206        Hmax = SwapIndx(2,Hmax)
1207        for h in range(Hmax[0]+1):
1208            for k in range(-Hmax[1],Hmax[1]+1):
1209                lmin = 0
1210                if k < 0:lmin = 1
1211                for l in range(lmin,Hmax[2]+1):
1212                    [h,k,l] = SwapIndx(-2,[h,k,l])
1213                    H = []
1214                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
1215                    if H:
1216                        rdsq = calc_rDsq(H,A)
1217                        if 0 < rdsq <= dminsq:
1218                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
1219                    [h,k,l] = SwapIndx(2,[h,k,l])
1220    elif Bravais in [7,8,9,10,11,12]:            #orthorhombic
1221        for h in range(Hmax[0]+1):
1222            for k in range(Hmax[1]+1):
1223                for l in range(Hmax[2]+1):
1224                    H = []
1225                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
1226                    if H:
1227                        rdsq = calc_rDsq(H,A)
1228                        if 0 < rdsq <= dminsq:
1229                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
1230    elif Bravais in [5,6]:                  #tetragonal
1231        for l in range(Hmax[2]+1):
1232            for k in range(Hmax[1]+1):
1233                for h in range(k,Hmax[0]+1):
1234                    H = []
1235                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
1236                    if H:
1237                        rdsq = calc_rDsq(H,A)
1238                        if 0 < rdsq <= dminsq:
1239                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
1240    elif Bravais in [3,4]:
1241        lmin = 0
1242        if Bravais == 3: lmin = -Hmax[2]                  #hexagonal/trigonal
1243        for l in range(lmin,Hmax[2]+1):
1244            for k in range(Hmax[1]+1):
1245                hmin = k
1246                if l < 0: hmin += 1
1247                for h in range(hmin,Hmax[0]+1):
1248                    H = []
1249                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
1250                    if H:
1251                        rdsq = calc_rDsq(H,A)
1252                        if 0 < rdsq <= dminsq:
1253                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
1254
1255    else:                                   #cubic
1256        for l in range(Hmax[2]+1):
1257            for k in range(l,Hmax[1]+1):
1258                for h in range(k,Hmax[0]+1):
1259                    H = []
1260                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
1261                    if H:
1262                        rdsq = calc_rDsq(H,A)
1263                        if 0 < rdsq <= dminsq:
1264                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
1265    return sortHKLd(HKL,True,False)
1266
1267def getHKLmax(dmin,SGData,A):
1268    'finds maximum allowed hkl for given A within dmin'
1269    SGLaue = SGData['SGLaue']
1270    if SGLaue in ['3R','3mR']:        #Rhombohedral axes
1271        Hmax = [0,0,0]
1272        cell = A2cell(A)
1273        aHx = cell[0]*math.sqrt(2.0*(1.0-cosd(cell[3])))
1274        cHx = cell[0]*math.sqrt(3.0*(1.0+2.0*cosd(cell[3])))
1275        Hmax[0] = Hmax[1] = int(round(aHx/dmin))
1276        Hmax[2] = int(round(cHx/dmin))
1277        #print Hmax,aHx,cHx
1278    else:                           # all others
1279        Hmax = MaxIndex(dmin,A)
1280    return Hmax
1281
1282def GenHLaue(dmin,SGData,A):
1283    """Generate the crystallographically unique powder diffraction reflections
1284    for a lattice and Bravais type
1285
1286    :param dmin: minimum d-spacing
1287    :param SGData: space group dictionary with at least
1288
1289        * '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'
1290        * 'SGLatt': lattice centering: one of 'P','A','B','C','I','F'
1291        * 'SGUniq': code for unique monoclinic axis one of 'a','b','c' (only if 'SGLaue' is '2/m') otherwise an empty string
1292
1293    :param A: reciprocal metric tensor elements as [G11,G22,G33,2*G12,2*G13,2*G23]
1294    :return: HKL = list of [h,k,l,d] sorted with largest d first and is unique
1295            part of reciprocal space ignoring anomalous dispersion
1296
1297    """
1298    import math
1299    SGLaue = SGData['SGLaue']
1300    SGLatt = SGData['SGLatt']
1301    SGUniq = SGData['SGUniq']
1302    #finds maximum allowed hkl for given A within dmin
1303    Hmax = getHKLmax(dmin,SGData,A)
1304
1305    dminsq = 1./(dmin**2)
1306    HKL = []
1307    if SGLaue == '-1':                       #triclinic
1308        for l in range(-Hmax[2],Hmax[2]+1):
1309            for k in range(-Hmax[1],Hmax[1]+1):
1310                hmin = 0
1311                if (k < 0) or (k ==0 and l < 0): hmin = 1
1312                for h in range(hmin,Hmax[0]+1):
1313                    H = []
1314                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
1315                    if H:
1316                        rdsq = calc_rDsq(H,A)
1317                        if 0 < rdsq <= dminsq:
1318                            HKL.append([h,k,l,1./math.sqrt(rdsq)])
1319    elif SGLaue == '2/m':                #monoclinic
1320        axisnum = 1 + ['a','b','c'].index(SGUniq)
1321        Hmax = SwapIndx(axisnum,Hmax)
1322        for h in range(Hmax[0]+1):
1323            for k in range(-Hmax[1],Hmax[1]+1):
1324                lmin = 0
1325                if k < 0:lmin = 1
1326                for l in range(lmin,Hmax[2]+1):
1327                    [h,k,l] = SwapIndx(-axisnum,[h,k,l])
1328                    H = []
1329                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
1330                    if H:
1331                        rdsq = calc_rDsq(H,A)
1332                        if 0 < rdsq <= dminsq:
1333                            HKL.append([h,k,l,1./math.sqrt(rdsq)])
1334                    [h,k,l] = SwapIndx(axisnum,[h,k,l])
1335    elif SGLaue in ['mmm','4/m','6/m']:            #orthorhombic
1336        for l in range(Hmax[2]+1):
1337            for h in range(Hmax[0]+1):
1338                kmin = 1
1339                if SGLaue == 'mmm' or h ==0: kmin = 0
1340                for k in range(kmin,Hmax[1]+1):
1341                    H = []
1342                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
1343                    if H:
1344                        rdsq = calc_rDsq(H,A)
1345                        if 0 < rdsq <= dminsq:
1346                            HKL.append([h,k,l,1./math.sqrt(rdsq)])
1347    elif SGLaue in ['4/mmm','6/mmm']:                  #tetragonal & hexagonal
1348        for l in range(Hmax[2]+1):
1349            for h in range(Hmax[0]+1):
1350                for k in range(h+1):
1351                    H = []
1352                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
1353                    if H:
1354                        rdsq = calc_rDsq(H,A)
1355                        if 0 < rdsq <= dminsq:
1356                            HKL.append([h,k,l,1./math.sqrt(rdsq)])
1357    elif SGLaue in ['3m1','31m','3','3R','3mR']:                  #trigonals
1358        for l in range(-Hmax[2],Hmax[2]+1):
1359            hmin = 0
1360            if l < 0: hmin = 1
1361            for h in range(hmin,Hmax[0]+1):
1362                if SGLaue in ['3R','3']:
1363                    kmax = h
1364                    kmin = -int((h-1.)/2.)
1365                else:
1366                    kmin = 0
1367                    kmax = h
1368                    if SGLaue in ['3m1','3mR'] and l < 0: kmax = h-1
1369                    if SGLaue == '31m' and l < 0: kmin = 1
1370                for k in range(kmin,kmax+1):
1371                    H = []
1372                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
1373                    if SGLaue in ['3R','3mR']:
1374                        H = Hx2Rh(H)
1375                    if H:
1376                        rdsq = calc_rDsq(H,A)
1377                        if 0 < rdsq <= dminsq:
1378                            HKL.append([H[0],H[1],H[2],1./math.sqrt(rdsq)])
1379    else:                                   #cubic
1380        for h in range(Hmax[0]+1):
1381            for k in range(h+1):
1382                lmin = 0
1383                lmax = k
1384                if SGLaue =='m3':
1385                    lmax = h-1
1386                    if h == k: lmax += 1
1387                for l in range(lmin,lmax+1):
1388                    H = []
1389                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
1390                    if H:
1391                        rdsq = calc_rDsq(H,A)
1392                        if 0 < rdsq <= dminsq:
1393                            HKL.append([h,k,l,1./math.sqrt(rdsq)])
1394    return sortHKLd(HKL,True,True)
1395
1396def GenPfHKLs(nMax,SGData,A):
1397    """Generate the unique pole figure reflections for a lattice and Bravais type.
1398    Min d-spacing=1.0A & no more than nMax returned
1399
1400    :param nMax: maximum number of hkls returned
1401    :param SGData: space group dictionary with at least
1402
1403        * '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'
1404        * 'SGLatt': lattice centering: one of 'P','A','B','C','I','F'
1405        * 'SGUniq': code for unique monoclinic axis one of 'a','b','c' (only if 'SGLaue' is '2/m') otherwise an empty string
1406
1407    :param A: reciprocal metric tensor elements as [G11,G22,G33,2*G12,2*G13,2*G23]
1408    :return: HKL = list of 'h k l' strings sorted with largest d first; no duplicate zones
1409
1410    """
1411    HKL = np.array(GenHLaue(1.0,SGData,A)).T[:3].T     #strip d-spacings
1412    N = min(nMax,len(HKL))
1413    return ['%d %d %d'%(h[0],h[1],h[2]) for h in HKL[:N]]
1414
1415def GenSSHLaue(dmin,SGData,SSGData,Vec,maxH,A):
1416    'needs a doc string'
1417    ifMag = False
1418    if 'MagSpGrp' in SGData:
1419        ifMag = True
1420    HKLs = []
1421    vec = np.array(Vec)
1422    vstar = np.sqrt(calc_rDsq(vec,A))     #find extra needed for -n SS reflections
1423    dvec = 1./(maxH*vstar+1./dmin)
1424    HKL = GenHLaue(dvec,SGData,A)
1425    SSdH = [vec*h for h in range(-maxH,maxH+1)]
1426    SSdH = dict(zip(range(-maxH,maxH+1),SSdH))
1427    for h,k,l,d in HKL:
1428        ext = G2spc.GenHKLf([h,k,l],SGData)[0]  #h,k,l must be integral values here
1429        if not ext and d >= dmin:
1430            HKLs.append([h,k,l,0,d])
1431        for dH in SSdH:
1432            if dH:
1433                DH = SSdH[dH]
1434                H = [h+DH[0],k+DH[1],l+DH[2]]
1435                d = 1./np.sqrt(calc_rDsq(H,A))
1436                if d >= dmin:
1437                    HKLM = np.array([h,k,l,dH])
1438                    if (G2spc.checkSSLaue([h,k,l,dH],SGData,SSGData) and G2spc.checkSSextc(HKLM,SSGData)) or ifMag:
1439                        HKLs.append([h,k,l,dH,d])
1440    return HKLs
1441
1442def LaueUnique2(SGData,refList):
1443    ''' Impose Laue symmetry on hkl
1444
1445    :param SGData: space group data from 'P '+Laue
1446    :param HKLF: np.array([[h,k,l,...]]) reflection set to be converted
1447
1448    :return: HKLF new reflection array with imposed Laue symmetry
1449    '''
1450    for ref in refList:
1451        H = ref[:3]
1452        Uniq = G2spc.GenHKLf(H,SGData)[2]
1453        Uniq = G2mth.sortArray(G2mth.sortArray(G2mth.sortArray(Uniq,2),1),0)
1454        ref[:3] = Uniq[-1]
1455    return refList
1456
1457def LaueUnique(Laue,HKLF):
1458    ''' Impose Laue symmetry on hkl
1459
1460    :param str Laue: Laue symbol, as below
1461
1462      centrosymmetric Laue groups::
1463
1464            ['-1','2/m','112/m','2/m11','mmm','-42m','-4m2','4/mmm','-3',
1465            '-31m','-3m1','6/m','6/mmm','m3','m3m']
1466
1467      noncentrosymmetric Laue groups::
1468
1469           ['1','2','211','112','m','m11','11m','222','mm2','m2m','2mm',
1470           '4','-4','422','4mm','3','312','321','31m','3m1','6','-6',
1471           '622','6mm','-62m','-6m2','23','432','-43m']
1472
1473    :param HKLF: np.array([[h,k,l,...]]) reflection set to be converted
1474
1475    :returns: HKLF new reflection array with imposed Laue symmetry
1476    '''
1477
1478    HKLFT = HKLF.T
1479    mat41 = np.array([[0,1,0],[-1,0,0],[0,0,1]])    #hkl -> k,-h,l
1480    mat43 = np.array([[0,-1,0],[1,0,0],[0,0,1]])    #hkl -> -k,h,l
1481    mat4bar = np.array([[0,-1,0],[1,0,0],[0,0,-1]]) #hkl -> k,-h,-l
1482    mat31 = np.array([[-1,-1,0],[1,0,0],[0,0,1]])   #hkl -> ihl = -h-k,h,l
1483    mat32 = np.array([[0,1,0],[-1,-1,0],[0,0,1]])   #hkl -> kil = k,-h-k,l
1484    matd3 = np.array([[0,1,0],[0,0,1],[1,0,0]])     #hkl -> k,l,h
1485    matd3q = np.array([[0,0,-1],[-1,0,0],[0,1,0]])  #hkl -> -l,-h,k
1486    matd3t = np.array([[0,0,-1],[1,0,0],[0,-1,0]])  #hkl -> -l,h,-k
1487    mat6 = np.array([[1,1,0],[-1,0,0],[0,0,1]])     #hkl -> h+k,-h,l really 65
1488    matdm = np.array([[0,1,0],[1,0,0],[0,0,1]])     #hkl -> k,h,l
1489    matdmp = np.array([[-1,-1,0],[0,1,0],[0,0,1]])  #hkl -> -h-k,k,l
1490    matkm = np.array([[-1,0,0],[1,1,0],[0,0,1]])    #hkl -> -h,h+k,l
1491    matd2 = np.array([[0,1,0],[1,0,0],[0,0,-1]])    #hkl -> k,h,-l
1492    matdm3 = np.array([[1,0,0],[0,0,1],[0,1,0]])    #hkl -> h,l,k
1493    mat2d43 = np.array([[0,1,0],[1,0,0],[0,0,1]])   #hkl -> k,-h,l
1494    matk2 = np.array([[-1,0,0],[1,1,0],[0,0,-1]])   #hkl -> -h,-i,-l
1495    #triclinic
1496    if Laue == '1': #ok
1497        pass
1498    elif Laue == '-1':  #ok
1499        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
1500        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]<0),HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
1501        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
1502    #monoclinic
1503    #noncentrosymmetric - all ok
1504    elif Laue == '2':
1505        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3])
1506        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3])
1507    elif Laue == '1 1 2':
1508        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1509        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]<0),HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1510    elif Laue == '2 1 1':
1511        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1512        HKLFT[:3] = np.where((HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1513    elif Laue == 'm':
1514        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1515    elif Laue == 'm 1 1':
1516        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1517    elif Laue == '1 1 m':
1518        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1519    #centrosymmetric - all ok
1520    elif Laue == '2/m 1 1':
1521        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1522        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1523        HKLFT[:3] = np.where((HKLFT[2]*HKLFT[0]==0)&(HKLFT[1]<0),HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1524    elif Laue == '2/m':
1525        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3])
1526        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1527        HKLFT[:3] = np.where((HKLFT[0]*HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1528    elif Laue == '1 1 2/m':
1529        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1530        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1531        HKLFT[:3] = np.where((HKLFT[1]*HKLFT[2]==0)&(HKLFT[0]<0),HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1532    #orthorhombic
1533    #noncentrosymmetric - all OK
1534    elif Laue == '2 2 2':
1535        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1536        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1537        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1538        HKLFT[:3] = np.where((HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1539    elif Laue == 'm m 2':
1540        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1541        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1542    elif Laue == '2 m m':
1543        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1544        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1545    elif Laue == 'm 2 m':
1546        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
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 == 'm m m':
1550        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1551        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1552        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1553    #tetragonal
1554    #noncentrosymmetric - all ok
1555    elif Laue == '4':
1556        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1557        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1558        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]>0),np.squeeze(np.inner(HKLF[:,:3],mat41[nxs,:,:])).T,HKLFT[:3])
1559    elif Laue == '-4':
1560        HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1561        HKLFT[:3] = np.where(HKLFT[0]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1562        HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1563        HKLFT[:3] = np.where(HKLFT[1]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1564        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1565    elif Laue == '4 2 2':
1566        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1567        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1568        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1569        HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]<HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1570        HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])   #in lieu od 2-fold
1571    elif Laue == '4 m m':
1572        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1573        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1574        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1575        HKLFT[:3] = np.where(HKLFT[0]<HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1576    elif Laue == '-4 2 m':
1577        HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1578        HKLFT[:3] = np.where(HKLFT[0]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
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,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1581        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1582        HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1583        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1584    elif Laue == '-4 m 2':
1585        HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1586        HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1587        HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]<=0),np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1588        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]<0),HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1589        HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]==0),np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1590        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1591        HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[0]>HKLFT[1]),np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1592    #centrosymmetric - all ok
1593    elif Laue == '4/m':
1594        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1595        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1596        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1597        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]>0),np.squeeze(np.inner(HKLF[:,:3],mat41[nxs,:,:])).T,HKLFT[:3])
1598    elif Laue == '4/m m m':
1599        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1600        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1601        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1602        HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],mat41[nxs,:,:])).T,HKLFT[:3])
1603        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1604    #trigonal - all hex cell
1605    #noncentrosymmetric - all ok
1606    elif Laue == '3':
1607        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1608        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1609        HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1610    elif Laue == '3 1 2':
1611        HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matk2[nxs,:,:])).T,HKLFT[:3])
1612        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1613        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1614        HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1615        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],matk2[nxs,:,:])).T,HKLFT[:3])
1616    elif Laue == '3 2 1':
1617        HKLFT[:3] = np.where(HKLFT[0]<=-2*HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1618        HKLFT[:3] = np.where(HKLFT[1]<-2*HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1619        HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1620        HKLFT[:3] = np.where((HKLFT[2]>0)&(HKLFT[1]==HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1621        HKLFT[:3] = np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T
1622        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])
1623    elif Laue == '3 1 m':
1624        HKLFT[:3] = np.where(HKLFT[0]>=HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1625        HKLFT[:3] = np.where(2*HKLFT[1]<-HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1626        HKLFT[:3] = np.where(HKLFT[1]>-2*HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matdmp[nxs,:,:])).T,HKLFT[:3])
1627        HKLFT[:3] = np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T
1628    elif Laue == '3 m 1':
1629        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1630        HKLFT[:3] = np.where((HKLFT[1]+HKLFT[0])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1631        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],matkm[nxs,:,:])).T,HKLFT[:3])
1632    #centrosymmetric
1633    elif Laue == '-3':  #ok
1634        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
1635        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1636        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1637        HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1638        HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[0]<0),-np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1639        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],-mat31[nxs,:,:])).T,HKLFT[:3])
1640    elif Laue == '-3 m 1':  #ok
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[1]+HKLFT[0])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1643        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],matkm[nxs,:,:])).T,HKLFT[:3])
1644        HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1645        HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]<HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1646    elif Laue == '-3 1 m':  #ok
1647        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
1648        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1649        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1650        HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1651        HKLFT[:3] = np.where(HKLFT[0]<=0,np.squeeze(np.inner(HKLF[:,:3],-mat31[nxs,:,:])).T,HKLFT[:3])
1652        HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1653    #hexagonal
1654    #noncentrosymmetric
1655    elif Laue == '6':   #ok
1656        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1657        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1658        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1659        HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1660    elif Laue == '-6':  #ok
1661        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1662        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1663        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1664        HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1665    elif Laue == '6 2 2':   #ok
1666        HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1667        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1668        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1669        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1670        HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1671        HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[0]>HKLFT[1]),np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1672    elif Laue == '6 m m':   #ok
1673        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1674        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1675        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1676        HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1677        HKLFT[:3] = np.where(HKLFT[0]>HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1678    elif Laue == '-6 m 2':  #ok
1679        HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matk2[nxs,:,:])).T,HKLFT[:3])
1680        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1681        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1682        HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1683        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],matk2[nxs,:,:])).T,HKLFT[:3])
1684        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1685    elif Laue == '-6 2 m':  #ok
1686        HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1687        HKLFT[:3] = np.where(HKLFT[0]<=-2*HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1688        HKLFT[:3] = np.where(HKLFT[1]<-2*HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1689        HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1690        HKLFT[:3] = np.where((HKLFT[2]>0)&(HKLFT[1]==HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1691        HKLFT[:3] = np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T
1692        HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1693        HKLFT[:3] = np.where(HKLFT[0]>HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1694    #centrosymmetric
1695    elif Laue == '6/m': #ok
1696        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1697        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1698        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1699        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1700        HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1701    elif Laue == '6/m m m': #ok
1702        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1703        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1704        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[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.T[nxs,:,:])).T,HKLFT[:3])
1707    #cubic - all ok
1708    #noncentrosymmetric -
1709    elif Laue == '2 3':
1710        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1711        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1712        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1713        HKLFT[:3] = np.where((HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1714        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])
1715        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])
1716        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])
1717        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])
1718        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3])
1719    elif Laue == '4 3 2':
1720        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1721        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1722        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1723        HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]<HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1724        HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])   #in lieu od 2-fold
1725        HKLFT[:3] = np.where((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2]),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
1726        HKLFT[:3] = np.where((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2]),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
1727        HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat2d43[nxs,:,:])).T,HKLFT[:3])
1728    elif Laue == '-4 3 m':
1729        HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1730        HKLFT[:3] = np.where(HKLFT[0]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1731        HKLFT[:3] = np.where(HKLFT[0]<=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],mat4bar[nxs,:,:])).T,HKLFT[:3])
1733        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1734        HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1735        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1736        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])
1737        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])
1738        HKLFT[:3] = np.where((HKLFT[2]>=0)&(HKLFT[1]<HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1739        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3])
1740        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])
1741        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])
1742    #centrosymmetric
1743    elif Laue == 'm 3':
1744        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1745        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1746        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1747        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])
1748        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])
1749    elif Laue == 'm 3 m':
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,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1752        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1753        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])
1754        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])
1755        HKLFT[:3] = np.where(HKLFT[0]>HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1756    return HKLFT.T
1757
1758
1759#Spherical harmonics routines
1760def OdfChk(SGLaue,L,M):
1761    'needs doc string'
1762    if not L%2 and abs(M) <= L:
1763        if SGLaue == '0':                      #cylindrical symmetry
1764            if M == 0: return True
1765        elif SGLaue == '-1':
1766            return True
1767        elif SGLaue == '2/m':
1768            if not abs(M)%2: return True
1769        elif SGLaue == 'mmm':
1770            if not abs(M)%2 and M >= 0: return True
1771        elif SGLaue == '4/m':
1772            if not abs(M)%4: return True
1773        elif SGLaue == '4/mmm':
1774            if not abs(M)%4 and M >= 0: return True
1775        elif SGLaue in ['3R','3']:
1776            if not abs(M)%3: return True
1777        elif SGLaue in ['3mR','3m1','31m']:
1778            if not abs(M)%3 and M >= 0: return True
1779        elif SGLaue == '6/m':
1780            if not abs(M)%6: return True
1781        elif SGLaue == '6/mmm':
1782            if not abs(M)%6 and M >= 0: return True
1783        elif SGLaue == 'm3':
1784            if M > 0:
1785                if L%12 == 2:
1786                    if M <= L//12: return True
1787                else:
1788                    if M <= L//12+1: return True
1789        elif SGLaue == 'm3m':
1790            if M > 0:
1791                if L%12 == 2:
1792                    if M <= L//12: return True
1793                else:
1794                    if M <= L//12+1: return True
1795    return False
1796
1797def GenSHCoeff(SGLaue,SamSym,L,IfLMN=True):
1798    'needs doc string'
1799    coeffNames = []
1800    for iord in [2*i+2 for i in range(L//2)]:
1801        for m in [i-iord for i in range(2*iord+1)]:
1802            if OdfChk(SamSym,iord,m):
1803                for n in [i-iord for i in range(2*iord+1)]:
1804                    if OdfChk(SGLaue,iord,n):
1805                        if IfLMN:
1806                            coeffNames.append('C(%d,%d,%d)'%(iord,m,n))
1807                        else:
1808                            coeffNames.append('C(%d,%d)'%(iord,n))
1809    return coeffNames
1810
1811def CrsAng(H,cell,SGData):
1812    'needs doc string'
1813    a,b,c,al,be,ga = cell
1814    SQ3 = 1.732050807569
1815    H1 = np.array([1,0,0])
1816    H2 = np.array([0,1,0])
1817    H3 = np.array([0,0,1])
1818    H4 = np.array([1,1,1])
1819    G,g = cell2Gmat(cell)
1820    Laue = SGData['SGLaue']
1821    Naxis = SGData['SGUniq']
1822    if len(H.shape) == 1:
1823        DH = np.inner(H,np.inner(G,H))
1824    else:
1825        DH = np.array([np.inner(h,np.inner(G,h)) for h in H])
1826    if Laue == '2/m':
1827        if Naxis == 'a':
1828            DR = np.inner(H1,np.inner(G,H1))
1829            DHR = np.inner(H,np.inner(G,H1))
1830        elif Naxis == 'b':
1831            DR = np.inner(H2,np.inner(G,H2))
1832            DHR = np.inner(H,np.inner(G,H2))
1833        else:
1834            DR = np.inner(H3,np.inner(G,H3))
1835            DHR = np.inner(H,np.inner(G,H3))
1836    elif Laue in ['R3','R3m']:
1837        DR = np.inner(H4,np.inner(G,H4))
1838        DHR = np.inner(H,np.inner(G,H4))
1839    else:
1840        DR = np.inner(H3,np.inner(G,H3))
1841        DHR = np.inner(H,np.inner(G,H3))
1842    DHR /= np.sqrt(DR*DH)
1843    phi = np.where(DHR <= 1.0,acosd(DHR),0.0)
1844    if Laue == '-1':
1845        BA = H.T[1]*a/(b-H.T[0]*cosd(ga))
1846        BB = H.T[0]*sind(ga)**2
1847    elif Laue == '2/m':
1848        if Naxis == 'a':
1849            BA = H.T[2]*b/(c-H.T[1]*cosd(al))
1850            BB = H.T[1]*sind(al)**2
1851        elif Naxis == 'b':
1852            BA = H.T[0]*c/(a-H.T[2]*cosd(be))
1853            BB = H.T[2]*sind(be)**2
1854        else:
1855            BA = H.T[1]*a/(b-H.T[0]*cosd(ga))
1856            BB = H.T[0]*sind(ga)**2
1857    elif Laue in ['mmm','4/m','4/mmm']:
1858        BA = H.T[1]*a
1859        BB = H.T[0]*b
1860    elif Laue in ['3R','3mR']:
1861        BA = H.T[0]+H.T[1]-2.0*H.T[2]
1862        BB = SQ3*(H.T[0]-H.T[1])
1863    elif Laue in ['m3','m3m']:
1864        BA = H.T[1]
1865        BB = H.T[0]
1866    else:
1867        BA = H.T[0]+2.0*H.T[1]
1868        BB = SQ3*H.T[0]
1869    beta = atan2d(BA,BB)
1870    return phi,beta
1871
1872def SamAng(Tth,Gangls,Sangl,IFCoup):
1873    """Compute sample orientation angles vs laboratory coord. system
1874
1875    :param Tth:        Signed theta
1876    :param Gangls:     Sample goniometer angles phi,chi,omega,azmuth
1877    :param Sangl:      Sample angle zeros om-0, chi-0, phi-0
1878    :param IFCoup:     True if omega & 2-theta coupled in CW scan
1879    :returns:
1880        psi,gam:    Sample odf angles
1881        dPSdA,dGMdA:    Angle zero derivatives
1882    """
1883
1884    if IFCoup:
1885        GSomeg = sind(Gangls[2]+Tth)
1886        GComeg = cosd(Gangls[2]+Tth)
1887    else:
1888        GSomeg = sind(Gangls[2])
1889        GComeg = cosd(Gangls[2])
1890    GSTth = sind(Tth)
1891    GCTth = cosd(Tth)
1892    GSazm = sind(Gangls[3])
1893    GCazm = cosd(Gangls[3])
1894    GSchi = sind(Gangls[1])
1895    GCchi = cosd(Gangls[1])
1896    GSphi = sind(Gangls[0]+Sangl[2])
1897    GCphi = cosd(Gangls[0]+Sangl[2])
1898    SSomeg = sind(Sangl[0])
1899    SComeg = cosd(Sangl[0])
1900    SSchi = sind(Sangl[1])
1901    SCchi = cosd(Sangl[1])
1902    AT = -GSTth*GComeg+GCTth*GCazm*GSomeg
1903    BT = GSTth*GSomeg+GCTth*GCazm*GComeg
1904    CT = -GCTth*GSazm*GSchi
1905    DT = -GCTth*GSazm*GCchi
1906
1907    BC1 = -AT*GSphi+(CT+BT*GCchi)*GCphi
1908    BC2 = DT-BT*GSchi
1909    BC3 = AT*GCphi+(CT+BT*GCchi)*GSphi
1910
1911    BC = BC1*SComeg*SCchi+BC2*SComeg*SSchi-BC3*SSomeg
1912    psi = acosd(BC)
1913
1914    BD = 1.0-BC**2
1915    C = np.where(BD>1.e-6,rpd/np.sqrt(BD),0.)
1916    dPSdA = [-C*(-BC1*SSomeg*SCchi-BC2*SSomeg*SSchi-BC3*SComeg),
1917        -C*(-BC1*SComeg*SSchi+BC2*SComeg*SCchi),
1918        -C*(-BC1*SSomeg-BC3*SComeg*SCchi)]
1919
1920    BA = -BC1*SSchi+BC2*SCchi
1921    BB = BC1*SSomeg*SCchi+BC2*SSomeg*SSchi+BC3*SComeg
1922    gam = atan2d(BB,BA)
1923
1924    BD = (BA**2+BB**2)/rpd
1925
1929
1930    dBBdO = BC1*SComeg*SCchi+BC2*SComeg*SSchi-BC3*SSomeg
1931    dBBdC = -BC1*SSomeg*SSchi+BC2*SSomeg*SCchi
1932    dBBdF = BC1*SComeg-BC3*SSomeg*SCchi
1933
1936
1937    return psi,gam,dPSdA,dGMdA
1938
1939BOH = {
1940'L=2':[[],[],[]],
1941'L=4':[[0.30469720,0.36418281],[],[]],
1942'L=6':[[-0.14104740,0.52775103],[],[]],
1943'L=8':[[0.28646862,0.21545346,0.32826995],[],[]],
1944'L=10':[[-0.16413497,0.33078546,0.39371345],[],[]],
1945'L=12':[[0.26141975,0.27266871,0.03277460,0.32589402],
1946    [0.09298802,-0.23773812,0.49446631,0.0],[]],
1947'L=14':[[-0.17557309,0.25821932,0.27709173,0.33645360],[],[]],
1948'L=16':[[0.24370673,0.29873515,0.06447688,0.00377,0.32574495],
1949    [0.12039646,-0.25330128,0.23950998,0.40962508,0.0],[]],
1950'L=18':[[-0.16914245,0.17017340,0.34598142,0.07433932,0.32696037],
1951    [-0.06901768,0.16006562,-0.24743528,0.47110273,0.0],[]],
1952'L=20':[[0.23067026,0.31151832,0.09287682,0.01089683,0.00037564,0.32573563],
1953    [0.13615420,-0.25048007,0.12882081,0.28642879,0.34620433,0.0],[]],
1954'L=22':[[-0.16109560,0.10244188,0.36285175,0.13377513,0.01314399,0.32585583],
1955    [-0.09620055,0.20244115,-0.22389483,0.17928946,0.42017231,0.0],[]],
1956'L=24':[[0.22050742,0.31770654,0.11661736,0.02049853,0.00150861,0.00003426,0.32573505],
1957    [0.13651722,-0.21386648,0.00522051,0.33939435,0.10837396,0.32914497,0.0],
1958    [0.05378596,-0.11945819,0.16272298,-0.26449730,0.44923956,0.0,0.0]],
1959'L=26':[[-0.15435003,0.05261630,0.35524646,0.18578869,0.03259103,0.00186197,0.32574594],
1960    [-0.11306511,0.22072681,-0.18706142,0.05439948,0.28122966,0.35634355,0.0],[]],
1961'L=28':[[0.21225019,0.32031716,0.13604702,0.03132468,0.00362703,0.00018294,0.00000294,0.32573501],
1962    [0.13219496,-0.17206256,-0.08742608,0.32671661,0.17973107,0.02567515,0.32619598,0.0],
1963    [0.07989184,-0.16735346,0.18839770,-0.20705337,0.12926808,0.42715602,0.0,0.0]],
1964'L=30':[[-0.14878368,0.01524973,0.33628434,0.22632587,0.05790047,0.00609812,0.00022898,0.32573594],
1965    [-0.11721726,0.20915005,-0.11723436,-0.07815329,0.31318947,0.13655742,0.33241385,0.0],
1966    [-0.04297703,0.09317876,-0.11831248,0.17355132,-0.28164031,0.42719361,0.0,0.0]],
1967'L=32':[[0.20533892,0.32087437,0.15187897,0.04249238,0.00670516,0.00054977,0.00002018,0.00000024,0.32573501],
1968    [0.12775091,-0.13523423,-0.14935701,0.28227378,0.23670434,0.05661270,0.00469819,0.32578978,0.0],
1969    [0.09703829,-0.19373733,0.18610682,-0.14407046,0.00220535,0.26897090,0.36633402,0.0,0.0]],
1970'L=34':[[-0.14409234,-0.01343681,0.31248977,0.25557722,0.08571889,0.01351208,0.00095792,0.00002550,0.32573508],
1971    [-0.11527834,0.18472133,-0.04403280,-0.16908618,0.27227021,0.21086614,0.04041752,0.32688152,0.0],
1972    [-0.06773139,0.14120811,-0.15835721,0.18357456,-0.19364673,0.08377174,0.43116318,0.0,0.0]]
1973}
1974
1975Lnorm = lambda L: 4.*np.pi/(2.0*L+1.)
1976
1977def GetKcl(L,N,SGLaue,phi,beta):
1978    'needs doc string'
1979    import pytexture as ptx
1980    if SGLaue in ['m3','m3m']:
1981        if 'array' in str(type(phi)) and np.any(phi.shape):
1982            Kcl = np.zeros_like(phi)
1983        else:
1984            Kcl = 0.
1985        for j in range(0,L+1,4):
1986            im = j//4
1987            if 'array' in str(type(phi)) and np.any(phi.shape):
1988                pcrs = ptx.pyplmpsi(L,j,len(phi),phi)[0]
1989            else:
1990                pcrs = ptx.pyplmpsi(L,j,1,phi)[0]
1991            Kcl += BOH['L=%d'%(L)][N-1][im]*pcrs*cosd(j*beta)
1992    else:
1993        if 'array' in str(type(phi)) and np.any(phi.shape):
1994            pcrs = ptx.pyplmpsi(L,N,len(phi),phi)[0]
1995        else:
1996            pcrs = ptx.pyplmpsi(L,N,1,phi)[0]
1997        pcrs *= RSQ2PI
1998        if N:
1999            pcrs *= SQ2
2000        if SGLaue in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']:
2001            if SGLaue in ['3mR','3m1','31m']:
2002                if N%6 == 3:
2003                    Kcl = pcrs*sind(N*beta)
2004                else:
2005                    Kcl = pcrs*cosd(N*beta)
2006            else:
2007                Kcl = pcrs*cosd(N*beta)
2008        else:
2009            Kcl = pcrs*(cosd(N*beta)+sind(N*beta))
2010    return Kcl
2011
2012def GetKsl(L,M,SamSym,psi,gam):
2013    'needs doc string'
2014    import pytexture as ptx
2015    if 'array' in str(type(psi)) and np.any(psi.shape):
2016        psrs,dpdps = ptx.pyplmpsi(L,M,len(psi),psi)
2017    else:
2018        psrs,dpdps = ptx.pyplmpsi(L,M,1,psi)
2019    psrs *= RSQ2PI
2020    dpdps *= RSQ2PI
2021    if M:
2022        psrs *= SQ2
2023        dpdps *= SQ2
2024    if SamSym in ['mmm',]:
2025        dum = cosd(M*gam)
2026        Ksl = psrs*dum
2027        dKsdp = dpdps*dum
2028        dKsdg = -psrs*M*sind(M*gam)
2029    else:
2030        dum = cosd(M*gam)+sind(M*gam)
2031        Ksl = psrs*dum
2032        dKsdp = dpdps*dum
2033        dKsdg = psrs*M*(-sind(M*gam)+cosd(M*gam))
2034    return Ksl,dKsdp,dKsdg
2035
2036def GetKclKsl(L,N,SGLaue,psi,phi,beta):
2037    """
2038    This is used for spherical harmonics description of preferred orientation;
2039        cylindrical symmetry only (M=0) and no sample angle derivatives returned
2040    """
2041    import pytexture as ptx
2042    Ksl,x = ptx.pyplmpsi(L,0,1,psi)
2043    Ksl *= RSQ2PI
2044    if SGLaue in ['m3','m3m']:
2045        Kcl = 0.0
2046        for j in range(0,L+1,4):
2047            im = j//4
2048            pcrs,dum = ptx.pyplmpsi(L,j,1,phi)
2049            Kcl += BOH['L=%d'%(L)][N-1][im]*pcrs*cosd(j*beta)
2050    else:
2051        pcrs,dum = ptx.pyplmpsi(L,N,1,phi)
2052        pcrs *= RSQ2PI
2053        if N:
2054            pcrs *= SQ2
2055        if SGLaue in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']:
2056            if SGLaue in ['3mR','3m1','31m']:
2057                if N%6 == 3:
2058                    Kcl = pcrs*sind(N*beta)
2059                else:
2060                    Kcl = pcrs*cosd(N*beta)
2061            else:
2062                Kcl = pcrs*cosd(N*beta)
2063        else:
2064            Kcl = pcrs*(cosd(N*beta)+sind(N*beta))
2065    return Kcl*Ksl,Lnorm(L)
2066
2067def Glnh(Start,SHCoef,psi,gam,SamSym):
2068    'needs doc string'
2069    import pytexture as ptx
2070
2071    if Start:
2072        ptx.pyqlmninit()
2073        Start = False
2074    Fln = np.zeros(len(SHCoef))
2075    for i,term in enumerate(SHCoef):
2076        l,m,n = eval(term.strip('C'))
2077        pcrs,dum = ptx.pyplmpsi(l,m,1,psi)
2078        pcrs *= RSQPI
2079        if m == 0:
2080            pcrs /= SQ2
2081        if SamSym in ['mmm',]:
2082            Ksl = pcrs*cosd(m*gam)
2083        else:
2084            Ksl = pcrs*(cosd(m*gam)+sind(m*gam))
2085        Fln[i] = SHCoef[term]*Ksl*Lnorm(l)
2086    ODFln = dict(zip(SHCoef.keys(),list(zip(SHCoef.values(),Fln))))
2087    return ODFln
2088
2089def Flnh(Start,SHCoef,phi,beta,SGData):
2090    'needs doc string'
2091    import pytexture as ptx
2092
2093    if Start:
2094        ptx.pyqlmninit()
2095        Start = False
2096    Fln = np.zeros(len(SHCoef))
2097    for i,term in enumerate(SHCoef):
2098        l,m,n = eval(term.strip('C'))
2099        if SGData['SGLaue'] in ['m3','m3m']:
2100            Kcl = 0.0
2101            for j in range(0,l+1,4):
2102                im = j//4
2103                pcrs,dum = ptx.pyplmpsi(l,j,1,phi)
2104                Kcl += BOH['L='+str(l)][n-1][im]*pcrs*cosd(j*beta)
2105        else:                #all but cubic
2106            pcrs,dum = ptx.pyplmpsi(l,n,1,phi)
2107            pcrs *= RSQPI
2108            if n == 0:
2109                pcrs /= SQ2
2110            if SGData['SGLaue'] in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']:
2111               if SGData['SGLaue'] in ['3mR','3m1','31m']:
2112                   if n%6 == 3:
2113                       Kcl = pcrs*sind(n*beta)
2114                   else:
2115                       Kcl = pcrs*cosd(n*beta)
2116               else:
2117                   Kcl = pcrs*cosd(n*beta)
2118            else:
2119                Kcl = pcrs*(cosd(n*beta)+sind(n*beta))
2120        Fln[i] = SHCoef[term]*Kcl*Lnorm(l)
2121    ODFln = dict(zip(SHCoef.keys(),list(zip(SHCoef.values(),Fln))))
2122    return ODFln
2123
2124def polfcal(ODFln,SamSym,psi,gam):
2125    '''Perform a pole figure computation.
2126    Note that the the number of gam values must either be 1 or must
2127    match psi. Updated for numpy 1.8.0
2128    '''
2129    import pytexture as ptx
2130    PolVal = np.ones_like(psi)
2131    for term in ODFln:
2132        if abs(ODFln[term][1]) > 1.e-3:
2133            l,m,n = eval(term.strip('C'))
2134            psrs,dum = ptx.pyplmpsi(l,m,len(psi),psi)
2135            if SamSym in ['-1','2/m']:
2136                if m:
2137                    Ksl = RSQPI*psrs*(cosd(m*gam)+sind(m*gam))
2138                else:
2139                    Ksl = RSQPI*psrs/SQ2
2140            else:
2141                if m:
2142                    Ksl = RSQPI*psrs*cosd(m*gam)
2143                else:
2144                    Ksl = RSQPI*psrs/SQ2
2145            PolVal += ODFln[term][1]*Ksl
2146    return PolVal
2147
2148def invpolfcal(ODFln,SGData,phi,beta):
2149    'needs doc string'
2150    import pytexture as ptx
2151
2152    invPolVal = np.ones_like(beta)
2153    for term in ODFln:
2154        if abs(ODFln[term][1]) > 1.e-3:
2155            l,m,n = eval(term.strip('C'))
2156            if SGData['SGLaue'] in ['m3','m3m']:
2157                Kcl = 0.0
2158                for j in range(0,l+1,4):
2159                    im = j//4
2160                    pcrs,dum = ptx.pyplmpsi(l,j,len(beta),phi)
2161                    Kcl += BOH['L=%d'%(l)][n-1][im]*pcrs*cosd(j*beta)
2162            else:                #all but cubic
2163                pcrs,dum = ptx.pyplmpsi(l,n,len(beta),phi)
2164                pcrs *= RSQPI
2165                if n == 0:
2166                    pcrs /= SQ2
2167                if SGData['SGLaue'] in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']:
2168                   if SGData['SGLaue'] in ['3mR','3m1','31m']:
2169                       if n%6 == 3:
2170                           Kcl = pcrs*sind(n*beta)
2171                       else:
2172                           Kcl = pcrs*cosd(n*beta)
2173                   else:
2174                       Kcl = pcrs*cosd(n*beta)
2175                else:
2176                    Kcl = pcrs*(cosd(n*beta)+sind(n*beta))
2177            invPolVal += ODFln[term][1]*Kcl
2178    return invPolVal
2179
2180
2181def textureIndex(SHCoef):
2182    'needs doc string'
2183    Tindx = 1.0
2184    for term in SHCoef:
2185        l = eval(term.strip('C'))[0]
2186        Tindx += SHCoef[term]**2/(2.0*l+1.)
2187    return Tindx
2188
2189# self-test materials follow.
2190selftestlist = []
2191'''Defines a list of self-tests'''
2192selftestquiet = True
2193def _ReportTest():
2194    'Report name and doc string of current routine when selftestquiet is False'
2195    if not selftestquiet:
2196        import inspect
2197        caller = inspect.stack()[1][3]
2198        doc = eval(caller).__doc__
2199        if doc is not None:
2200            print('testing '+__file__+' with '+caller+' ('+doc+')')
2201        else:
2202            print('testing '+__file__()+" with "+caller)
2203NeedTestData = True
2204def TestData():
2205    array = np.array
2206    global NeedTestData
2207    NeedTestData = False
2208    global CellTestData
2209    # output from uctbx computed on platform darwin on 2010-05-28
2210    CellTestData = [
2211# cell, g, G, cell*, V, V*
2212  [(4, 4, 4, 90, 90, 90),
2213   array([[  1.60000000e+01,   9.79717439e-16,   9.79717439e-16],
2214       [  9.79717439e-16,   1.60000000e+01,   9.79717439e-16],
2215       [  9.79717439e-16,   9.79717439e-16,   1.60000000e+01]]), array([[  6.25000000e-02,   3.82702125e-18,   3.82702125e-18],
2216       [  3.82702125e-18,   6.25000000e-02,   3.82702125e-18],
2217       [  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],
2218# cell, g, G, cell*, V, V*
2219  [(4.0999999999999996, 5.2000000000000002, 6.2999999999999998, 100, 80, 130),
2220   array([[ 16.81      , -13.70423184,   4.48533243],
2221       [-13.70423184,  27.04      ,  -5.6887143 ],
2222       [  4.48533243,  -5.6887143 ,  39.69      ]]), array([[ 0.10206349,  0.05083339, -0.00424823],
2223       [ 0.05083339,  0.06344997,  0.00334956],
2224       [-0.00424823,  0.00334956,  0.02615544]]), (0.31947376387537696, 0.25189277536327803, 0.16172643497798223, 85.283666420376008, 94.716333579624006, 50.825714168082683), 100.98576357983838, 0.0099023858863968445],
2225# cell, g, G, cell*, V, V*
2226  [(3.5, 3.5, 6, 90, 90, 120),
2227   array([[  1.22500000e+01,  -6.12500000e+00,   1.28587914e-15],
2228       [ -6.12500000e+00,   1.22500000e+01,   1.28587914e-15],
2229       [  1.28587914e-15,   1.28587914e-15,   3.60000000e+01]]), array([[  1.08843537e-01,   5.44217687e-02,   3.36690552e-18],
2230       [  5.44217687e-02,   1.08843537e-01,   3.36690552e-18],
2231       [  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],
2232  ]
2233    global CoordTestData
2234    CoordTestData = [
2235# cell, ((frac, ortho),...)
2236  ((4,4,4,90,90,90,), [
2237 ((0.10000000000000001, 0.0, 0.0),(0.40000000000000002, 0.0, 0.0)),
2238 ((0.0, 0.10000000000000001, 0.0),(2.4492935982947065e-17, 0.40000000000000002, 0.0)),
2239 ((0.0, 0.0, 0.10000000000000001),(2.4492935982947065e-17, -2.4492935982947065e-17, 0.40000000000000002)),
2240 ((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(0.40000000000000013, 0.79999999999999993, 1.2)),
2241 ((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(0.80000000000000016, 1.2, 0.40000000000000002)),
2242 ((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(1.2, 0.80000000000000004, 0.40000000000000002)),
2243 ((0.5, 0.5, 0.5),(2.0, 1.9999999999999998, 2.0)),
2244]),
2245# cell, ((frac, ortho),...)
2246  ((4.1,5.2,6.3,100,80,130,), [
2247 ((0.10000000000000001, 0.0, 0.0),(0.40999999999999998, 0.0, 0.0)),
2248 ((0.0, 0.10000000000000001, 0.0),(-0.33424955703700043, 0.39834311042186865, 0.0)),
2249 ((0.0, 0.0, 0.10000000000000001),(0.10939835193016617, -0.051013289294572106, 0.6183281045774256)),
2250 ((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(0.069695941716497567, 0.64364635296002093, 1.8549843137322766)),
2251 ((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(-0.073350319180835066, 1.1440160419710339, 0.6183281045774256)),
2252 ((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(0.67089923785616512, 0.74567293154916525, 0.6183281045774256)),
2253 ((0.5, 0.5, 0.5),(0.92574397446582857, 1.7366491056364828, 3.0916405228871278)),
2254]),
2255# cell, ((frac, ortho),...)
2256  ((3.5,3.5,6,90,90,120,), [
2257 ((0.10000000000000001, 0.0, 0.0),(0.35000000000000003, 0.0, 0.0)),
2258 ((0.0, 0.10000000000000001, 0.0),(-0.17499999999999993, 0.3031088913245536, 0.0)),
2259 ((0.0, 0.0, 0.10000000000000001),(3.6739403974420595e-17, -3.6739403974420595e-17, 0.60000000000000009)),
2260 ((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(2.7675166561703527e-16, 0.60621778264910708, 1.7999999999999998)),
2261 ((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(0.17500000000000041, 0.90932667397366063, 0.60000000000000009)),
2262 ((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(0.70000000000000018, 0.6062177826491072, 0.60000000000000009)),
2263 ((0.5, 0.5, 0.5),(0.87500000000000067, 1.5155444566227676, 3.0)),
2264]),
2265]
2266    global LaueTestData             #generated by GSAS
2267    LaueTestData = {
2268    '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),
2269        (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),
2270        (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))],
2271    '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),
2272        (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),
2273        (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),
2274        (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))],
2275    '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),
2276        (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),
2277        (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),
2278        (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),
2279        (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),
2280        (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),
2281        (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),
2282        (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),
2283        (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),
2284        (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),
2285        (2,1,5,6),(2,1,-5,6),(3,-1,-5,6))],
2286    '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),
2287        (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),
2288        (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),
2289        (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),
2290        (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),
2291        (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),
2292        (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),
2293        (3,0,4,6),(3,1,-2,12),(3,0,-4,6),(1,1,6,12),(2,2,3,12))],
2294    '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),
2295        (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),
2296        (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),
2297        (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),
2298        (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),
2299        (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),
2300        (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))],
2301    }
2302
2303    global FLnhTestData
2304    FLnhTestData = [{
2305    'C(4,0,0)': (0.965, 0.42760447),
2306    'C(2,0,0)': (1.0122, -0.80233610),
2307    'C(2,0,2)': (0.0061, 8.37491546E-03),
2308    'C(6,0,4)': (-0.0898, 4.37985696E-02),
2309    'C(6,0,6)': (-0.1369, -9.04081762E-02),
2310    'C(6,0,0)': (0.5935, -0.18234928),
2311    'C(4,0,4)': (0.1872, 0.16358127),
2312    'C(6,0,2)': (0.6193, 0.27573633),
2313    'C(4,0,2)': (-0.1897, 0.12530720)},[1,0,0]]
2314def test0():
2315    if NeedTestData: TestData()
2316    msg = 'test cell2Gmat, fillgmat, Gmat2cell'
2317    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2318        G, g = cell2Gmat(cell)
2319        assert np.allclose(G,tG),msg
2320        assert np.allclose(g,tg),msg
2321        tcell = Gmat2cell(g)
2322        assert np.allclose(cell,tcell),msg
2323        tcell = Gmat2cell(G)
2324        assert np.allclose(tcell,trcell),msg
2325if __name__ == '__main__': selftestlist.append(test0)
2326
2327def test1():
2328    'test cell2A and A2Gmat'
2329    _ReportTest()
2330    if NeedTestData: TestData()
2331    msg = 'test cell2A and A2Gmat'
2332    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2333        G, g = A2Gmat(cell2A(cell))
2334        assert np.allclose(G,tG),msg
2335        assert np.allclose(g,tg),msg
2336if __name__ == '__main__': selftestlist.append(test1)
2337
2338def test2():
2339    'test Gmat2A, A2cell, A2Gmat, Gmat2cell'
2340    _ReportTest()
2341    if NeedTestData: TestData()
2342    msg = 'test Gmat2A, A2cell, A2Gmat, Gmat2cell'
2343    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2344        G, g = cell2Gmat(cell)
2345        tcell = A2cell(Gmat2A(G))
2346        assert np.allclose(cell,tcell),msg
2347if __name__ == '__main__': selftestlist.append(test2)
2348
2349def test3():
2350    'test invcell2Gmat'
2351    _ReportTest()
2352    if NeedTestData: TestData()
2353    msg = 'test invcell2Gmat'
2354    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2355        G, g = invcell2Gmat(trcell)
2356        assert np.allclose(G,tG),msg
2357        assert np.allclose(g,tg),msg
2358if __name__ == '__main__': selftestlist.append(test3)
2359
2360def test4():
2361    'test calc_rVsq, calc_rV, calc_V'
2362    _ReportTest()
2363    if NeedTestData: TestData()
2364    msg = 'test calc_rVsq, calc_rV, calc_V'
2365    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2366        assert np.allclose(calc_rV(cell2A(cell)),trV), msg
2367        assert np.allclose(calc_V(cell2A(cell)),tV), msg
2368if __name__ == '__main__': selftestlist.append(test4)
2369
2370def test5():
2371    'test A2invcell'
2372    _ReportTest()
2373    if NeedTestData: TestData()
2374    msg = 'test A2invcell'
2375    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2376        rcell = A2invcell(cell2A(cell))
2377        assert np.allclose(rcell,trcell),msg
2378if __name__ == '__main__': selftestlist.append(test5)
2379
2380def test6():
2381    'test cell2AB'
2382    _ReportTest()
2383    if NeedTestData: TestData()
2384    msg = 'test cell2AB'
2385    for (cell,coordlist) in CoordTestData:
2386        A,B = cell2AB(cell)
2387        for (frac,ortho) in coordlist:
2388            to = np.inner(A,frac)
2389            tf = np.inner(B,to)
2390            assert np.allclose(ortho,to), msg
2391            assert np.allclose(frac,tf), msg
2392            to = np.sum(A*frac,axis=1)
2393            tf = np.sum(B*to,axis=1)
2394            assert np.allclose(ortho,to), msg
2395            assert np.allclose(frac,tf), msg
2396if __name__ == '__main__': selftestlist.append(test6)
2397
2398def test7():
2399    'test GetBraviasNum(...) and GenHBravais(...)'
2400    _ReportTest()
2401    import os.path
2402    import sys
2403    import GSASIIspc as spc
2404    testdir = os.path.join(os.path.split(os.path.abspath( __file__ ))[0],'testinp')
2405    if os.path.exists(testdir):
2406        if testdir not in sys.path: sys.path.insert(0,testdir)
2407    import sgtbxlattinp
2408    derror = 1e-4
2409    def indexmatch(hklin, hkllist, system):
2410        for hklref in hkllist:
2411            hklref = list(hklref)
2412            # these permutations are far from complete, but are sufficient to
2413            # allow the test to complete
2414            if system == 'cubic':
2415                permlist = [(1,2,3),(1,3,2),(2,1,3),(2,3,1),(3,1,2),(3,2,1),]
2416            elif system == 'monoclinic':
2417                permlist = [(1,2,3),(-1,2,-3)]
2418            else:
2419                permlist = [(1,2,3)]
2420
2421            for perm in permlist:
2422                hkl = [abs(i) * hklin[abs(i)-1] / i for i in perm]
2423                if hkl == hklref: return True
2424                if [-i for i in hkl] == hklref: return True
2425        else:
2426            return False
2427
2428    for key in sgtbxlattinp.sgtbx7:
2429        spdict = spc.SpcGroup(key)
2430        cell = sgtbxlattinp.sgtbx7[key][0]
2431        system = spdict[1]['SGSys']
2432        center = spdict[1]['SGLatt']
2433
2434        bravcode = GetBraviasNum(center, system)
2435
2436        g2list = GenHBravais(sgtbxlattinp.dmin, bravcode, cell2A(cell))
2437
2438        assert len(sgtbxlattinp.sgtbx7[key][1]) == len(g2list), 'Reflection lists differ for %s' % key
2439        for h,k,l,d,num in g2list:
2440            for hkllist,dref in sgtbxlattinp.sgtbx7[key][1]:
2441                if abs(d-dref) < derror:
2442                    if indexmatch((h,k,l,), hkllist, system):
2443                        break
2444            else:
2445                assert 0,'No match for %s at %s (%s)' % ((h,k,l),d,key)
2446if __name__ == '__main__': selftestlist.append(test7)
2447
2448def test8():
2449    'test GenHLaue'
2450    _ReportTest()
2451    import GSASIIspc as spc
2452    import sgtbxlattinp
2453    derror = 1e-4
2454    dmin = sgtbxlattinp.dmin
2455
2456    def indexmatch(hklin, hklref, system, axis):
2457        # these permutations are far from complete, but are sufficient to
2458        # allow the test to complete
2459        if system == 'cubic':
2460            permlist = [(1,2,3),(1,3,2),(2,1,3),(2,3,1),(3,1,2),(3,2,1),]
2461        elif system == 'monoclinic' and axis=='b':
2462            permlist = [(1,2,3),(-1,2,-3)]
2463        elif system == 'monoclinic' and axis=='a':
2464            permlist = [(1,2,3),(1,-2,-3)]
2465        elif system == 'monoclinic' and axis=='c':
2466            permlist = [(1,2,3),(-1,-2,3)]
2467        elif system == 'trigonal':
2468            permlist = [(1,2,3),(2,1,3),(-1,-2,3),(-2,-1,3)]
2469        elif system == 'rhombohedral':
2470            permlist = [(1,2,3),(2,3,1),(3,1,2)]
2471        else:
2472            permlist = [(1,2,3)]
2473
2474        hklref = list(hklref)
2475        for perm in permlist:
2476            hkl = [abs(i) * hklin[abs(i)-1] / i for i in perm]
2477            if hkl == hklref: return True
2478            if [-i for i in hkl] == hklref: return True
2479        return False
2480
2481    for key in sgtbxlattinp.sgtbx8:
2482        spdict = spc.SpcGroup(key)[1]
2483        cell = sgtbxlattinp.sgtbx8[key][0]
2484        Axis = spdict['SGUniq']
2485        system = spdict['SGSys']
2486
2487        g2list = GenHLaue(dmin,spdict,cell2A(cell))
2488        #if len(g2list) != len(sgtbxlattinp.sgtbx8[key][1]):
2489        #    print 'failed',key,':' ,len(g2list),'vs',len(sgtbxlattinp.sgtbx8[key][1])
2490        #    print 'GSAS-II:'
2491        #    for h,k,l,d in g2list: print '  ',(h,k,l),d
2492        #    print 'SGTBX:'
2493        #    for hkllist,dref in sgtbxlattinp.sgtbx8[key][1]: print '  ',hkllist,dref
2494        assert len(g2list) == len(sgtbxlattinp.sgtbx8[key][1]), (
2495            'Reflection lists differ for %s' % key
2496            )
2497        #match = True
2498        for h,k,l,d in g2list:
2499            for hkllist,dref in sgtbxlattinp.sgtbx8[key][1]:
2500                if abs(d-dref) < derror:
2501                    if indexmatch((h,k,l,), hkllist, system, Axis): break
2502            else:
2503                assert 0,'No match for %s at %s (%s)' % ((h,k,l),d,key)
2504                #match = False
2505        #if not match:
2506            #for hkllist,dref in sgtbxlattinp.sgtbx8[key][1]: print '  ',hkllist,dref
2507            #print center, Laue, Axis, system
2508if __name__ == '__main__': selftestlist.append(test8)
2509
2510def test9():
2511    'test GenHLaue'
2512    _ReportTest()
2513    import GSASIIspc as G2spc
2514    if NeedTestData: TestData()
2515    for spc in LaueTestData:
2516        data = LaueTestData[spc]
2517        cell = data[0]
2518        hklm = np.array(data[1])
2519        H = hklm[-1][:3]
2520        hklO = hklm.T[:3].T
2521        A = cell2A(cell)
2522        dmin = 1./np.sqrt(calc_rDsq(H,A))
2523        SGData = G2spc.SpcGroup(spc)[1]
2524        hkls = np.array(GenHLaue(dmin,SGData,A))
2525        hklN = hkls.T[:3].T
2526        #print spc,hklO.shape,hklN.shape
2527        err = True
2528        for H in hklO:
2529            if H not in hklN:
2530                print ('%d %s'%(H,' missing from hkl from GSASII'))
2531                err = False
2532        assert(err)
2533if __name__ == '__main__': selftestlist.append(test9)
2534
2535
2536
2537
2538if __name__ == '__main__':
2539    # run self-tests
2540    selftestquiet = False
2541    for test in selftestlist:
2542        test()
2543    print ("OK")
Note: See TracBrowser for help on using the repository browser.