source: trunk/GSASIIlattice.py @ 3686

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

install new use of Bilbao site - call SUBGROUPS for making subgroup phase selections from parent structure
some modifications to mag subgroup stuff to accommodate

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