source: trunk/GSASIIlattice.py @ 4519

Last change on this file since 4519 was 4519, checked in by vondreele, 3 years ago

G2fpaGUI - minor cleanup
G2lattice - test on 'T' instead of 'C' so that 'B' PWDR types properly handled
G2math - add getPinkalpha, getPinkbeta % deriv routines for pink beam peak shapes

modify setPeakparams to use them

G2plot - test on 'T' instead of 'C' so that 'B' PWDR types properly handled
G2pwd - add pink beam peak shape function - same as TOF peak shape; scale sig & gam to be in centidegrees

  • add wtFactor to FitPeaks? fix bad PWDR weights for peak fitting so reduced chi2 is more appropriate

G2pwdGUI - test on 'T' instead of 'C' so that 'B' PWDR types properly handled

  • add wtFactor to DoPeakFit? to fix bad PWDR weights for peak fitting so reduced chi2 is more appropriate

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