source: trunk/GSASIIlattice.py @ 3783

Last change on this file since 3783 was 3783, checked in by vondreele, 4 years ago

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