source: trunk/GSASIIlattice.py @ 3774

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

fix super indexing problem in transposeHKLF
fix reflection generation for incommensurate mag case in G2lattice & G2pwd
clean up non Fourier modulation calcs & remove analytic derivative stuff (now numeric)
fix uij derivative bug
work on incommensurate magnetic sturcture factors - not working yet
clean up of testDeriv - better choices for delt & force reflection regeneration before each derivative test

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