source: trunk/GSASIIlattice.py @ 4195

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

add a find utility to GSASIIfiles for finding a file within a directory tree
faster version of ExpandCell? & define SwapItems? in GSASIIlattice
changes to RMCProfile interface

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