source: trunk/GSASIIlattice.py @ 4268

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

fixes to RMCProfile - GSAS-II stuff.
works with x-ray F(Q), etc.

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