source: trunk/GSASIIlattice.py @ 4248

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

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